← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/examples/poisson_test/mmio.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/sandbox/pysparse/examples/poisson_test/mmio.c
1
 
/* 
2
 
*   Matrix Market I/O library for ANSI C
3
 
*
4
 
*   See http://math.nist.gov/MatrixMarket for details.
5
 
*
6
 
*
7
 
*/
8
 
 
9
 
 
10
 
#include <stdio.h>
11
 
#include <string.h>
12
 
#include <stdlib.h>
13
 
#include <ctype.h>
14
 
 
15
 
#include "mmio.h"
16
 
 
17
 
int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
18
 
                double **val_, int **I_, int **J_)
19
 
{
20
 
    FILE *f;
21
 
    MM_typecode matcode;
22
 
    int M, N, nz;
23
 
    int i;
24
 
    double *val;
25
 
    int *I, *J;
26
 
 
27
 
    if ((f = fopen(fname, "r")) == NULL)
28
 
            return -1;
29
 
 
30
 
 
31
 
    if (mm_read_banner(f, &matcode) != 0)
32
 
    {
33
 
        printf("mm_read_unsymetric: Could not process Matrix Market banner ");
34
 
        printf(" in file [%s]\n", fname);
35
 
        return -1;
36
 
    }
37
 
 
38
 
 
39
 
 
40
 
    if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
41
 
            mm_is_sparse(matcode)))
42
 
    {
43
 
        fprintf(stderr, "Sorry, this application does not support ");
44
 
        fprintf(stderr, "Market Market type: [%s]\n",
45
 
                mm_typecode_to_str(matcode));
46
 
        return -1;
47
 
    }
48
 
 
49
 
    /* find out size of sparse matrix: M, N, nz .... */
50
 
 
51
 
    if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
52
 
    {
53
 
        fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
54
 
        return -1;
55
 
    }
56
 
 
57
 
    *M_ = M;
58
 
    *N_ = N;
59
 
    *nz_ = nz;
60
 
 
61
 
    /* reseve memory for matrices */
62
 
 
63
 
    I = (int *) malloc(nz * sizeof(int));
64
 
    J = (int *) malloc(nz * sizeof(int));
65
 
    val = (double *) malloc(nz * sizeof(double));
66
 
 
67
 
    *val_ = val;
68
 
    *I_ = I;
69
 
    *J_ = J;
70
 
 
71
 
    /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
72
 
    /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
73
 
    /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
74
 
 
75
 
    for (i=0; i<nz; i++)
76
 
    {
77
 
        fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
78
 
        I[i]--;  /* adjust from 1-based to 0-based */
79
 
        J[i]--;
80
 
    }
81
 
    fclose(f);
82
 
 
83
 
    return 0;
84
 
}
85
 
 
86
 
int mm_is_valid(MM_typecode matcode)
87
 
{
88
 
    if (!mm_is_matrix(matcode)) return 0;
89
 
    if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
90
 
    if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
91
 
    if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || 
92
 
                mm_is_skew(matcode))) return 0;
93
 
    return 1;
94
 
}
95
 
 
96
 
int mm_read_banner(FILE *f, MM_typecode *matcode)
97
 
{
98
 
    char line[MM_MAX_LINE_LENGTH];
99
 
    char banner[MM_MAX_TOKEN_LENGTH];
100
 
    char mtx[MM_MAX_TOKEN_LENGTH]; 
101
 
    char crd[MM_MAX_TOKEN_LENGTH];
102
 
    char data_type[MM_MAX_TOKEN_LENGTH];
103
 
    char storage_scheme[MM_MAX_TOKEN_LENGTH];
104
 
    char *p;
105
 
 
106
 
 
107
 
    mm_clear_typecode(matcode);  
108
 
 
109
 
    if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) 
110
 
        return MM_PREMATURE_EOF;
111
 
 
112
 
    if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, 
113
 
        storage_scheme) != 5)
114
 
        return MM_PREMATURE_EOF;
115
 
 
116
 
    for (p=mtx; *p!='\0'; *p=tolower(*p),p++);  /* convert to lower case */
117
 
    for (p=crd; *p!='\0'; *p=tolower(*p),p++);  
118
 
    for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
119
 
    for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
120
 
 
121
 
    /* check for banner */
122
 
    if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
123
 
        return MM_NO_HEADER;
124
 
 
125
 
    /* first field should be "mtx" */
126
 
    if (strcmp(mtx, MM_MTX_STR) != 0)
127
 
        return  MM_UNSUPPORTED_TYPE;
128
 
    mm_set_matrix(matcode);
129
 
 
130
 
 
131
 
    /* second field describes whether this is a sparse matrix (in coordinate
132
 
            storgae) or a dense array */
133
 
 
134
 
 
135
 
    if (strcmp(crd, MM_SPARSE_STR) == 0)
136
 
        mm_set_sparse(matcode);
137
 
    else
138
 
    if (strcmp(crd, MM_DENSE_STR) == 0)
139
 
            mm_set_dense(matcode);
140
 
    else
141
 
        return MM_UNSUPPORTED_TYPE;
142
 
    
143
 
 
144
 
    /* third field */
145
 
 
146
 
    if (strcmp(data_type, MM_REAL_STR) == 0)
147
 
        mm_set_real(matcode);
148
 
    else
149
 
    if (strcmp(data_type, MM_COMPLEX_STR) == 0)
150
 
        mm_set_complex(matcode);
151
 
    else
152
 
    if (strcmp(data_type, MM_PATTERN_STR) == 0)
153
 
        mm_set_pattern(matcode);
154
 
    else
155
 
    if (strcmp(data_type, MM_INT_STR) == 0)
156
 
        mm_set_integer(matcode);
157
 
    else
158
 
        return MM_UNSUPPORTED_TYPE;
159
 
    
160
 
 
161
 
    /* fourth field */
162
 
 
163
 
    if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
164
 
        mm_set_general(matcode);
165
 
    else
166
 
    if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
167
 
        mm_set_symmetric(matcode);
168
 
    else
169
 
    if (strcmp(storage_scheme, MM_HERM_STR) == 0)
170
 
        mm_set_hermitian(matcode);
171
 
    else
172
 
    if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
173
 
        mm_set_skew(matcode);
174
 
    else
175
 
        return MM_UNSUPPORTED_TYPE;
176
 
        
177
 
 
178
 
    return 0;
179
 
}
180
 
 
181
 
int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
182
 
{
183
 
    if (fprintf(f, "%d %d %d\n", M, N, nz) != 3)
184
 
        return MM_COULD_NOT_WRITE_FILE;
185
 
    else 
186
 
        return 0;
187
 
}
188
 
 
189
 
int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
190
 
{
191
 
    char line[MM_MAX_LINE_LENGTH];
192
 
    int num_items_read;
193
 
 
194
 
    /* set return null parameter values, in case we exit with errors */
195
 
    *M = *N = *nz = 0;
196
 
 
197
 
    /* now continue scanning until you reach the end-of-comments */
198
 
    do 
199
 
    {
200
 
        if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
201
 
            return MM_PREMATURE_EOF;
202
 
    }while (line[0] == '%');
203
 
 
204
 
    /* line[] is either blank or has M,N, nz */
205
 
    if (sscanf(line, "%d %d %d", M, N, nz) == 3)
206
 
        return 0;
207
 
        
208
 
    else
209
 
    do
210
 
    { 
211
 
        num_items_read = fscanf(f, "%d %d %d", M, N, nz); 
212
 
        if (num_items_read == EOF) return MM_PREMATURE_EOF;
213
 
    }
214
 
    while (num_items_read != 3);
215
 
 
216
 
    return 0;
217
 
}
218
 
 
219
 
 
220
 
int mm_read_mtx_array_size(FILE *f, int *M, int *N)
221
 
{
222
 
    char line[MM_MAX_LINE_LENGTH];
223
 
    int num_items_read;
224
 
    /* set return null parameter values, in case we exit with errors */
225
 
    *M = *N = 0;
226
 
        
227
 
    /* now continue scanning until you reach the end-of-comments */
228
 
    do 
229
 
    {
230
 
        if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
231
 
            return MM_PREMATURE_EOF;
232
 
    }while (line[0] == '%');
233
 
 
234
 
    /* line[] is either blank or has M,N, nz */
235
 
    if (sscanf(line, "%d %d", M, N) == 2)
236
 
        return 0;
237
 
        
238
 
    else /* we have a blank line */
239
 
    do
240
 
    { 
241
 
        num_items_read = fscanf(f, "%d %d", M, N); 
242
 
        if (num_items_read == EOF) return MM_PREMATURE_EOF;
243
 
    }
244
 
    while (num_items_read != 2);
245
 
 
246
 
    return 0;
247
 
}
248
 
 
249
 
int mm_write_mtx_array_size(FILE *f, int M, int N)
250
 
{
251
 
    if (fprintf(f, "%d %d\n", M, N) != 2)
252
 
        return MM_COULD_NOT_WRITE_FILE;
253
 
    else 
254
 
        return 0;
255
 
}
256
 
 
257
 
 
258
 
 
259
 
/*-------------------------------------------------------------------------*/
260
 
 
261
 
/******************************************************************/
262
 
/* use when I[], J[], and val[]J, and val[] are already allocated */
263
 
/******************************************************************/
264
 
 
265
 
int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
266
 
        double val[], MM_typecode matcode)
267
 
{
268
 
    int i;
269
 
    if (mm_is_complex(matcode))
270
 
    {
271
 
        for (i=0; i<nz; i++)
272
 
            if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
273
 
                != 4) return MM_PREMATURE_EOF;
274
 
    }
275
 
    else if (mm_is_real(matcode))
276
 
    {
277
 
        for (i=0; i<nz; i++)
278
 
        {
279
 
            if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
280
 
                != 3) return MM_PREMATURE_EOF;
281
 
 
282
 
        }
283
 
    }
284
 
 
285
 
    else if (mm_is_pattern(matcode))
286
 
    {
287
 
        for (i=0; i<nz; i++)
288
 
            if (fscanf(f, "%d %d", &I[i], &J[i])
289
 
                != 2) return MM_PREMATURE_EOF;
290
 
    }
291
 
    else
292
 
        return MM_UNSUPPORTED_TYPE;
293
 
 
294
 
    return 0;
295
 
        
296
 
}
297
 
 
298
 
int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
299
 
        double *real, double *imag, MM_typecode matcode)
300
 
{
301
 
    if (mm_is_complex(matcode))
302
 
    {
303
 
            if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
304
 
                != 4) return MM_PREMATURE_EOF;
305
 
    }
306
 
    else if (mm_is_real(matcode))
307
 
    {
308
 
            if (fscanf(f, "%d %d %lg\n", I, J, real)
309
 
                != 3) return MM_PREMATURE_EOF;
310
 
 
311
 
    }
312
 
 
313
 
    else if (mm_is_pattern(matcode))
314
 
    {
315
 
            if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
316
 
    }
317
 
    else
318
 
        return MM_UNSUPPORTED_TYPE;
319
 
 
320
 
    return 0;
321
 
        
322
 
}
323
 
 
324
 
 
325
 
/************************************************************************
326
 
    mm_read_mtx_crd()  fills M, N, nz, array of values, and return
327
 
                        type code, e.g. 'MCRS'
328
 
 
329
 
                        if matrix is complex, values[] is of size 2*nz,
330
 
                            (nz pairs of real/imaginary values)
331
 
************************************************************************/
332
 
 
333
 
int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, 
334
 
        double **val, MM_typecode *matcode)
335
 
{
336
 
    int ret_code;
337
 
    FILE *f;
338
 
 
339
 
    if (strcmp(fname, "stdin") == 0) f=stdin;
340
 
    else
341
 
    if ((f = fopen(fname, "r")) == NULL)
342
 
        return MM_COULD_NOT_READ_FILE;
343
 
 
344
 
 
345
 
    if ((ret_code = mm_read_banner(f, matcode)) != 0)
346
 
        return ret_code;
347
 
 
348
 
    if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && 
349
 
            mm_is_matrix(*matcode)))
350
 
        return MM_UNSUPPORTED_TYPE;
351
 
 
352
 
    if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
353
 
        return ret_code;
354
 
 
355
 
 
356
 
    *I = (int *)  malloc(*nz * sizeof(int));
357
 
    *J = (int *)  malloc(*nz * sizeof(int));
358
 
    *val = NULL;
359
 
 
360
 
    if (mm_is_complex(*matcode))
361
 
    {
362
 
        *val = (double *) malloc(*nz * 2 * sizeof(double));
363
 
        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
364
 
                *matcode);
365
 
        if (ret_code != 0) return ret_code;
366
 
    }
367
 
    else if (mm_is_real(*matcode))
368
 
    {
369
 
        *val = (double *) malloc(*nz * sizeof(double));
370
 
        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
371
 
                *matcode);
372
 
        if (ret_code != 0) return ret_code;
373
 
    }
374
 
 
375
 
    else if (mm_is_pattern(*matcode))
376
 
    {
377
 
        ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
378
 
                *matcode);
379
 
        if (ret_code != 0) return ret_code;
380
 
    }
381
 
 
382
 
    if (f != stdin) fclose(f);
383
 
    return 0;
384
 
}
385
 
 
386
 
int mm_write_banner(FILE *f, MM_typecode matcode)
387
 
{
388
 
    char *str = mm_typecode_to_str(matcode);
389
 
    int ret_code;
390
 
 
391
 
    ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, str);
392
 
    free(str);
393
 
    if (ret_code !=2 )
394
 
        return MM_COULD_NOT_WRITE_FILE;
395
 
    else
396
 
        return 0;
397
 
}
398
 
 
399
 
int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
400
 
        double val[], MM_typecode matcode)
401
 
{
402
 
    FILE *f;
403
 
    int i;
404
 
 
405
 
    if (strcmp(fname, "stdout") == 0) 
406
 
        f = stdout;
407
 
    else
408
 
    if ((f = fopen(fname, "w")) == NULL)
409
 
        return MM_COULD_NOT_WRITE_FILE;
410
 
    
411
 
    /* print banner followed by typecode */
412
 
    fprintf(f, "%s ", MatrixMarketBanner);
413
 
    fprintf(f, "%s\n", mm_typecode_to_str(matcode));
414
 
 
415
 
    /* print matrix sizes and nonzeros */
416
 
    fprintf(f, "%d %d %d\n", M, N, nz);
417
 
 
418
 
    /* print values */
419
 
    if (mm_is_pattern(matcode))
420
 
        for (i=0; i<nz; i++)
421
 
            fprintf(f, "%d %d\n", I[i], J[i]);
422
 
    else
423
 
    if (mm_is_real(matcode))
424
 
        for (i=0; i<nz; i++)
425
 
            fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
426
 
    else
427
 
    if (mm_is_complex(matcode))
428
 
        for (i=0; i<nz; i++)
429
 
            fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i], 
430
 
                        val[2*i+1]);
431
 
    else
432
 
    {
433
 
        if (f != stdout) fclose(f);
434
 
        return MM_UNSUPPORTED_TYPE;
435
 
    }
436
 
 
437
 
    if (f !=stdout) fclose(f);
438
 
 
439
 
    return 0;
440
 
}
441
 
    
442
 
 
443
 
char  *mm_typecode_to_str(MM_typecode matcode)
444
 
{
445
 
    char buffer[MM_MAX_LINE_LENGTH];
446
 
    char *types[4];
447
 
    int error =0;
448
 
 
449
 
    /* check for MTX type */
450
 
    if (mm_is_matrix(matcode)) 
451
 
        types[0] = MM_MTX_STR;
452
 
    else
453
 
        error=1;
454
 
 
455
 
    /* check for CRD or ARR matrix */
456
 
    if (mm_is_sparse(matcode))
457
 
        types[1] = MM_SPARSE_STR;
458
 
    else
459
 
    if (mm_is_dense(matcode))
460
 
        types[1] = MM_DENSE_STR;
461
 
    else
462
 
        return NULL;
463
 
 
464
 
    /* check for element data type */
465
 
    if (mm_is_real(matcode))
466
 
        types[2] = MM_REAL_STR;
467
 
    else
468
 
    if (mm_is_complex(matcode))
469
 
        types[2] = MM_COMPLEX_STR;
470
 
    else
471
 
    if (mm_is_pattern(matcode))
472
 
        types[2] = MM_PATTERN_STR;
473
 
    else
474
 
    if (mm_is_integer(matcode))
475
 
        types[2] = MM_INT_STR;
476
 
    else
477
 
        return NULL;
478
 
 
479
 
 
480
 
    /* check for symmetry type */
481
 
    if (mm_is_general(matcode))
482
 
        types[3] = MM_GENERAL_STR;
483
 
    else
484
 
    if (mm_is_symmetric(matcode))
485
 
        types[3] = MM_SYMM_STR;
486
 
    else 
487
 
    if (mm_is_hermitian(matcode))
488
 
        types[3] = MM_HERM_STR;
489
 
    else 
490
 
    if (mm_is_skew(matcode))
491
 
        types[3] = MM_SKEW_STR;
492
 
    else
493
 
        return NULL;
494
 
 
495
 
    sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
496
 
    return strdup(buffer);
497
 
 
498
 
}