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, int M, int N, int nz)
00235 {
00236   fprintf(f, "%d %d %d\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 
00271 int mm_read_mtx_array_size(FILE *f, int *M, int *N)
00272 {
00273     char line[MM_MAX_LINE_LENGTH];
00274     int num_items_read;
00275     /* set return null parameter values, in case we exit with errors */
00276     *M = *N = 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, "%d %d", M, N) == 2)
00287         return 0;
00288         
00289     else /* we have a blank line */
00290     do
00291     { 
00292         num_items_read = fscanf(f, "%d %d", M, N); 
00293         if (num_items_read == EOF) return MM_PREMATURE_EOF;
00294     }
00295     while (num_items_read != 2);
00296 
00297     return 0;
00298 }
00299 
00300 int mm_write_mtx_array_size(FILE *f, int M, int N)
00301 {
00302     fprintf(f, "%d %d\n", M, N);
00303     return 0;
00304 }
00305 
00306 
00307 
00308 /*-------------------------------------------------------------------------*/
00309 
00310 /******************************************************************/
00311 /* use when I[], J[], and val[]J, and val[] are already allocated */
00312 /******************************************************************/
00313 
00314 int mm_read_mtx_crd_data(FILE *f, int M, int N, int nz, int I[], int J[],
00315         double val[], MM_typecode matcode)
00316 {
00317     int i;
00318     if (mm_is_complex(matcode))
00319     {
00320         for (i=0; i<nz; i++)
00321             if (fscanf(f, "%d %d %lg %lg", &I[i], &J[i], &val[2*i], &val[2*i+1])
00322                 != 4) return MM_PREMATURE_EOF;
00323     }
00324     else if (mm_is_real(matcode))
00325     {
00326         for (i=0; i<nz; i++)
00327         {
00328             if (fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i])
00329                 != 3) return MM_PREMATURE_EOF;
00330 
00331         }
00332     }
00333 
00334     else if (mm_is_pattern(matcode))
00335     {
00336         for (i=0; i<nz; i++)
00337             if (fscanf(f, "%d %d", &I[i], &J[i])
00338                 != 2) return MM_PREMATURE_EOF;
00339     }
00340     else
00341         return MM_UNSUPPORTED_TYPE;
00342 
00343     return 0;
00344         
00345 }
00346 
00347 int mm_read_mtx_crd_entry(FILE *f, int *I, int *J,
00348         double *real, double *imag, MM_typecode matcode)
00349 {
00350     if (mm_is_complex(matcode))
00351     {
00352             if (fscanf(f, "%d %d %lg %lg", I, J, real, imag)
00353                 != 4) return MM_PREMATURE_EOF;
00354     }
00355     else if (mm_is_real(matcode))
00356     {
00357             if (fscanf(f, "%d %d %lg\n", I, J, real)
00358                 != 3) return MM_PREMATURE_EOF;
00359 
00360     }
00361 
00362     else if (mm_is_pattern(matcode))
00363     {
00364             if (fscanf(f, "%d %d", I, J) != 2) return MM_PREMATURE_EOF;
00365     }
00366     else
00367         return MM_UNSUPPORTED_TYPE;
00368 
00369     return 0;
00370         
00371 }
00372 
00373 
00374 /************************************************************************
00375     mm_read_mtx_crd()  fills M, N, nz, array of values, and return
00376                         type code, e.g. 'MCRS'
00377 
00378                         if matrix is complex, values[] is of size 2*nz,
00379                             (nz pairs of real/imaginary values)
00380 ************************************************************************/
00381 
00382 int mm_read_mtx_crd(char *fname, int *M, int *N, int *nz, int **I, int **J, 
00383         double **val, MM_typecode *matcode)
00384 {
00385     int ret_code;
00386     FILE *f;
00387 
00388     if (strcmp(fname, "stdin") == 0) f=stdin;
00389     else
00390     if ((f = fopen(fname, "r")) == NULL)
00391         return MM_COULD_NOT_READ_FILE;
00392 
00393 
00394     if ((ret_code = mm_read_banner(f, matcode)) != 0)
00395         return ret_code;
00396 
00397     if (!(mm_is_valid(*matcode) && mm_is_sparse(*matcode) && 
00398             mm_is_matrix(*matcode)))
00399         return MM_UNSUPPORTED_TYPE;
00400 
00401     if ((ret_code = mm_read_mtx_crd_size(f, M, N, nz)) != 0)
00402         return ret_code;
00403 
00404 
00405     //*I = (int *)  malloc(*nz * sizeof(int));
00406     //*J = (int *)  malloc(*nz * sizeof(int));
00407     //*val = NULL;
00408 
00409     *I = new int[*nz];
00410     *J = new int[*nz];
00411     *val = 0;
00412 
00413     if (mm_is_complex(*matcode))
00414     {
00415         //*val = (double *) malloc(*nz * 2 * sizeof(double));
00416         *val = new double[2*(*nz)];
00417         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00418                 *matcode);
00419         if (ret_code != 0) return ret_code;
00420     }
00421     else if (mm_is_real(*matcode))
00422     {
00423         //*val = (double *) malloc(*nz * sizeof(double));
00424         *val = new double[*nz];
00425         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00426                 *matcode);
00427         if (ret_code != 0) return ret_code;
00428     }
00429 
00430     else if (mm_is_pattern(*matcode))
00431     {
00432         ret_code = mm_read_mtx_crd_data(f, *M, *N, *nz, *I, *J, *val, 
00433                 *matcode);
00434         if (ret_code != 0) return ret_code;
00435     }
00436 
00437     if (f != stdin) fclose(f);
00438     return 0;
00439 }
00440 
00441 int mm_write_banner(FILE *f, MM_typecode matcode)
00442 {
00443   char buffer[MM_MAX_LINE_LENGTH];
00444 
00445   mm_typecode_to_str(matcode, buffer);
00446   int ret_code;
00447   
00448   ret_code = fprintf(f, "%s %s\n", MatrixMarketBanner, buffer);
00449   //free(str);
00450   //delete [] str;
00451     return 0;
00452 }
00453 
00454 int mm_write_mtx_crd(char fname[], int M, int N, int nz, int I[], int J[],
00455         double val[], MM_typecode matcode)
00456 {
00457     FILE *f;
00458     int i;
00459 
00460     if (strcmp(fname, "stdout") == 0) 
00461         f = stdout;
00462     else
00463     if ((f = fopen(fname, "w")) == NULL)
00464         return MM_COULD_NOT_WRITE_FILE;
00465     
00466     /* print banner followed by typecode */
00467     char buffer[MM_MAX_LINE_LENGTH];
00468     mm_typecode_to_str(matcode, buffer);
00469     fprintf(f, "%s ", MatrixMarketBanner);
00470     fprintf(f, "%s\n", buffer);
00471 
00472     /* print matrix sizes and nonzeros */
00473     fprintf(f, "%d %d %d\n", M, N, nz);
00474 
00475     /* print values */
00476     if (mm_is_pattern(matcode))
00477         for (i=0; i<nz; i++)
00478             fprintf(f, "%d %d\n", I[i], J[i]);
00479     else
00480     if (mm_is_real(matcode))
00481         for (i=0; i<nz; i++)
00482             fprintf(f, "%d %d %20.16g\n", I[i], J[i], val[i]);
00483     else
00484     if (mm_is_complex(matcode))
00485         for (i=0; i<nz; i++)
00486             fprintf(f, "%d %d %20.16g %20.16g\n", I[i], J[i], val[2*i], 
00487                         val[2*i+1]);
00488     else
00489     {
00490         if (f != stdout) fclose(f);
00491         return MM_UNSUPPORTED_TYPE;
00492     }
00493 
00494     if (f !=stdout) fclose(f);
00495 
00496     return 0;
00497 }
00498     
00499 
00500  void mm_typecode_to_str(MM_typecode matcode, char * buffer)
00501 {
00502     char *types[4];
00503     int error =0;
00504 
00505     /* check for MTX type */
00506     if (mm_is_matrix(matcode)) 
00507         types[0] = MM_MTX_STR;
00508     else
00509         error=1;
00510 
00511     /* check for CRD or ARR matrix */
00512     if (mm_is_sparse(matcode))
00513         types[1] = MM_SPARSE_STR;
00514     else
00515     if (mm_is_dense(matcode))
00516         types[1] = MM_DENSE_STR;
00517     else
00518         return;
00519 
00520     /* check for element data type */
00521     if (mm_is_real(matcode))
00522         types[2] = MM_REAL_STR;
00523     else
00524     if (mm_is_complex(matcode))
00525         types[2] = MM_COMPLEX_STR;
00526     else
00527     if (mm_is_pattern(matcode))
00528         types[2] = MM_PATTERN_STR;
00529     else
00530     if (mm_is_integer(matcode))
00531         types[2] = MM_INT_STR;
00532     else
00533         return;
00534 
00535 
00536     /* check for symmetry type */
00537     if (mm_is_general(matcode))
00538         types[3] = MM_GENERAL_STR;
00539     else
00540     if (mm_is_symmetric(matcode))
00541         types[3] = MM_SYMM_STR;
00542     else 
00543     if (mm_is_hermitian(matcode))
00544         types[3] = MM_HERM_STR;
00545     else 
00546     if (mm_is_skew(matcode))
00547         types[3] = MM_SKEW_STR;
00548     else
00549         return;
00550 
00551     sprintf(buffer,"%s %s %s %s", types[0], types[1], types[2], types[3]);
00552     return;
00553 
00554 }
00555 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines