Mercurial > forge

comparison main/sparse/SuperLU/SRC/zreadhb.c @ 0:6b33357c7561 octave-forge

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