EpetraExt Development
EpetraExt_OperatorOut.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 #include "EpetraExt_OperatorOut.h"
00042 #include "EpetraExt_mmio.h"
00043 #include "Epetra_Comm.h"
00044 #include "Epetra_Map.h"
00045 #include "Epetra_Vector.h"
00046 #include "Epetra_MultiVector.h"
00047 #include "Epetra_IntVector.h"
00048 #include "Epetra_SerialDenseVector.h"
00049 #include "Epetra_IntSerialDenseVector.h"
00050 #include "Epetra_Import.h"
00051 #include "Epetra_CrsMatrix.h"
00052 
00053 using namespace EpetraExt;
00054 namespace EpetraExt {
00055 
00056 int OperatorToMatlabFile( const char *filename, const Epetra_Operator & A) {
00057 
00058   // Simple wrapper to make it clear what can be used to write to Matlab format
00059   EPETRA_CHK_ERR(OperatorToMatrixMarketFile(filename, A, 0, 0, false));
00060   return(0);
00061 }
00062 
00063 int OperatorToMatrixMarketFile( const char *filename, const Epetra_Operator & A, 
00064          const char * matrixName,
00065          const char *matrixDescription, 
00066          bool writeHeader) {
00067 
00068   const Epetra_Map & domainMap = A.OperatorDomainMap();
00069   const Epetra_Map & rangeMap = A.OperatorRangeMap();
00070 
00071   if (!domainMap.UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00072   if (!rangeMap.UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00073   
00074   long long M = rangeMap.NumGlobalElements64();
00075   long long N = domainMap.NumGlobalElements64();
00076   long long nz = 0;
00077 
00078   FILE * handle = 0;
00079 
00080   // To get count of nonzero terms we do multiplies ...
00081   if (get_nz(A, nz)) {EPETRA_CHK_ERR(-1);}
00082 
00083   if (domainMap.Comm().MyPID()==0) { // Only PE 0 does this section
00084 
00085     handle = fopen(filename,"w");
00086     if (!handle) {EPETRA_CHK_ERR(-1);}
00087     MM_typecode matcode;
00088     mm_initialize_typecode(&matcode);
00089     mm_set_matrix(&matcode);
00090     mm_set_coordinate(&matcode);
00091     mm_set_real(&matcode);
00092     
00093 
00094     if (writeHeader==true) { // Only write header if requested (true by default)
00095     
00096       if (mm_write_banner(handle, matcode)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00097       
00098       if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
00099       if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
00100       
00101       if (mm_write_mtx_crd_size(handle, M, N, nz)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00102     }
00103   }
00104     
00105   if (OperatorToHandle(handle, A)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}// Everybody calls this routine
00106 
00107   if (handle!=0) fclose(handle);
00108   return(0);
00109 }
00110 
00111 int OperatorToHandle(FILE * handle, const Epetra_Operator & A) {
00112 
00113   const Epetra_Map & domainMap = A.OperatorDomainMap();
00114   const Epetra_Map & rangeMap = A.OperatorRangeMap();
00115   long long N = domainMap.NumGlobalElements64();
00116 
00117   //cout << "rangeMap = " << rangeMap << endl;
00118   Epetra_Map rootRangeMap = Epetra_Util::Create_Root_Map(rangeMap);
00119   //cout << "rootRangeMap = " << rootRangeMap << endl;
00120   Epetra_Map rootDomainMap = Epetra_Util::Create_Root_Map(domainMap, -1); // Replicate on all processors
00121   Epetra_Import importer(rootRangeMap, rangeMap);
00122 
00123   int chunksize = 5; // Let's do multiple RHS at a time
00124   long long numchunks = N/chunksize;
00125   int rem = N%chunksize;
00126 
00127   if (rem>0) {
00128     Epetra_MultiVector xrem(domainMap, rem);
00129     Epetra_MultiVector yrem(rangeMap, rem);
00130     Epetra_MultiVector yrem1(rootRangeMap, rem);
00131     // Put 1's in slots
00132     for (int j=0; j<rem; j++) {
00133       long long curGlobalCol = rootDomainMap.GID64(j); // Should return same value on all processors
00134       if (domainMap.MyGID(curGlobalCol)) {
00135   int curCol = domainMap.LID(curGlobalCol);
00136   xrem[j][curCol] = 1.0;
00137       }
00138     }
00139     EPETRA_CHK_ERR(A.Apply(xrem, yrem));
00140     EPETRA_CHK_ERR(yrem1.Import(yrem, importer, Insert));
00141     EPETRA_CHK_ERR(writeOperatorStrip(handle, yrem1, rootDomainMap, rootRangeMap, 0));
00142   }
00143 
00144   if (numchunks>0) {
00145     Epetra_MultiVector x(domainMap, chunksize);
00146     Epetra_MultiVector y(rangeMap, chunksize);
00147     Epetra_MultiVector y1(rootRangeMap, chunksize);
00148     for (long long ichunk = 0; ichunk<numchunks; ichunk++) {
00149       long long startCol = ichunk*chunksize+rem;
00150       // Put 1's in slots
00151       for (int j=0; j<chunksize; j++) {
00152   long long curGlobalCol = rootDomainMap.GID64(startCol+j); // Should return same value on all processors
00153   if (domainMap.MyGID(curGlobalCol)){
00154     int curCol = domainMap.LID(curGlobalCol);
00155     x[j][curCol] = 1.0;
00156   }
00157       }
00158       EPETRA_CHK_ERR(A.Apply(x, y));
00159       EPETRA_CHK_ERR(y1.Import(y, importer, Insert));
00160       EPETRA_CHK_ERR(writeOperatorStrip(handle, y1, rootDomainMap, rootRangeMap, startCol));
00161       // Put 0's in slots
00162       for (int j=0; j<chunksize; j++) {
00163   long long curGlobalCol = rootDomainMap.GID64(startCol+j); // Should return same value on all processors
00164   if (domainMap.MyGID(curGlobalCol)){
00165     int curCol = domainMap.LID(curGlobalCol);
00166     x[j][curCol] = 0.0;
00167   }
00168       }
00169     }
00170   }
00171 
00172   return(0);
00173 }
00174 int writeOperatorStrip(FILE * handle, const Epetra_MultiVector & y, const Epetra_Map & rootDomainMap, const Epetra_Map & rootRangeMap, long long startColumn) {
00175 
00176   long long numRows = y.GlobalLength64();
00177   int numCols = y.NumVectors();
00178   long long ioffset = 1 - rootRangeMap.IndexBase64(); // Matlab indices start at 1
00179   long long joffset = 1 - rootDomainMap.IndexBase64(); // Matlab indices start at 1
00180   if (y.Comm().MyPID()!=0) {
00181     if (y.MyLength()!=0) {EPETRA_CHK_ERR(-1);}
00182   }
00183   else {
00184     if (numRows!=y.MyLength()) {EPETRA_CHK_ERR(-1);}
00185     for (int j=0; j<numCols; j++) {
00186       long long J = rootDomainMap.GID64(j + startColumn) + joffset;
00187       for (long long i=0; i<numRows; i++) {
00188   double val = y[j][i];
00189   if (val!=0.0) {
00190     long long I = rootRangeMap.GID64(i) + ioffset;
00191     fprintf(handle, "%lld %lld %22.16e\n", I, J, val);
00192   }
00193       }
00194     }
00195   }
00196   return(0);
00197 }
00198 int get_nz(const Epetra_Operator & A, long long & nz) {
00199   
00200   const Epetra_Map & domainMap = A.OperatorDomainMap();
00201   const Epetra_Map & rangeMap = A.OperatorRangeMap();
00202     
00203   long long N = domainMap.NumGlobalElements64();
00204 
00205   Epetra_Map rootDomainMap = Epetra_Util::Create_Root_Map(domainMap, -1); // Replicate on all processors
00206 
00207 
00208   int chunksize = 5; // Let's do multiple RHS at a time
00209   long long numchunks = N/chunksize;
00210   int rem = N%chunksize;
00211 
00212   long long lnz = 0;
00213   if (rem>0) {
00214     Epetra_MultiVector xrem(domainMap, rem);
00215     Epetra_MultiVector yrem(rangeMap, rem);
00216     // Put 1's in slots
00217     for (int j=0; j<rem; j++) {
00218       long long curGlobalCol = rootDomainMap.GID64(j);
00219       if (domainMap.MyGID(curGlobalCol)) xrem[j][domainMap.LID(curGlobalCol)] = 1.0;
00220     }
00221     EPETRA_CHK_ERR(A.Apply(xrem, yrem));
00222     for (int j=0; j<rem; j++) {
00223       int mylength = yrem.MyLength();
00224       for (int i=0; i<mylength; i++) 
00225   if (yrem[j][i]!=0.0) lnz++;
00226     }
00227   }
00228 
00229   if (numchunks>0) {
00230     Epetra_MultiVector x(domainMap, chunksize);
00231     Epetra_MultiVector y(rangeMap, chunksize);
00232     for (long long ichunk = 0; ichunk<numchunks; ichunk++) {
00233       long long startCol = ichunk*chunksize+rem;
00234       // Put 1's in slots
00235       for (int j=0; j<chunksize; j++) {
00236   long long curGlobalCol = rootDomainMap.GID64(startCol+j);
00237   if (domainMap.MyGID(curGlobalCol)) x[j][domainMap.LID(curGlobalCol)] = 1.0;
00238       }
00239       EPETRA_CHK_ERR(A.Apply(x, y));
00240       for (int j=0; j<chunksize; j++) {
00241   int mylength = y.MyLength();
00242   for (int i=0; i<mylength; i++) 
00243     if (y[j][i]!=0.0) lnz++;
00244       }
00245       // Put 0's in slots
00246       for (int j=0; j<chunksize; j++) {
00247   long long curGlobalCol = rootDomainMap.GID64(startCol+j);
00248   if (domainMap.MyGID(curGlobalCol)) x[j][domainMap.LID(curGlobalCol)] = 0.0;
00249       }
00250     }
00251   }
00252     
00253   // Sum up nonzero counts
00254   EPETRA_CHK_ERR(A.Comm().SumAll(&lnz, &nz, 1));
00255   return(0);
00256 }
00257 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines