EpetraExt_RowMatrixOut.cpp

Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 //@HEADER
00028 #include "EpetraExt_RowMatrixOut.h"
00029 #include "EpetraExt_mmio.h"
00030 #include "Epetra_Comm.h"
00031 #include "Epetra_Map.h"
00032 #include "Epetra_Vector.h"
00033 #include "Epetra_IntVector.h"
00034 #include "Epetra_SerialDenseVector.h"
00035 #include "Epetra_IntSerialDenseVector.h"
00036 #include "Epetra_Import.h"
00037 #include "Epetra_CrsMatrix.h"
00038 
00039 using namespace EpetraExt;
00040 namespace EpetraExt {
00041 
00042 int RowMatrixToMatlabFile( const char *filename, const Epetra_RowMatrix & A) {
00043 
00044   // Simple wrapper to make it clear what can be used to write to Matlab format
00045   return(RowMatrixToMatrixMarketFile(filename, A, 0, 0, false));
00046 }
00047 
00048 int RowMatrixToMatrixMarketFile( const char *filename, const Epetra_RowMatrix & A, 
00049          const char * matrixName,
00050          const char *matrixDescription, 
00051          bool writeHeader) {
00052   int M = A.NumGlobalRows();
00053   int N = A.NumGlobalCols();
00054   int nz = A.NumGlobalNonzeros();
00055 
00056   FILE * handle = 0;
00057 
00058   if (A.RowMatrixRowMap().Comm().MyPID()==0) { // Only PE 0 does this section
00059 
00060     handle = fopen(filename,"w");
00061     if (!handle) {EPETRA_CHK_ERR(-1);}
00062     MM_typecode matcode;
00063     mm_initialize_typecode(&matcode);
00064     mm_set_matrix(&matcode);
00065     mm_set_coordinate(&matcode);
00066     mm_set_real(&matcode);
00067 
00068     if (writeHeader==true) { // Only write header if requested (true by default)
00069     
00070       if (mm_write_banner(handle, matcode)!=0) {EPETRA_CHK_ERR(-1);}
00071       
00072       if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
00073       if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
00074       
00075       if (mm_write_mtx_crd_size(handle, M, N, nz)!=0) {EPETRA_CHK_ERR(-1);}
00076     }
00077   }
00078     
00079   if (RowMatrixToHandle(handle, A)!=0) {EPETRA_CHK_ERR(-1);}// Everybody calls this routine
00080 
00081   if (A.RowMatrixRowMap().Comm().MyPID()==0) // Only PE 0 opened a file
00082     if (fclose(handle)!=0) {EPETRA_CHK_ERR(-1);}
00083   return(0);
00084 }
00085 
00086 int RowMatrixToHandle(FILE * handle, const Epetra_RowMatrix & A) {
00087 
00088   Epetra_Map map = A.RowMatrixRowMap();
00089   const Epetra_Comm & comm = map.Comm();
00090   int numProc = comm.NumProc();
00091 
00092   if (numProc==1)
00093     writeRowMatrix(handle, A);
00094   else {
00095     int numRows = map.NumMyElements();
00096     
00097     Epetra_Map allGidsMap(-1, numRows, 0,comm);
00098     
00099     Epetra_IntVector allGids(allGidsMap);
00100     for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
00101     
00102     // Now construct a RowMatrix on PE 0 by strip-mining the rows of the input matrix A.
00103     int numChunks = numProc;
00104     int stripSize = allGids.GlobalLength()/numChunks;
00105     int remainder = allGids.GlobalLength()%numChunks;
00106     int curStart = 0;
00107     int curStripSize = 0;
00108     Epetra_IntSerialDenseVector importGidList;
00109     if (comm.MyPID()==0) 
00110       importGidList.Size(stripSize+1); // Set size of vector to max needed
00111     for (int i=0; i<numChunks; i++) {
00112       if (comm.MyPID()==0) { // Only PE 0 does this part
00113   curStripSize = stripSize;
00114   if (i<remainder) curStripSize++; // handle leftovers
00115   for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
00116   curStart += curStripSize;
00117       }
00118       // The following import map will be non-trivial only on PE 0.
00119       if (comm.MyPID()>0) assert(curStripSize==0);
00120       Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
00121       Epetra_Import gidImporter(importGidMap, allGidsMap);
00122       Epetra_IntVector importGids(importGidMap);
00123       if (importGids.Import(allGids, gidImporter, Insert)!=0) {EPETRA_CHK_ERR(-1); }
00124 
00125       // importGids now has a list of GIDs for the current strip of matrix rows.
00126       // Use these values to build another importer that will get rows of the matrix.
00127 
00128       // The following import map will be non-trivial only on PE 0.
00129       Epetra_Map importMap(-1, importGids.MyLength(), importGids.Values(), map.IndexBase(), comm);
00130       Epetra_Import importer(importMap, map);
00131       Epetra_CrsMatrix importA(Copy, importMap, 0);
00132       if (importA.Import(A, importer, Insert)!=0) {EPETRA_CHK_ERR(-1); }
00133       if (importA.FillComplete(A.OperatorDomainMap(), importMap)!=0) {EPETRA_CHK_ERR(-1);}
00134 
00135       // Finally we are ready to write this strip of the matrix to ostream
00136       if (writeRowMatrix(handle, importA)!=0) {EPETRA_CHK_ERR(-1);}
00137     }
00138   }
00139   return(0);
00140 }
00141 int writeRowMatrix(FILE * handle, const Epetra_RowMatrix & A) {
00142 
00143   int numRows = A.NumGlobalRows();
00144   Epetra_Map rowMap = A.RowMatrixRowMap();
00145   Epetra_Map colMap = A.RowMatrixColMap();
00146   const Epetra_Comm & comm = rowMap.Comm();
00147   int ioffset = 1 - rowMap.IndexBase(); // Matlab indices start at 1
00148   int joffset = 1 - colMap.IndexBase(); // Matlab indices start at 1
00149   if (comm.MyPID()!=0) {
00150     if (A.NumMyRows()!=0) {EPETRA_CHK_ERR(-1);}
00151     if (A.NumMyCols()!=0) {EPETRA_CHK_ERR(-1);}
00152   }
00153   else {
00154     if (numRows!=A.NumMyRows()) {EPETRA_CHK_ERR(-1);}
00155     Epetra_SerialDenseVector values(A.MaxNumEntries());
00156     Epetra_IntSerialDenseVector indices(A.MaxNumEntries());
00157     for (int i=0; i<numRows; i++) {
00158       int I = rowMap.GID(i) + ioffset;
00159       int numEntries;
00160       if (A.ExtractMyRowCopy(i, values.Length(), numEntries, 
00161            values.Values(), indices.Values())!=0) {EPETRA_CHK_ERR(-1);}
00162       for (int j=0; j<numEntries; j++) {
00163   int J = colMap.GID(indices[j]) + joffset;
00164   double val = values[j];
00165   fprintf(handle, "%d %d %22.16e\n", I, J, val);
00166       }
00167     }
00168   }
00169   return(0);
00170 }
00171 } // namespace EpetraExt

Generated on Thu Sep 18 12:31:44 2008 for EpetraExt by doxygen 1.3.9.1