← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to Lib/linsolve/SuperLU/SRC/zreadhb.c

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/linsolve/SuperLU/SRC/zreadhb.c
1
 
 
2
 
 
3
 
/*
4
 
 * -- SuperLU routine (version 2.0) --
5
 
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
6
 
 * and Lawrence Berkeley National Lab.
7
 
 * November 15, 1997
8
 
 *
9
 
 */
10
 
#include <stdio.h>
11
 
#include <stdlib.h>
12
 
#include "zsp_defs.h"
13
 
 
14
 
 
15
 
/* Eat up the rest of the current line */
16
 
int zDumpLine(FILE *fp)
17
 
{
18
 
    register int c;
19
 
    while ((c = fgetc(fp)) != '\n') ;
20
 
    return 0;
21
 
}
22
 
 
23
 
int zParseIntFormat(char *buf, int *num, int *size)
24
 
{
25
 
    char *tmp;
26
 
 
27
 
    tmp = buf;
28
 
    while (*tmp++ != '(') ;
29
 
    sscanf(tmp, "%d", num);
30
 
    while (*tmp != 'I' && *tmp != 'i') ++tmp;
31
 
    ++tmp;
32
 
    sscanf(tmp, "%d", size);
33
 
    return 0;
34
 
}
35
 
 
36
 
int zParseFloatFormat(char *buf, int *num, int *size)
37
 
{
38
 
    char *tmp, *period;
39
 
    
40
 
    tmp = buf;
41
 
    while (*tmp++ != '(') ;
42
 
    *num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
43
 
    while (*tmp != 'E' && *tmp != 'e' && *tmp != 'D' && *tmp != 'd'
44
 
           && *tmp != 'F' && *tmp != 'f') {
45
 
        /* May find kP before nE/nD/nF, like (1P6F13.6). In this case the
46
 
           num picked up refers to P, which should be skipped. */
47
 
        if (*tmp=='p' || *tmp=='P') {
48
 
           ++tmp;
49
 
           *num = atoi(tmp); /*sscanf(tmp, "%d", num);*/
50
 
        } else {
51
 
           ++tmp;
52
 
        }
53
 
    }
54
 
    ++tmp;
55
 
    period = tmp;
56
 
    while (*period != '.' && *period != ')') ++period ;
57
 
    *period = '\0';
58
 
    *size = atoi(tmp); /*sscanf(tmp, "%2d", size);*/
59
 
 
60
 
    return 0;
61
 
}
62
 
 
63
 
int zReadVector(FILE *fp, int n, int *where, int perline, int persize)
64
 
{
65
 
    register int i, j, item;
66
 
    char tmp, buf[100];
67
 
    
68
 
    i = 0;
69
 
    while (i < n) {
70
 
        fgets(buf, 100, fp);    /* read a line at a time */
71
 
        for (j=0; j<perline && i<n; j++) {
72
 
            tmp = buf[(j+1)*persize];     /* save the char at that place */
73
 
            buf[(j+1)*persize] = 0;       /* null terminate */
74
 
            item = atoi(&buf[j*persize]); 
75
 
            buf[(j+1)*persize] = tmp;     /* recover the char at that place */
76
 
            where[i++] = item - 1;
77
 
        }
78
 
    }
79
 
 
80
 
    return 0;
81
 
}
82
 
 
83
 
/* Read complex numbers as pairs of (real, imaginary) */
84
 
int zReadValues(FILE *fp, int n, doublecomplex *destination, int perline, int persize)
85
 
{
86
 
    register int i, j, k, s, pair;
87
 
    register double realpart;
88
 
    char tmp, buf[100];
89
 
    
90
 
    i = pair = 0;
91
 
    while (i < n) {
92
 
        fgets(buf, 100, fp);    /* read a line at a time */
93
 
        for (j=0; j<perline && i<n; j++) {
94
 
            tmp = buf[(j+1)*persize];     /* save the char at that place */
95
 
            buf[(j+1)*persize] = 0;       /* null terminate */
96
 
            s = j*persize;
97
 
            for (k = 0; k < persize; ++k) /* No D_ format in C */
98
 
                if ( buf[s+k] == 'D' || buf[s+k] == 'd' ) buf[s+k] = 'E';
99
 
            if ( pair == 0 ) {
100
 
                /* The value is real part */
101
 
                realpart = atof(&buf[s]);
102
 
                pair = 1;
103
 
            } else {
104
 
                /* The value is imaginary part */
105
 
                destination[i].r = realpart;
106
 
                destination[i++].i = atof(&buf[s]);
107
 
                pair = 0;
108
 
            }
109
 
            buf[(j+1)*persize] = tmp;     /* recover the char at that place */
110
 
        }
111
 
    }
112
 
 
113
 
    return 0;
114
 
}
115
 
 
116
 
 
117
 
void
118
 
zreadhb(int *nrow, int *ncol, int *nonz,
119
 
        doublecomplex **nzval, int **rowind, int **colptr)
120
 
{
121
 
/* 
122
 
 * Purpose
123
 
 * =======
124
 
 * 
125
 
 * Read a DOUBLE COMPLEX PRECISION matrix stored in Harwell-Boeing format 
126
 
 * as described below.
127
 
 * 
128
 
 * Line 1 (A72,A8) 
129
 
 *      Col. 1 - 72   Title (TITLE) 
130
 
 *      Col. 73 - 80  Key (KEY) 
131
 
 * 
132
 
 * Line 2 (5I14) 
133
 
 *      Col. 1 - 14   Total number of lines excluding header (TOTCRD) 
134
 
 *      Col. 15 - 28  Number of lines for pointers (PTRCRD) 
135
 
 *      Col. 29 - 42  Number of lines for row (or variable) indices (INDCRD) 
136
 
 *      Col. 43 - 56  Number of lines for numerical values (VALCRD) 
137
 
 *      Col. 57 - 70  Number of lines for right-hand sides (RHSCRD) 
138
 
 *                    (including starting guesses and solution vectors 
139
 
 *                     if present) 
140
 
 *                    (zero indicates no right-hand side data is present) 
141
 
 *
142
 
 * Line 3 (A3, 11X, 4I14) 
143
 
 *      Col. 1 - 3    Matrix type (see below) (MXTYPE) 
144
 
 *      Col. 15 - 28  Number of rows (or variables) (NROW) 
145
 
 *      Col. 29 - 42  Number of columns (or elements) (NCOL) 
146
 
 *      Col. 43 - 56  Number of row (or variable) indices (NNZERO) 
147
 
 *                    (equal to number of entries for assembled matrices) 
148
 
 *      Col. 57 - 70  Number of elemental matrix entries (NELTVL) 
149
 
 *                    (zero in the case of assembled matrices) 
150
 
 * Line 4 (2A16, 2A20) 
151
 
 *      Col. 1 - 16   Format for pointers (PTRFMT) 
152
 
 *      Col. 17 - 32  Format for row (or variable) indices (INDFMT) 
153
 
 *      Col. 33 - 52  Format for numerical values of coefficient matrix (VALFMT) 
154
 
 *      Col. 53 - 72 Format for numerical values of right-hand sides (RHSFMT) 
155
 
 *
156
 
 * Line 5 (A3, 11X, 2I14) Only present if there are right-hand sides present 
157
 
 *      Col. 1        Right-hand side type: 
158
 
 *                        F for full storage or M for same format as matrix 
159
 
 *      Col. 2        G if a starting vector(s) (Guess) is supplied. (RHSTYP) 
160
 
 *      Col. 3        X if an exact solution vector(s) is supplied. 
161
 
 *      Col. 15 - 28  Number of right-hand sides (NRHS) 
162
 
 *      Col. 29 - 42  Number of row indices (NRHSIX) 
163
 
 *                    (ignored in case of unassembled matrices) 
164
 
 *
165
 
 * The three character type field on line 3 describes the matrix type. 
166
 
 * The following table lists the permitted values for each of the three 
167
 
 * characters. As an example of the type field, RSA denotes that the matrix 
168
 
 * is real, symmetric, and assembled. 
169
 
 *
170
 
 * First Character: 
171
 
 *      R Real matrix 
172
 
 *      C Complex matrix 
173
 
 *      P Pattern only (no numerical values supplied) 
174
 
 *
175
 
 * Second Character: 
176
 
 *      S Symmetric 
177
 
 *      U Unsymmetric 
178
 
 *      H Hermitian 
179
 
 *      Z Skew symmetric 
180
 
 *      R Rectangular 
181
 
 *
182
 
 * Third Character: 
183
 
 *      A Assembled 
184
 
 *      E Elemental matrices (unassembled) 
185
 
 *
186
 
 */
187
 
 
188
 
    register int i, numer_lines = 0, rhscrd = 0;
189
 
    int tmp, colnum, colsize, rownum, rowsize, valnum, valsize;
190
 
    char buf[100], type[4], key[10];
191
 
    FILE *fp;
192
 
 
193
 
    fp = stdin;
194
 
 
195
 
    /* Line 1 */
196
 
    fgets(buf, 100, fp);
197
 
    fputs(buf, stdout);
198
 
#if 0
199
 
    fscanf(fp, "%72c", buf); buf[72] = 0;
200
 
    printf("Title: %s", buf);
201
 
    fscanf(fp, "%8c", key);  key[8] = 0;
202
 
    printf("Key: %s\n", key);
203
 
    zDumpLine(fp);
204
 
#endif
205
 
 
206
 
    /* Line 2 */
207
 
    for (i=0; i<5; i++) {
208
 
        fscanf(fp, "%14c", buf); buf[14] = 0;
209
 
        sscanf(buf, "%d", &tmp);
210
 
        if (i == 3) numer_lines = tmp;
211
 
        if (i == 4 && tmp) rhscrd = tmp;
212
 
    }
213
 
    zDumpLine(fp);
214
 
 
215
 
    /* Line 3 */
216
 
    fscanf(fp, "%3c", type);
217
 
    fscanf(fp, "%11c", buf); /* pad */
218
 
    type[3] = 0;
219
 
#ifdef DEBUG
220
 
    printf("Matrix type %s\n", type);
221
 
#endif
222
 
    
223
 
    fscanf(fp, "%14c", buf); sscanf(buf, "%d", nrow);
224
 
    fscanf(fp, "%14c", buf); sscanf(buf, "%d", ncol);
225
 
    fscanf(fp, "%14c", buf); sscanf(buf, "%d", nonz);
226
 
    fscanf(fp, "%14c", buf); sscanf(buf, "%d", &tmp);
227
 
    
228
 
    if (tmp != 0)
229
 
          printf("This is not an assembled matrix!\n");
230
 
    if (*nrow != *ncol)
231
 
        printf("Matrix is not square.\n");
232
 
    zDumpLine(fp);
233
 
 
234
 
    /* Allocate storage for the three arrays ( nzval, rowind, colptr ) */
235
 
    zallocateA(*ncol, *nonz, nzval, rowind, colptr);
236
 
 
237
 
    /* Line 4: format statement */
238
 
    fscanf(fp, "%16c", buf);
239
 
    zParseIntFormat(buf, &colnum, &colsize);
240
 
    fscanf(fp, "%16c", buf);
241
 
    zParseIntFormat(buf, &rownum, &rowsize);
242
 
    fscanf(fp, "%20c", buf);
243
 
    zParseFloatFormat(buf, &valnum, &valsize);
244
 
    fscanf(fp, "%20c", buf);
245
 
    zDumpLine(fp);
246
 
 
247
 
    /* Line 5: right-hand side */    
248
 
    if ( rhscrd ) zDumpLine(fp); /* skip RHSFMT */
249
 
    
250
 
#ifdef DEBUG
251
 
    printf("%d rows, %d nonzeros\n", *nrow, *nonz);
252
 
    printf("colnum %d, colsize %d\n", colnum, colsize);
253
 
    printf("rownum %d, rowsize %d\n", rownum, rowsize);
254
 
    printf("valnum %d, valsize %d\n", valnum, valsize);
255
 
#endif
256
 
    
257
 
    zReadVector(fp, *ncol+1, *colptr, colnum, colsize);
258
 
    zReadVector(fp, *nonz, *rowind, rownum, rowsize);
259
 
    if ( numer_lines ) {
260
 
        zReadValues(fp, *nonz, *nzval, valnum, valsize);
261
 
    }
262
 
    
263
 
    fclose(fp);
264
 
 
265
 
}
266