EpetraExt Development
EpetraExt_mmio.cpp
Go to the documentation of this file.
00001 /* 
00002 *   Matrix Market I/O library for ANSI C
00003 *
00004 *   See http://math.nist.gov/MatrixMarket for details.
00005 *
00006 *
00007 */
00008 
00009 /* JW - This is essentially a C file that was "converted" to C++, but still
00010    contains C function calls.  Since it is independent of the rest of
00011    the package, we will not include the central ConfigDefs file, but
00012    will rather #include the C headers that we need.  This should fix the
00013    build problem on sass2889.
00014 #include <Epetra_ConfigDefs.h> */
00015 #include <string.h>
00016 #include <stdio.h>
00017 #include <ctype.h>
00018 #include "EpetraExt_mmio.h"
00019 
00020 using namespace EpetraExt;
00021 namespace EpetraExt {
00022 int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
00023                 double **val_, int **I_, int **J_)
00024 {
00025     FILE *f;
00026     MM_typecode matcode;
00027     int M, N, nz;
00028     int i;
00029     double *val;
00030     int *I, *J;
00031  
00032     if ((f = fopen(fname, "r")) == NULL)
00033             return -1;
00034  
00035  
00036     if (mm_read_banner(f, &matcode) != 0)
00037     {
00038         printf("mm_read_unsymetric: Could not process Matrix Market banner ");
00039         printf(" in file [%s]\n", fname);
00040         return -1;
00041     }
00042  
00043  
00044  
00045     if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
00046             mm_is_sparse(matcode)))
00047     {
00048       char buffer[MM_MAX_LINE_LENGTH];
00049       mm_typecode_to_str(matcode, buffer);
00050       fprintf(stderr, "Sorry, this application does not support ");
00051       fprintf(stderr, "Market Market type: [%s]\n",buffer);
00052         return -1;
00053     }
00054  
00055     /* find out size of sparse matrix: M, N, nz .... */
00056  
00057     if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
00058     {
00059         fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
00060         return -1;
00061     }
00062  
00063     *M_ = M;
00064     *N_ = N;
00065     *nz_ = nz;
00066  
00067     /* reseve memory for matrices */
00068  
00069     //I = (int *) malloc(nz * sizeof(int));
00070     //J = (int *) malloc(nz * sizeof(int));
00071     //val = (double *) malloc(nz * sizeof(double));
00072  
00073     I = new int[nz];
00074     J = new int[nz];
00075     val = new double[nz];
00076  
00077     *val_ = val;
00078     *I_ = I;
00079     *J_ = J;
00080  
00081     /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
00082     /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
00083     /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
00084  
00085     for (i=0; i<nz; i++)
00086     {
00087         fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
00088         I[i]--;  /* adjust from 1-based to 0-based */
00089         J[i]--;
00090     }
00091     fclose(f);
00092  
00093     return 0;
00094 }
00095 
00096 int mm_is_valid(MM_typecode matcode)
00097 {
00098     if (!mm_is_matrix(matcode)) return 0;
00099     if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
00100     if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
00101     if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || 
00102                 mm_is_skew(matcode))) return 0;
00103     return 1;
00104 }
00105 
00106 int mm_read_banner(FILE *f, MM_typecode *matcode)
00107 {
00108     char line[MM_MAX_LINE_LENGTH];
00109     char banner[MM_MAX_TOKEN_LENGTH];
00110     char mtx[MM_MAX_TOKEN_LENGTH]; 
00111     char crd[MM_MAX_TOKEN_LENGTH];
00112     char data_type[MM_MAX_TOKEN_LENGTH];
00113     char storage_scheme[MM_MAX_TOKEN_LENGTH];
00114     char *p;
00115 
00116 
00117     mm_clear_typecode(matcode);  
00118 
00119     if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) 
00120         return MM_PREMATURE_EOF;
00121 
00122     if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, 
00123         storage_scheme) != 5)
00124         return MM_PREMATURE_EOF;
00125 
00126     for (p=mtx; *p!='\0'; *p=tolower(*p),p++);  /* convert to lower case */
00127     for (p=crd; *p!='\0'; *p=tolower(*p),p++);  
00128     for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
00129     for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
00130 
00131     /* check for banner */
00132     if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
00133         return MM_NO_HEADER;
00134 
00135     /* first field should be "mtx" */
00136     if (strcmp(mtx, MM_MTX_STR) != 0)
00137         return  MM_UNSUPPORTED_TYPE;
00138     mm_set_matrix(matcode);
00139 
00140 
00141     /* second field describes whether this is a sparse matrix (in coordinate
00142             storgae) or a dense array */
00143 
00144 
00145     if (strcmp(crd, MM_SPARSE_STR) == 0)
00146         mm_set_sparse(matcode);
00147     else
00148     if (strcmp(crd, MM_DENSE_STR) == 0)
00149             mm_set_dense(matcode);
00150     else
00151         return MM_UNSUPPORTED_TYPE;
00152     
00153 
00154     /* third field */
00155 
00156     if (strcmp(data_type, MM_REAL_STR) == 0)
00157         mm_set_real(matcode);
00158     else
00159     if (strcmp(data_type, MM_COMPLEX_STR) == 0)
00160         mm_set_complex(matcode);
00161     else
00162     if (strcmp(data_type, MM_PATTERN_STR) == 0)
00163         mm_set_pattern(matcode);
00164     else
00165     if (strcmp(data_type, MM_INT_STR) == 0)
00166         mm_set_integer(matcode);
00167     else
00168         return MM_UNSUPPORTED_TYPE;
00169     
00170 
00171     /* fourth field */
00172 
00173     if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
00174         mm_set_general(matcode);
00175     else
00176     if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
00177         mm_set_symmetric(matcode);
00178     else
00179     if (strcmp(storage_scheme, MM_HERM_STR) == 0)
00180         mm_set_hermitian(matcode);
00181     else
00182     if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
00183         mm_set_skew(matcode);
00184     else
00185         return MM_UNSUPPORTED_TYPE;
00186         
00187 
00188     return 0;
00189 }
00190 
00191 int mm_write_mtx_crd_size(FILE *f, int M, int N, int nz)
00192 {
00193   fprintf(f, "%d %d %d\n", M, N, nz);
00194   return 0;
00195 }
00196 
00197 int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
00198 {
00199     char line[MM_MAX_LINE_LENGTH];
00200     int num_items_read;
00201 
00202     /* set return null parameter values, in case we exit with errors */
00203     *M = *N = *nz = 0;
00204 
00205     /* now continue scanning until you reach the end-of-comments */
00206     do 
00207     {
00208         if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
00209             return MM_PREMATURE_EOF;
00210     }while (line[0] == '%');
00211 
00212     /* line[] is either blank or has M,N, nz */
00213     if (sscanf(line, "%d %d %d", M, N, nz) == 3)
00214         return 0;
00215         
00216     else
00217     do
00218     { 
00219         num_items_read = fscanf(f, "%d %d %d", M, N, nz); 
00220         if (num_items_read == EOF) return MM_PREMATURE_EOF;
00221     }
00222     while (num_items_read != 3);
00223 
00224     return 0;
00225 }
00226 
00227 
00228 int mm_read_mtx_array_size(FILE *f, int *M, int *N)
00229 {
00230     char line[MM_MAX_LINE_LENGTH];
00231     int num_items_read;
00232     /* set return null parameter values, in case we exit with errors */
00233     *M = *N = 0;
00234   
00235     /* now continue scanning until you reach the end-of-comments */
00236     do 
00237     {
00238         if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
00239             return MM_PREMATURE_EOF;
00240     }while (line[0] == '%');
00241 
00242     /* line[] is either blank or has M,N, nz */
00243     if (sscanf(line, "%d %d", M, N) == 2)
00244         return 0;
00245         
00246     else /* we have a blank line */
00247     do
00248     { 
00249         num_items_read = fscanf(f, "%d %d", M, N); 
00250         if (num_items_read == EOF) return MM_PREMATURE_EOF;
00251     }
00252     while (num_items_read != 2);
00253 
00254     return 0;
00255 }
00256 
00257 int mm_write_mtx_array_size(FILE *f, int M, int N)
00258 {
00259     fprintf(f, "%d %d\n", M, N);
00260     return 0;
00261 }
00262 
00263 
00264 
00265 /*-------------------------------------------------------------------------*/
00266 
00267 /******************************************************************/
00268 /* use when I[], J[], and val[]J, and val[] are already allocated */
00269 /******************************************************************/
00270 
00271 int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
00272         double val[], MM_typecode matcode)
00273 {
00274     int i;
00275     if (mm_is_complex(matcode))
00276     {
00277         for (i=0; i<nz; i++)
00278             if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
00279                 != 4) return MM_PREMATURE_EOF;
00280     }
00281     else if (mm_is_real(matcode))
00282     {
00283         for (i=0; i<nz; i++)
00284         {
00285             if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
00286                 != 3) return MM_PREMATURE_EOF;
00287 
00288         }
00289     }
00290 
00291     else if (mm_is_pattern(matcode))
00292     {
00293         for (i=0; i<nz; i++)
00294             if (fscanf(f, "%d %d", &I[i], &J[i])
00295                 != 2) return MM_PREMATURE_EOF;
00296     }
00297     else
00298         return MM_UNSUPPORTED_TYPE;
00299 
00300     return 0;
00301         
00302 }
00303 
00304 int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
00305         double *real, double *imag, MM_typecode matcode)
00306 {
00307     if (mm_is_complex(matcode))
00308     {
00309             if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
00310                 != 4) return MM_PREMATURE_EOF;
00311     }
00312     else if (mm_is_real(matcode))
00313     {
00314             if (fscanf(f, "%d %d %lg\n", I, J, real)
00315                 != 3) return MM_PREMATURE_EOF;
00316 
00317     }
00318 
00319     else if (mm_is_pattern(matcode))
00320     {
00321             if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
00322     }
00323     else
00324         return MM_UNSUPPORTED_TYPE;
00325 
00326     return 0;
00327         
00328 }
00329 
00330 
00331 /************************************************************************
00332     mm_read_mtx_crd()  fills M, N, nz, array of values, and return
00333                         type code, e.g. 'MCRS'
00334 
00335                         if matrix is complex, values[] is of size 2*nz,
00336                             (nz pairs of real/imaginary values)
00337 ************************************************************************/
00338 
00339 int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, 
00340         double **val, MM_typecode *matcode)
00341 {
00342     int ret_code;
00343     FILE *f;
00344 
00345     if (strcmp(fname, "stdin") == 0) f=stdin;
00346     else
00347     if ((f = fopen(fname, "r")) == NULL)
00348         return MM_COULD_NOT_READ_FILE;
00349 
00350 
00351     if ((ret_code = mm_read_banner(f, matcode)) != 0)
00352         return ret_code;
00353 
00354     if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && 
00355             mm_is_matrix(*matcode)))
00356         return MM_UNSUPPORTED_TYPE;
00357 
00358     if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
00359         return ret_code;
00360 
00361 
00362     //*I = (int *)  malloc(*nz * sizeof(int));
00363     //*J = (int *)  malloc(*nz * sizeof(int));
00364     //*val = NULL;
00365 
00366     *I = new int[*nz];
00367     *J = new int[*nz];
00368     *val = 0;
00369 
00370     if (mm_is_complex(*matcode))
00371     {
00372         //*val = (double *) malloc(*nz * 2 * sizeof(double));
00373         *val = new double[2*(*nz)];
00374         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00375                 *matcode);
00376         if (ret_code != 0) return ret_code;
00377     }
00378     else if (mm_is_real(*matcode))
00379     {
00380         //*val = (double *) malloc(*nz * sizeof(double));
00381         *val = new double[*nz];
00382         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00383                 *matcode);
00384         if (ret_code != 0) return ret_code;
00385     }
00386 
00387     else if (mm_is_pattern(*matcode))
00388     {
00389         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00390                 *matcode);
00391         if (ret_code != 0) return ret_code;
00392     }
00393 
00394     if (f != stdin) fclose(f);
00395     return 0;
00396 }
00397 
00398 int mm_write_banner(FILE *f, MM_typecode matcode)
00399 {
00400   char buffer[MM_MAX_LINE_LENGTH];
00401 
00402   mm_typecode_to_str(matcode, buffer);
00403   int ret_code;
00404   
00405   ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, buffer);
00406   //free(str);
00407   //delete [] str;
00408     return 0;
00409 }
00410 
00411 int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
00412         double val[], MM_typecode matcode)
00413 {
00414     FILE *f;
00415     int i;
00416 
00417     if (strcmp(fname, "stdout") == 0) 
00418         f = stdout;
00419     else
00420     if ((f = fopen(fname, "w")) == NULL)
00421         return MM_COULD_NOT_WRITE_FILE;
00422     
00423     /* print banner followed by typecode */
00424     char buffer[MM_MAX_LINE_LENGTH];
00425     mm_typecode_to_str(matcode, buffer);
00426     fprintf(f, "%s ", MatrixMarketBanner);
00427     fprintf(f, "%s\n", buffer);
00428 
00429     /* print matrix sizes and nonzeros */
00430     fprintf(f, "%d %d %d\n", M, N, nz);
00431 
00432     /* print values */
00433     if (mm_is_pattern(matcode))
00434         for (i=0; i<nz; i++)
00435             fprintf(f, "%d %d\n", I[i], J[i]);
00436     else
00437     if (mm_is_real(matcode))
00438         for (i=0; i<nz; i++)
00439             fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
00440     else
00441     if (mm_is_complex(matcode))
00442         for (i=0; i<nz; i++)
00443             fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i], 
00444                         val[2*i+1]);
00445     else
00446     {
00447         if (f != stdout) fclose(f);
00448         return MM_UNSUPPORTED_TYPE;
00449     }
00450 
00451     if (f !=stdout) fclose(f);
00452 
00453     return 0;
00454 }
00455     
00456 
00457  void mm_typecode_to_str(MM_typecode matcode, char * buffer)
00458 {
00459     char *types[4];
00460     int error =0;
00461 
00462     /* check for MTX type */
00463     if (mm_is_matrix(matcode)) 
00464         types[0] = MM_MTX_STR;
00465     else
00466         error=1;
00467 
00468     /* check for CRD or ARR matrix */
00469     if (mm_is_sparse(matcode))
00470         types[1] = MM_SPARSE_STR;
00471     else
00472     if (mm_is_dense(matcode))
00473         types[1] = MM_DENSE_STR;
00474     else
00475         return;
00476 
00477     /* check for element data type */
00478     if (mm_is_real(matcode))
00479         types[2] = MM_REAL_STR;
00480     else
00481     if (mm_is_complex(matcode))
00482         types[2] = MM_COMPLEX_STR;
00483     else
00484     if (mm_is_pattern(matcode))
00485         types[2] = MM_PATTERN_STR;
00486     else
00487     if (mm_is_integer(matcode))
00488         types[2] = MM_INT_STR;
00489     else
00490         return;
00491 
00492 
00493     /* check for symmetry type */
00494     if (mm_is_general(matcode))
00495         types[3] = MM_GENERAL_STR;
00496     else
00497     if (mm_is_symmetric(matcode))
00498         types[3] = MM_SYMM_STR;
00499     else 
00500     if (mm_is_hermitian(matcode))
00501         types[3] = MM_HERM_STR;
00502     else 
00503     if (mm_is_skew(matcode))
00504         types[3] = MM_SKEW_STR;
00505     else
00506         return;
00507 
00508     sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
00509     return;
00510 
00511 }
00512 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines