1
2 #include <stdio.h>
3 #include <string.h>
4 #include "mmio.h"
5 #include "hbio.h"
6 #include "lusolio.h"
7 #include <R.h>
8
9 /* Utility routines to read matrix files in the Coordinate Text File format*/
10
ctf_read_A(char * filename,int maxm,int maxn,int maxnz,int * m,int * n,int * nnzero,int * iA,int * jA,REAL * Aij)11 MYBOOL ctf_read_A(char *filename, int maxm, int maxn, int maxnz,
12 int *m, int *n, int *nnzero, int *iA, int *jA, REAL *Aij)
13 {
14 FILE *iofile;
15 int eof;
16 char buffer[100];
17 int k, i, j;
18 REAL Ak;
19 MYBOOL filldata;
20
21 *nnzero = 0;
22 *m = 0;
23 *n = 0;
24
25 iofile = fopen(filename, "r" );
26 if(iofile == NULL) {
27 Rprintf("A file %s does not exist\n", filename);
28 return( FALSE );
29 }
30
31 filldata = (MYBOOL) !((iA == NULL) && (jA == NULL) && (Aij == NULL));
32 eof = TRUE;
33 for (k = 1; k <= maxnz; k++) {
34 eof = feof(iofile);
35 if(eof)
36 break;
37 eof = fscanf(iofile, "%d %d %s", &i, &j, buffer);
38 if(eof == 0 || eof == EOF || i <= 0 || j <= 0 || strlen(buffer) == 0)
39 break;
40 Ak = atof(buffer);
41 (*nnzero)++;
42 if (filldata) {
43 iA[k] = i;
44 jA[k] = j;
45 Aij[k] = Ak;
46 }
47 if (i > *m) *m = i;
48 if (j > *n) *n = j;
49 }
50 fclose( iofile );
51 if(!eof) {
52 Rprintf("Too much data in A file. Increase maxnz\n");
53 Rprintf("Current maxnz = %d\n", maxnz);
54 return( FALSE );
55 }
56 Rprintf("A read successfully\n");
57 Rprintf("m = %d n = %d nnzero = %d\n",
58 *m, *n, *nnzero);
59 if (*m > maxm || *n > maxn) {
60 Rprintf("However, matrix dimensions exceed maxm or maxn\n");
61 return( FALSE );
62 }
63 return( TRUE );
64 }
65
ctf_size_A(char * filename,int * m,int * n,int * nnzero)66 MYBOOL ctf_size_A(char *filename, int *m, int *n, int *nnzero)
67 {
68 int maxint = 200000000;
69
70 return( ctf_read_A(filename, maxint, maxint, maxint,
71 m, n, nnzero, NULL, NULL, NULL) );
72 }
73
ctf_read_b(char * filename,int m,REAL * b)74 MYBOOL ctf_read_b(char *filename, int m, REAL *b)
75 {
76 FILE *iofile;
77 int eof;
78 char buffer[100];
79 int i;
80
81 iofile = fopen(filename, "r");
82 if(iofile == NULL) {
83 Rprintf("b file %s does not exist\n", filename);
84 return( FALSE );
85 }
86
87 for (i = 1; i <= m; i++) {
88 if(feof(iofile))
89 goto x350;
90 eof = fscanf(iofile, "%s", buffer);
91 if(eof == 0 || eof == EOF)
92 goto x350;
93 b[i] = atof(buffer);
94 }
95
96 fclose( iofile );
97 Rprintf("b read successfully\n");
98 return( TRUE );
99
100 x350:
101 fclose( iofile );
102 Rprintf("Not enough data in b file.\n");
103 return( FALSE );
104 }
105
106
107 /* Utility routines to read matrix files in the MatrixMarket format*/
108 #define mmf_recsize 255
mmf_read_A(char * filename,int maxM,int maxN,int maxnz,int * M,int * N,int * nz,int * iA,int * jA,REAL * Aij)109 MYBOOL mmf_read_A(char *filename, int maxM, int maxN, int maxnz,
110 int *M, int *N, int *nz, int *iA, int *jA, REAL *Aij)
111 {
112 MM_typecode matcode;
113 FILE *f;
114 int i, k, ret_code, ival, jval;
115 REAL Aval;
116 MYBOOL status = FALSE, ispat, filldata;
117 char buf[mmf_recsize];
118
119 f = fopen(filename, "r");
120 if(f == NULL)
121 return( status );
122
123 if(mm_read_banner(f, &matcode) != 0) {
124 Rprintf("Could not process Matrix Market banner.\n");
125 goto x900;
126 }
127
128 /* Screen matrix types since LUSOL only supports a
129 subset of the Matrix Market data types. */
130 if(mm_is_complex(matcode) || mm_is_pattern(matcode)) {
131 Rprintf("Sorry, this application does not support ");
132 Rprintf("Market Market type: [%s]\n", mm_typecode_to_str(matcode));
133 goto x900;
134 }
135
136 /* Verify that we have sufficient array storage */
137 filldata = (MYBOOL) !((iA == NULL) && (jA == NULL) && (Aij == NULL));
138 if(filldata && maxN > 1 && jA == NULL) {
139 Rprintf("Market Market insufficient array storage specified\n");
140 goto x900;
141 }
142
143 /* Find out size of sparse matrix .... */
144 ret_code = mm_read_mtx_crd_size(f, M, N, nz);
145 if(ret_code != 0 || !filldata || (*M > maxM) || (*N > maxN) || (*nz > maxnz)) {
146 status = !filldata;
147 goto x900;
148 }
149
150
151 /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
152 /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
153 /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
154
155 /* Read dense matrix in column order */
156 ispat = (MYBOOL) mm_is_pattern(matcode);
157 k = 1;
158 if(mm_is_dense(matcode)) {
159 maxN = MIN(maxN, *N);
160 for (jval = 1; jval <= maxN; jval++) {
161 for (i = 1; i <= *M; i++) {
162 if(fgets(buf, mmf_recsize-1, f) == NULL)
163 break;
164 if(sscanf(buf, "%lg\n", &Aval) == 0)
165 break;
166 if(Aval != 0) {
167 if(iA != NULL)
168 iA[k] = i;
169 if(jA != NULL)
170 jA[k] = jval;
171
172 /* Make sure we handle dense vector reading properly */
173 if(iA == NULL && jA == NULL)
174 Aij[i] = Aval;
175 else
176 Aij[k] = Aval;
177 k++;
178 }
179 }
180 }
181 }
182 /* Read sparse matrix by coordinate */
183 else {
184 for (i = 1; i <= *nz; i++) {
185 if(fgets(buf, mmf_recsize-1, f) == NULL)
186 break;
187 if(buf[0] == '%')
188 continue;
189 if(ispat) {
190 if(sscanf(buf, "%d %d\n", &ival, &jval) == 0)
191 continue;
192 Aij[k] = 1.0;
193 }
194 else
195 if(sscanf(buf, "%d %d %lg\n", &ival, &jval, &Aval) == 0)
196 continue;
197
198 /* Check if it is a nonzero and we are within column dimension */
199 if(Aval != 0 && jval <= maxN) {
200 Aij[k] = Aval;
201 if(iA != NULL)
202 iA[k] = ival;
203 if(jA != NULL)
204 jA[k] = jval;
205 k++;
206 }
207 }
208 }
209 *nz = k - 1;
210
211 /* Handle case where only the lower triangular parts are given */
212 if(!mm_is_general(matcode)) {
213 if((M != N) || (maxN != maxM) || (2*(*nz) > maxnz)) {
214 Rprintf("Market Market cannot fill in symmetry data\n");
215 goto x900;
216 }
217 ispat = mm_is_skew(matcode);
218 for(i = 1; i <= *nz; i++) {
219 iA[k] = jA[i];
220 jA[k] = iA[i];
221 if(ispat)
222 Aij[k] = -Aij[i];
223 else
224 Aij[k] = Aij[i];
225 k++;
226 }
227 *nz = k - 1;
228 }
229 status = TRUE;
230
231 /* Finish up */
232 x900:
233 fclose( f );
234 return( status );
235 }
236
mmf_size_A(char * filename,int * M,int * N,int * nz)237 MYBOOL mmf_size_A(char *filename, int *M, int *N, int *nz)
238 {
239 int maxint = 200000000;
240
241 return( mmf_read_A(filename, maxint, maxint, maxint,
242 M, N, nz, NULL, NULL, NULL) );
243 }
244
mmf_read_b(char * filename,int m,REAL * b)245 MYBOOL mmf_read_b(char *filename, int m, REAL *b)
246 {
247 int im, jn, nnzero;
248
249 return( mmf_read_A(filename, m, 1, m,
250 &im, &jn, &nnzero, NULL, NULL, b));
251 }
252
253
254 /* Utility routines to read matrix files in Harwell-Boeing format*/
255
hbf_read_A(char * filename,int maxM,int maxN,int maxnz,int * M,int * N,int * nz,int * iA,int * jA,REAL * Aij)256 MYBOOL hbf_read_A(char *filename, int maxM, int maxN, int maxnz,
257 int *M, int *N, int *nz, int *iA, int *jA, REAL *Aij)
258 {
259 MYBOOL success;
260
261 success = hbf_size_A(filename, M, N, nz);
262 if(!success)
263 return( success );
264
265 Aij[1] = 0;
266 success = (MYBOOL) readHB_mat_double(filename, jA, iA-1, Aij-1);
267
268 /* Test if we have a pattern matrix and fill it with all zeros */
269 if(Aij[1] == 0) {
270 int i;
271 for(i = 1; i <= *nz; i++)
272 Aij[i] = 1;
273 }
274
275 /* Expand the column nz counts to triplet format */
276 if(success) {
277 int i, j, ii, k;
278 k = *nz;
279 for(j = *N; j > 0; j--) {
280 ii = jA[j];
281 for(i = jA[j-1]; i < ii; i++, k--)
282 jA[k] = j;
283 }
284 }
285 return( success );
286 }
287
hbf_size_A(char * filename,int * M,int * N,int * nz)288 MYBOOL hbf_size_A(char *filename, int *M, int *N, int *nz)
289 {
290 int Nrhs;
291 char *Type;
292
293 return( (MYBOOL) readHB_info(filename, M, N, nz, &Type, &Nrhs) );
294 }
295
hbf_read_b(char * filename,int m,REAL * b)296 MYBOOL hbf_read_b(char *filename, int m, REAL *b)
297 {
298 return( (MYBOOL) readHB_aux_double(filename, 'F', b) ); /* Same format as matrix */
299 }
300