EpetraExt 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 type0[20];
00557     char type1[20];
00558     char type2[20];
00559     char type3[20];
00560     int error =0;
00561 
00562     /* check for MTX type */
00563     if (mm_is_matrix(matcode))
00564         strcpy(type0, MM_MTX_STR);
00565     else
00566         error=1;
00567 
00568     /* check for CRD or ARR matrix */
00569     if (mm_is_sparse(matcode))
00570         strcpy(type1, MM_SPARSE_STR);
00571     else
00572     if (mm_is_dense(matcode))
00573         strcpy(type1, MM_DENSE_STR);
00574     else
00575         return;
00576 
00577     /* check for element data type */
00578     if (mm_is_real(matcode))
00579         strcpy(type2, MM_REAL_STR);
00580     else
00581     if (mm_is_complex(matcode))
00582         strcpy(type2, MM_COMPLEX_STR);
00583     else
00584     if (mm_is_pattern(matcode))
00585         strcpy(type2, MM_PATTERN_STR);
00586     else
00587     if (mm_is_integer(matcode))
00588         strcpy(type2, MM_INT_STR);
00589     else
00590         return;
00591 
00592     /* check for symmetry type */
00593     if (mm_is_general(matcode))
00594         strcpy(type3, MM_GENERAL_STR);
00595     else
00596     if (mm_is_symmetric(matcode))
00597         strcpy(type3, MM_SYMM_STR);
00598     else
00599     if (mm_is_hermitian(matcode))
00600         strcpy(type3, MM_HERM_STR);
00601     else
00602     if (mm_is_skew(matcode))
00603         strcpy(type3, MM_SKEW_STR);
00604     else
00605         return;
00606 
00607     sprintf(buffer,"%s %s %s %s", type0, type1, type2, type3);
00608     return;
00609 }
00610 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines