← Back to branch summary

~ubuntu-branches/ubuntu/lucid/python-scipy/lucid

« back to all changes in this revision

Viewing changes to Lib/sandbox/pysparse/examples/poisson_test/mmio.c

  • Committer: Bazaar Package Importer
  • Author(s): Matthias Klose
  • Date: 2007-01-07 14:12:12 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20070107141212-mm0ebkh5b37hcpzn
http://bugs.debian.org/402783
* Remove build dependency on python-numpy-dev.
* python-scipy: Depend on python-numpy instead of python-numpy-dev.
* Package builds on other archs than i386. Closes: #402783.

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
}