EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_mmio.cpp
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ***********************************************************************
00004 //
00005 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00006 //                 Copyright (2011) Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ***********************************************************************
00041 //@HEADER
00042 */
00043 
00044 /* 
00045 *   Matrix Market I/O library for ANSI C
00046 *
00047 *   See http://math.nist.gov/MatrixMarket for details.
00048 *
00049 *
00050 */
00051 
00052 /* JW - This is essentially a C file that was "converted" to C++, but still
00053    contains C function calls.  Since it is independent of the rest of
00054    the package, we will not include the central ConfigDefs file, but
00055    will rather #include the C headers that we need.  This should fix the
00056    build problem on sass2889.
00057 #include <Epetra_ConfigDefs.h> */
00058 #include <string.h>
00059 #include <stdio.h>
00060 #include <ctype.h>
00061 #include "EpetraExt_mmio.h"
00062 
00063 using namespace EpetraExt;
00064 namespace EpetraExt {
00065 int mm_read_unsymmetric_sparse(const char *fname, int *M_, int *N_, int *nz_,
00066                 double **val_, int **I_, int **J_)
00067 {
00068     FILE *f;
00069     MM_typecode matcode;
00070     int M, N, nz;
00071     int i;
00072     double *val;
00073     int *I, *J;
00074  
00075     if ((f = fopen(fname, "r")) == NULL)
00076             return -1;
00077  
00078  
00079     if (mm_read_banner(f, &matcode) != 0)
00080     {
00081         printf("mm_read_unsymetric: Could not process Matrix Market banner ");
00082         printf(" in file [%s]\n", fname);
00083         return -1;
00084     }
00085  
00086  
00087  
00088     if ( !(mm_is_real(matcode) && mm_is_matrix(matcode) &&
00089             mm_is_sparse(matcode)))
00090     {
00091       char buffer[MM_MAX_LINE_LENGTH];
00092       mm_typecode_to_str(matcode, buffer);
00093       fprintf(stderr, "Sorry, this application does not support ");
00094       fprintf(stderr, "Market Market type: [%s]\n",buffer);
00095         return -1;
00096     }
00097  
00098     /* find out size of sparse matrix: M, N, nz .... */
00099  
00100     if (mm_read_mtx_crd_size(f, &M, &N, &nz) !=0)
00101     {
00102         fprintf(stderr, "read_unsymmetric_sparse(): could not parse matrix size.\n");
00103         return -1;
00104     }
00105  
00106     *M_ = M;
00107     *N_ = N;
00108     *nz_ = nz;
00109  
00110     /* reseve memory for matrices */
00111  
00112     //I = (int *) malloc(nz * sizeof(int));
00113     //J = (int *) malloc(nz * sizeof(int));
00114     //val = (double *) malloc(nz * sizeof(double));
00115  
00116     I = new int[nz];
00117     J = new int[nz];
00118     val = new double[nz];
00119  
00120     *val_ = val;
00121     *I_ = I;
00122     *J_ = J;
00123  
00124     /* NOTE: when reading in doubles, ANSI C requires the use of the "l"  */
00125     /*   specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
00126     /*  (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15)            */
00127  
00128     for (i=0; i<nz; i++)
00129     {
00130         fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
00131         I[i]--;  /* adjust from 1-based to 0-based */
00132         J[i]--;
00133     }
00134     fclose(f);
00135  
00136     return 0;
00137 }
00138 
00139 int mm_is_valid(MM_typecode matcode)
00140 {
00141     if (!mm_is_matrix(matcode)) return 0;
00142     if (mm_is_dense(matcode) && mm_is_pattern(matcode)) return 0;
00143     if (mm_is_real(matcode) && mm_is_hermitian(matcode)) return 0;
00144     if (mm_is_pattern(matcode) && (mm_is_hermitian(matcode) || 
00145                 mm_is_skew(matcode))) return 0;
00146     return 1;
00147 }
00148 
00149 int mm_read_banner(FILE *f, MM_typecode *matcode)
00150 {
00151     char line[MM_MAX_LINE_LENGTH];
00152     char banner[MM_MAX_TOKEN_LENGTH];
00153     char mtx[MM_MAX_TOKEN_LENGTH]; 
00154     char crd[MM_MAX_TOKEN_LENGTH];
00155     char data_type[MM_MAX_TOKEN_LENGTH];
00156     char storage_scheme[MM_MAX_TOKEN_LENGTH];
00157     char *p;
00158 
00159 
00160     mm_clear_typecode(matcode);  
00161 
00162     if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL) 
00163         return MM_PREMATURE_EOF;
00164 
00165     if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd, data_type, 
00166         storage_scheme) != 5)
00167         return MM_PREMATURE_EOF;
00168 
00169     for (p=mtx; *p!='\0'; *p=tolower(*p),p++);  /* convert to lower case */
00170     for (p=crd; *p!='\0'; *p=tolower(*p),p++);  
00171     for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
00172     for (p=storage_scheme; *p!='\0'; *p=tolower(*p),p++);
00173 
00174     /* check for banner */
00175     if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
00176         return MM_NO_HEADER;
00177 
00178     /* first field should be "mtx" */
00179     if (strcmp(mtx, MM_MTX_STR) != 0)
00180         return  MM_UNSUPPORTED_TYPE;
00181     mm_set_matrix(matcode);
00182 
00183 
00184     /* second field describes whether this is a sparse matrix (in coordinate
00185             storgae) or a dense array */
00186 
00187 
00188     if (strcmp(crd, MM_SPARSE_STR) == 0)
00189         mm_set_sparse(matcode);
00190     else
00191     if (strcmp(crd, MM_DENSE_STR) == 0)
00192             mm_set_dense(matcode);
00193     else
00194         return MM_UNSUPPORTED_TYPE;
00195     
00196 
00197     /* third field */
00198 
00199     if (strcmp(data_type, MM_REAL_STR) == 0)
00200         mm_set_real(matcode);
00201     else
00202     if (strcmp(data_type, MM_COMPLEX_STR) == 0)
00203         mm_set_complex(matcode);
00204     else
00205     if (strcmp(data_type, MM_PATTERN_STR) == 0)
00206         mm_set_pattern(matcode);
00207     else
00208     if (strcmp(data_type, MM_INT_STR) == 0)
00209         mm_set_integer(matcode);
00210     else
00211         return MM_UNSUPPORTED_TYPE;
00212     
00213 
00214     /* fourth field */
00215 
00216     if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
00217         mm_set_general(matcode);
00218     else
00219     if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
00220         mm_set_symmetric(matcode);
00221     else
00222     if (strcmp(storage_scheme, MM_HERM_STR) == 0)
00223         mm_set_hermitian(matcode);
00224     else
00225     if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
00226         mm_set_skew(matcode);
00227     else
00228         return MM_UNSUPPORTED_TYPE;
00229         
00230 
00231     return 0;
00232 }
00233 
00234 int mm_write_mtx_crd_size(FILE *f, long long M, long long N, long long nz)
00235 {
00236   fprintf(f, "%lld %lld %lld\n", M, N, nz);
00237   return 0;
00238 }
00239 
00240 int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz )
00241 {
00242     char line[MM_MAX_LINE_LENGTH];
00243     int num_items_read;
00244 
00245     /* set return null parameter values, in case we exit with errors */
00246     *M = *N = *nz = 0;
00247 
00248     /* now continue scanning until you reach the end-of-comments */
00249     do 
00250     {
00251         if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
00252             return MM_PREMATURE_EOF;
00253     }while (line[0] == '%');
00254 
00255     /* line[] is either blank or has M,N, nz */
00256     if (sscanf(line, "%d %d %d", M, N, nz) == 3)
00257         return 0;
00258         
00259     else
00260     do
00261     { 
00262         num_items_read = fscanf(f, "%d %d %d", M, N, nz); 
00263         if (num_items_read == EOF) return MM_PREMATURE_EOF;
00264     }
00265     while (num_items_read != 3);
00266 
00267     return 0;
00268 }
00269 
00270 int mm_read_mtx_crd_size(FILE *f, long long *M, long long *N, long long *nz )
00271 {
00272     char line[MM_MAX_LINE_LENGTH];
00273     int num_items_read;
00274 
00275     /* set return null parameter values, in case we exit with errors */
00276     *M = *N = *nz = 0;
00277 
00278     /* now continue scanning until you reach the end-of-comments */
00279     do 
00280     {
00281         if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
00282             return MM_PREMATURE_EOF;
00283     }while (line[0] == '%');
00284 
00285     /* line[] is either blank or has M,N, nz */
00286     if (sscanf(line, "%lld %lld %lld", M, N, nz) == 3)
00287         return 0;
00288         
00289     else
00290     do
00291     { 
00292         num_items_read = fscanf(f, "%lld %lld %lld", M, N, nz); 
00293         if (num_items_read == EOF) return MM_PREMATURE_EOF;
00294     }
00295     while (num_items_read != 3);
00296 
00297     return 0;
00298 }
00299 
00300 int mm_read_mtx_array_size(FILE *f, int *M, int *N)
00301 {
00302     char line[MM_MAX_LINE_LENGTH];
00303     int num_items_read;
00304     /* set return null parameter values, in case we exit with errors */
00305     *M = *N = 0;
00306   
00307     /* now continue scanning until you reach the end-of-comments */
00308     do 
00309     {
00310         if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL) 
00311             return MM_PREMATURE_EOF;
00312     }while (line[0] == '%');
00313 
00314     /* line[] is either blank or has M,N, nz */
00315     if (sscanf(line, "%d %d", M, N) == 2)
00316         return 0;
00317         
00318     else /* we have a blank line */
00319     do
00320     { 
00321         num_items_read = fscanf(f, "%d %d", M, N); 
00322         if (num_items_read == EOF) return MM_PREMATURE_EOF;
00323     }
00324     while (num_items_read != 2);
00325 
00326     return 0;
00327 }
00328 
00329 int mm_write_mtx_array_size(FILE *f, long long M, long long N)
00330 {
00331     fprintf(f, "%lld %lld\n", M, N);
00332     return 0;
00333 }
00334 
00335 
00336 
00337 /*-------------------------------------------------------------------------*/
00338 
00339 /******************************************************************/
00340 /* use when I[], J[], and val[]J, and val[] are already allocated */
00341 /******************************************************************/
00342 
00343 int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
00344         double val[], MM_typecode matcode)
00345 {
00346     int i;
00347     if (mm_is_complex(matcode))
00348     {
00349         for (i=0; i<nz; i++)
00350             if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
00351                 != 4) return MM_PREMATURE_EOF;
00352     }
00353     else if (mm_is_real(matcode))
00354     {
00355         for (i=0; i<nz; i++)
00356         {
00357             if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
00358                 != 3) return MM_PREMATURE_EOF;
00359 
00360         }
00361     }
00362 
00363     else if (mm_is_pattern(matcode))
00364     {
00365         for (i=0; i<nz; i++)
00366             if (fscanf(f, "%d %d", &I[i], &J[i])
00367                 != 2) return MM_PREMATURE_EOF;
00368     }
00369     else
00370         return MM_UNSUPPORTED_TYPE;
00371 
00372     return 0;
00373         
00374 }
00375 
00376 int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
00377         double *real, double *imag, MM_typecode matcode)
00378 {
00379     if (mm_is_complex(matcode))
00380     {
00381             if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
00382                 != 4) return MM_PREMATURE_EOF;
00383     }
00384     else if (mm_is_real(matcode))
00385     {
00386             if (fscanf(f, "%d %d %lg\n", I, J, real)
00387                 != 3) return MM_PREMATURE_EOF;
00388 
00389     }
00390 
00391     else if (mm_is_pattern(matcode))
00392     {
00393             if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
00394     }
00395     else
00396         return MM_UNSUPPORTED_TYPE;
00397 
00398     return 0;
00399         
00400 }
00401 
00402 int mm_read_mtx_crd_entry(FILE *f, long long *I, long long *J,
00403         double *real, double *imag, MM_typecode matcode)
00404 {
00405     if (mm_is_complex(matcode))
00406     {
00407             if (fscanf(f, "%lld %lld %lg %lg", I, J, real, imag)
00408                 != 4) return MM_PREMATURE_EOF;
00409     }
00410     else if (mm_is_real(matcode))
00411     {
00412             if (fscanf(f, "%lld %lld %lg\n", I, J, real)
00413                 != 3) return MM_PREMATURE_EOF;
00414 
00415     }
00416 
00417     else if (mm_is_pattern(matcode))
00418     {
00419             if (fscanf(f, "%lld %lld", I, J) != 2) return MM_PREMATURE_EOF;
00420     }
00421     else
00422         return MM_UNSUPPORTED_TYPE;
00423 
00424     return 0;
00425         
00426 }
00427 
00428 /************************************************************************
00429     mm_read_mtx_crd()  fills M, N, nz, array of values, and return
00430                         type code, e.g. 'MCRS'
00431 
00432                         if matrix is complex, values[] is of size 2*nz,
00433                             (nz pairs of real/imaginary values)
00434 ************************************************************************/
00435 
00436 int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, 
00437         double **val, MM_typecode *matcode)
00438 {
00439     int ret_code;
00440     FILE *f;
00441 
00442     if (strcmp(fname, "stdin") == 0) f=stdin;
00443     else
00444     if ((f = fopen(fname, "r")) == NULL)
00445         return MM_COULD_NOT_READ_FILE;
00446 
00447 
00448     if ((ret_code = mm_read_banner(f, matcode)) != 0)
00449         return ret_code;
00450 
00451     if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && 
00452             mm_is_matrix(*matcode)))
00453         return MM_UNSUPPORTED_TYPE;
00454 
00455     if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
00456         return ret_code;
00457 
00458 
00459     //*I = (int *)  malloc(*nz * sizeof(int));
00460     //*J = (int *)  malloc(*nz * sizeof(int));
00461     //*val = NULL;
00462 
00463     *I = new int[*nz];
00464     *J = new int[*nz];
00465     *val = 0;
00466 
00467     if (mm_is_complex(*matcode))
00468     {
00469         //*val = (double *) malloc(*nz * 2 * sizeof(double));
00470         *val = new double[2*(*nz)];
00471         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00472                 *matcode);
00473         if (ret_code != 0) return ret_code;
00474     }
00475     else if (mm_is_real(*matcode))
00476     {
00477         //*val = (double *) malloc(*nz * sizeof(double));
00478         *val = new double[*nz];
00479         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00480                 *matcode);
00481         if (ret_code != 0) return ret_code;
00482     }
00483 
00484     else if (mm_is_pattern(*matcode))
00485     {
00486         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00487                 *matcode);
00488         if (ret_code != 0) return ret_code;
00489     }
00490 
00491     if (f != stdin) fclose(f);
00492     return 0;
00493 }
00494 
00495 int mm_write_banner(FILE *f, MM_typecode matcode)
00496 {
00497   char buffer[MM_MAX_LINE_LENGTH];
00498 
00499   mm_typecode_to_str(matcode, buffer);
00500   int ret_code;
00501   
00502   ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, buffer);
00503   //free(str);
00504   //delete [] str;
00505     return 0;
00506 }
00507 
00508 int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
00509         double val[], MM_typecode matcode)
00510 {
00511     FILE *f;
00512     int i;
00513 
00514     if (strcmp(fname, "stdout") == 0) 
00515         f = stdout;
00516     else
00517     if ((f = fopen(fname, "w")) == NULL)
00518         return MM_COULD_NOT_WRITE_FILE;
00519     
00520     /* print banner followed by typecode */
00521     char buffer[MM_MAX_LINE_LENGTH];
00522     mm_typecode_to_str(matcode, buffer);
00523     fprintf(f, "%s ", MatrixMarketBanner);
00524     fprintf(f, "%s\n", buffer);
00525 
00526     /* print matrix sizes and nonzeros */
00527     fprintf(f, "%d %d %d\n", M, N, nz);
00528 
00529     /* print values */
00530     if (mm_is_pattern(matcode))
00531         for (i=0; i<nz; i++)
00532             fprintf(f, "%d %d\n", I[i], J[i]);
00533     else
00534     if (mm_is_real(matcode))
00535         for (i=0; i<nz; i++)
00536             fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
00537     else
00538     if (mm_is_complex(matcode))
00539         for (i=0; i<nz; i++)
00540             fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i], 
00541                         val[2*i+1]);
00542     else
00543     {
00544         if (f != stdout) fclose(f);
00545         return MM_UNSUPPORTED_TYPE;
00546     }
00547 
00548     if (f !=stdout) fclose(f);
00549 
00550     return 0;
00551 }
00552     
00553 
00554  void mm_typecode_to_str(MM_typecode matcode, char * buffer)
00555 {
00556     char *types[4];
00557     int error =0;
00558 
00559     /* check for MTX type */
00560     if (mm_is_matrix(matcode)) 
00561         types[0] = MM_MTX_STR;
00562     else
00563         error=1;
00564 
00565     /* check for CRD or ARR matrix */
00566     if (mm_is_sparse(matcode))
00567         types[1] = MM_SPARSE_STR;
00568     else
00569     if (mm_is_dense(matcode))
00570         types[1] = MM_DENSE_STR;
00571     else
00572         return;
00573 
00574     /* check for element data type */
00575     if (mm_is_real(matcode))
00576         types[2] = MM_REAL_STR;
00577     else
00578     if (mm_is_complex(matcode))
00579         types[2] = MM_COMPLEX_STR;
00580     else
00581     if (mm_is_pattern(matcode))
00582         types[2] = MM_PATTERN_STR;
00583     else
00584     if (mm_is_integer(matcode))
00585         types[2] = MM_INT_STR;
00586     else
00587         return;
00588 
00589 
00590     /* check for symmetry type */
00591     if (mm_is_general(matcode))
00592         types[3] = MM_GENERAL_STR;
00593     else
00594     if (mm_is_symmetric(matcode))
00595         types[3] = MM_SYMM_STR;
00596     else 
00597     if (mm_is_hermitian(matcode))
00598         types[3] = MM_HERM_STR;
00599     else 
00600     if (mm_is_skew(matcode))
00601         types[3] = MM_SKEW_STR;
00602     else
00603         return;
00604 
00605     sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
00606     return;
00607 
00608 }
00609 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines