EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_MultiVectorOut.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_MultiVectorOut.h"
00029 #include "EpetraExt_mmio.h"
00030 #include "Epetra_Comm.h"
00031 #include "Epetra_BlockMap.h"
00032 #include "Epetra_Map.h"
00033 #include "Epetra_Vector.h"
00034 #include "Epetra_IntVector.h"
00035 #include "Epetra_SerialDenseVector.h"
00036 #include "Epetra_IntSerialDenseVector.h"
00037 #include "Epetra_Import.h"
00038 #include "Epetra_CrsMatrix.h"
00039 
00040 using namespace EpetraExt;
00041 namespace EpetraExt {
00042 
00043 int MultiVectorToMatlabFile( const char *filename, const Epetra_MultiVector & A) {
00044 
00045   std::FILE * handle = 0;
00046   if (A.Map().Comm().MyPID()==0) { // Only PE 0 does this section
00047     handle = fopen(filename,"w");
00048     if (!handle) return(-1);
00049   }
00050   if (MultiVectorToMatlabHandle(handle, A)) return(-1); // Everybody calls this routine
00051 
00052   if (A.Map().Comm().MyPID()==0) // Only PE 0 opened a file
00053     if (fclose(handle)) return(-1);
00054   return(0);
00055 }
00056 
00057 int MultiVectorToMatrixMarketFile( const char *filename, const Epetra_MultiVector & A, 
00058          const char * matrixName,
00059          const char *matrixDescription, 
00060          bool writeHeader) {
00061   int M = A.GlobalLength();
00062   int N = A.NumVectors();
00063 
00064   FILE * handle = 0;
00065 
00066   if (A.Map().Comm().MyPID()==0) { // Only PE 0 does this section
00067 
00068     handle = fopen(filename,"w");
00069     if (!handle) return(-1);
00070     MM_typecode matcode;
00071     mm_initialize_typecode(&matcode);
00072     mm_set_matrix(&matcode);
00073     mm_set_array(&matcode);
00074     mm_set_real(&matcode);
00075 
00076     if (writeHeader==true) { // Only write header if requested (true by default)
00077     
00078       if (mm_write_banner(handle, matcode)) return(-1);
00079       
00080       if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
00081       if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
00082       
00083       if (mm_write_mtx_array_size(handle, M, N)) return(-1);
00084     }
00085   }
00086     
00087   if (MultiVectorToMatrixMarketHandle(handle, A)) return(-1); // Everybody calls this routine
00088 
00089   if (A.Map().Comm().MyPID()==0) // Only PE 0 opened a file
00090     if (fclose(handle)) return(-1);
00091   return(0);
00092 }
00093 
00094 int MultiVectorToMatlabHandle(FILE * handle, const Epetra_MultiVector & A) {
00095   return(MultiVectorToHandle(handle, A, false));
00096 }
00097 int MultiVectorToMatrixMarketHandle(FILE * handle, const Epetra_MultiVector & A) {
00098   return(MultiVectorToHandle(handle, A, true));
00099 }
00100 int MultiVectorToHandle(FILE * handle, const Epetra_MultiVector & A, bool mmFormat) {
00101 
00102   Epetra_BlockMap bmap = A.Map();
00103   const Epetra_Comm & comm = bmap.Comm();
00104   int numProc = comm.NumProc();
00105 
00106   if (numProc==1)
00107     writeMultiVector(handle, A, mmFormat);
00108   else {
00109 
00110     // In the case of more than one column in the multivector, and writing to MatrixMarket
00111     // format, we call this function recursively, passing each vector of the multivector
00112     // individually so that we can get all of it written to file before going on to the next 
00113     // multivector
00114     if (A.NumVectors()>1 && mmFormat) {
00115       for (int i=0; i<A.NumVectors(); i++)
00116   if (MultiVectorToHandle(handle, *(A(i)), mmFormat)) return(-1);
00117       return(0);
00118     }
00119 
00120     Epetra_Map map(-1, bmap.NumMyPoints(), 0, comm);
00121     // Create a veiw of this multivector using a map (instead of block map)
00122     Epetra_MultiVector A1(View, map, A.Pointers(), A.NumVectors());
00123     int numRows = map.NumMyElements();
00124     
00125     Epetra_Map allGidsMap(-1, numRows, 0,comm);
00126     
00127     Epetra_IntVector allGids(allGidsMap);
00128     for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
00129     
00130     // Now construct a MultiVector on PE 0 by strip-mining the rows of the input matrix A.
00131     int numChunks = numProc;
00132     int stripSize = allGids.GlobalLength()/numChunks;
00133     int remainder = allGids.GlobalLength()%numChunks;
00134     int curStart = 0;
00135     int curStripSize = 0;
00136     Epetra_IntSerialDenseVector importGidList;
00137     if (comm.MyPID()==0) 
00138       importGidList.Size(stripSize+1); // Set size of vector to max needed
00139     for (int i=0; i<numChunks; i++) {
00140       if (comm.MyPID()==0) { // Only PE 0 does this part
00141   curStripSize = stripSize;
00142   if (i<remainder) curStripSize++; // handle leftovers
00143   for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
00144   curStart += curStripSize;
00145       }
00146       // The following import map will be non-trivial only on PE 0.
00147       Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
00148       Epetra_Import gidImporter(importGidMap, allGidsMap);
00149       Epetra_IntVector importGids(importGidMap);
00150       if (importGids.Import(allGids, gidImporter, Insert)) return(-1); 
00151 
00152       // importGids now has a list of GIDs for the current strip of matrix rows.
00153       // Use these values to build another importer that will get rows of the matrix.
00154 
00155       // The following import map will be non-trivial only on PE 0.
00156       Epetra_Map importMap(-1, importGids.MyLength(), importGids.Values(), 0, comm);
00157       Epetra_Import importer(importMap, map);
00158       Epetra_MultiVector importA(importMap, A1.NumVectors());
00159       if (importA.Import(A1, importer, Insert)) return(-1); 
00160 
00161       // Finally we are ready to write this strip of the matrix to ostream
00162       if (writeMultiVector(handle, importA, mmFormat)) return(-1);
00163     }
00164   }
00165   return(0);
00166 }
00167 int writeMultiVector(FILE * handle, const Epetra_MultiVector & A, bool mmFormat) {
00168 
00169   int ierr = 0;
00170   int length = A.GlobalLength();
00171   int numVectors = A.NumVectors();
00172   const Epetra_Comm & comm = A.Map().Comm();
00173   if (comm.MyPID()!=0) {
00174     if (A.MyLength()!=0) ierr = -1;
00175   }
00176   else {
00177     if (length!=A.MyLength()) ierr = -1;
00178     for (int j=0; j<numVectors; j++) {
00179       for (int i=0; i<length; i++) {
00180   double val = A[j][i];
00181   if (mmFormat)
00182     fprintf(handle, "%22.16e\n", val);
00183   else
00184     fprintf(handle, "%22.16e ", val);
00185       }
00186       if (!mmFormat) fprintf(handle, "%s", "\n");
00187     }
00188   }
00189   int ierrGlobal;
00190   comm.MinAll(&ierr, &ierrGlobal, 1); // If any processor has -1, all return -1
00191   return(ierrGlobal);
00192 }
00193 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines