← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/sparse/SuperLU/SRC/creadhb.c

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/sparse/SuperLU/SRC/creadhb.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 "csp_defs.h"
 
13
 
 
14
 
 
15
/* Eat up the rest of the current line */
 
16
int cDumpLine(FILE *fp)
 
17
{
 
18
    register int c;
 
19
    while ((c = fgetc(fp)) != '\n') ;
 
20
    return 0;
 
21
}
 
22
 
 
23
int cParseIntFormat(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 cParseFloatFormat(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 cReadVector(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 cReadValues(FILE *fp, int n, complex *destination, int perline, int persize)
 
85
{
 
86
    register int i, j, k, s, pair;
 
87
    register float 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
creadhb(int *nrow, int *ncol, int *nonz,
 
119
        complex **nzval, int **rowind, int **colptr)
 
120
{
 
121
/* 
 
122
 * Purpose
 
123
 * =======
 
124
 * 
 
125
 * Read a 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
    cDumpLine(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
    cDumpLine(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
    cDumpLine(fp);
 
233
 
 
234
    /* Allocate storage for the three arrays ( nzval, rowind, colptr ) */
 
235
    callocateA(*ncol, *nonz, nzval, rowind, colptr);
 
236
 
 
237
    /* Line 4: format statement */
 
238
    fscanf(fp, "%16c", buf);
 
239
    cParseIntFormat(buf, &colnum, &colsize);
 
240
    fscanf(fp, "%16c", buf);
 
241
    cParseIntFormat(buf, &rownum, &rowsize);
 
242
    fscanf(fp, "%20c", buf);
 
243
    cParseFloatFormat(buf, &valnum, &valsize);
 
244
    fscanf(fp, "%20c", buf);
 
245
    cDumpLine(fp);
 
246
 
 
247
    /* Line 5: right-hand side */    
 
248
    if ( rhscrd ) cDumpLine(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
    cReadVector(fp, *ncol+1, *colptr, colnum, colsize);
 
258
    cReadVector(fp, *nonz, *rowind, rownum, rowsize);
 
259
    if ( numer_lines ) {
 
260
        cReadValues(fp, *nonz, *nzval, valnum, valsize);
 
261
    }
 
262
    
 
263
    fclose(fp);
 
264
 
 
265
}
 
266