EpetraExt_OperatorOut.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_OperatorOut.h"
00029 #include "EpetraExt_mmio.h"
00030 #include "Epetra_Comm.h"
00031 #include "Epetra_Map.h"
00032 #include "Epetra_Vector.h"
00033 #include "Epetra_MultiVector.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 OperatorToMatlabFile( const char *filename, const Epetra_Operator & A) {
00044 
00045   // Simple wrapper to make it clear what can be used to write to Matlab format
00046   EPETRA_CHK_ERR(OperatorToMatrixMarketFile(filename, A, 0, 0, false));
00047   return(0);
00048 }
00049 
00050 int OperatorToMatrixMarketFile( const char *filename, const Epetra_Operator & A, 
00051          const char * matrixName,
00052          const char *matrixDescription, 
00053          bool writeHeader) {
00054 
00055   const Epetra_Map & domainMap = A.OperatorDomainMap();
00056   const Epetra_Map & rangeMap = A.OperatorRangeMap();
00057 
00058   if (!domainMap.UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00059   if (!rangeMap.UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
00060   
00061   int M = rangeMap.NumGlobalElements();
00062   int N = domainMap.NumGlobalElements();
00063 
00064   FILE * handle = 0;
00065 
00066   // To get count of nonzero terms we do multiplies ...
00067   int nz = 0;
00068   if (get_nz(A, nz)) {EPETRA_CHK_ERR(-1);}
00069 
00070   if (domainMap.Comm().MyPID()==0) { // Only PE 0 does this section
00071 
00072     handle = fopen(filename,"w");
00073     if (!handle) {EPETRA_CHK_ERR(-1);}
00074     MM_typecode matcode;
00075     mm_initialize_typecode(&matcode);
00076     mm_set_matrix(&matcode);
00077     mm_set_coordinate(&matcode);
00078     mm_set_real(&matcode);
00079     
00080 
00081     if (writeHeader==true) { // Only write header if requested (true by default)
00082     
00083       if (mm_write_banner(handle, matcode)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00084       
00085       if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
00086       if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
00087       
00088       if (mm_write_mtx_crd_size(handle, M, N, nz)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
00089     }
00090   }
00091     
00092   if (OperatorToHandle(handle, A)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}// Everybody calls this routine
00093 
00094   if (handle!=0) fclose(handle);
00095   return(0);
00096 }
00097 
00098 int OperatorToHandle(FILE * handle, const Epetra_Operator & A) {
00099 
00100   const Epetra_Map & domainMap = A.OperatorDomainMap();
00101   const Epetra_Map & rangeMap = A.OperatorRangeMap();
00102   int N = domainMap.NumGlobalElements();
00103 
00104   //cout << "rangeMap = " << rangeMap << endl;
00105   Epetra_Map rootRangeMap = Epetra_Util::Create_Root_Map(rangeMap);
00106   //cout << "rootRangeMap = " << rootRangeMap << endl;
00107   Epetra_Map rootDomainMap = Epetra_Util::Create_Root_Map(domainMap, -1); // Replicate on all processors
00108   Epetra_Import importer(rootRangeMap, rangeMap);
00109 
00110   int chunksize = 5; // Let's do multiple RHS at a time
00111   int numchunks = N/chunksize;
00112   int rem = N%chunksize;
00113 
00114   if (rem>0) {
00115     Epetra_MultiVector xrem(domainMap, rem);
00116     Epetra_MultiVector yrem(rangeMap, rem);
00117     Epetra_MultiVector yrem1(rootRangeMap, rem);
00118     // Put 1's in slots
00119     for (int j=0; j<rem; j++) {
00120       int curGlobalCol = rootDomainMap.GID(j); // Should return same value on all processors
00121       if (domainMap.MyGID(curGlobalCol)) {
00122   int curCol = domainMap.LID(curGlobalCol);
00123   xrem[j][curCol] = 1.0;
00124       }
00125     }
00126     EPETRA_CHK_ERR(A.Apply(xrem, yrem));
00127     EPETRA_CHK_ERR(yrem1.Import(yrem, importer, Insert));
00128     EPETRA_CHK_ERR(writeOperatorStrip(handle, yrem1, rootDomainMap, rootRangeMap, 0));
00129   }
00130 
00131   if (numchunks>0) {
00132     Epetra_MultiVector x(domainMap, chunksize);
00133     Epetra_MultiVector y(rangeMap, chunksize);
00134     Epetra_MultiVector y1(rootRangeMap, chunksize);
00135     for (int ichunk = 0; ichunk<numchunks; ichunk++) {
00136       int startCol = ichunk*chunksize+rem;
00137       // Put 1's in slots
00138       for (int j=0; j<chunksize; j++) {
00139   int curGlobalCol = rootDomainMap.GID(startCol+j); // Should return same value on all processors
00140   if (domainMap.MyGID(curGlobalCol)){
00141     int curCol = domainMap.LID(curGlobalCol);
00142     x[j][curCol] = 1.0;
00143   }
00144       }
00145       EPETRA_CHK_ERR(A.Apply(x, y));
00146       EPETRA_CHK_ERR(y1.Import(y, importer, Insert));
00147       EPETRA_CHK_ERR(writeOperatorStrip(handle, y1, rootDomainMap, rootRangeMap, startCol));
00148       // Put 0's in slots
00149       for (int j=0; j<chunksize; j++) {
00150   int curGlobalCol = rootDomainMap.GID(startCol+j); // Should return same value on all processors
00151   if (domainMap.MyGID(curGlobalCol)){
00152     int curCol = domainMap.LID(curGlobalCol);
00153     x[j][curCol] = 0.0;
00154   }
00155       }
00156     }
00157   }
00158 
00159   return(0);
00160 }
00161 int writeOperatorStrip(FILE * handle, const Epetra_MultiVector & y, const Epetra_Map & rootDomainMap, const Epetra_Map & rootRangeMap, int startColumn) {
00162 
00163   int numRows = y.GlobalLength();
00164   int numCols = y.NumVectors();
00165   int ioffset = 1 - rootRangeMap.IndexBase(); // Matlab indices start at 1
00166   int joffset = 1 - rootDomainMap.IndexBase(); // Matlab indices start at 1
00167   if (y.Comm().MyPID()!=0) {
00168     if (y.MyLength()!=0) {EPETRA_CHK_ERR(-1);}
00169   }
00170   else {
00171     if (numRows!=y.MyLength()) {EPETRA_CHK_ERR(-1);}
00172     for (int j=0; j<numCols; j++) {
00173       int J = rootDomainMap.GID(j + startColumn) + joffset;
00174       for (int i=0; i<numRows; i++) {
00175   double val = y[j][i];
00176   if (val!=0.0) {
00177     int I = rootRangeMap.GID(i) + ioffset;
00178     fprintf(handle, "%d %d %22.16e\n", I, J, val);
00179   }
00180       }
00181     }
00182   }
00183   return(0);
00184 }
00185 int get_nz(const Epetra_Operator & A, int & nz) {
00186   
00187   const Epetra_Map & domainMap = A.OperatorDomainMap();
00188   const Epetra_Map & rangeMap = A.OperatorRangeMap();
00189     
00190 
00191   int N = domainMap.NumGlobalElements();
00192   Epetra_Map rootDomainMap = Epetra_Util::Create_Root_Map(domainMap, -1); // Replicate on all processors
00193 
00194 
00195   int chunksize = 5; // Let's do multiple RHS at a time
00196   int numchunks = N/chunksize;
00197   int rem = N%chunksize;
00198 
00199   int lnz = 0;
00200   if (rem>0) {
00201     Epetra_MultiVector xrem(domainMap, rem);
00202     Epetra_MultiVector yrem(rangeMap, rem);
00203     // Put 1's in slots
00204     for (int j=0; j<rem; j++) {
00205       int curGlobalCol = rootDomainMap.GID(j);
00206       if (domainMap.MyGID(curGlobalCol)) xrem[j][domainMap.LID(curGlobalCol)] = 1.0;
00207     }
00208     EPETRA_CHK_ERR(A.Apply(xrem, yrem));
00209     for (int j=0; j<rem; j++) {
00210       int mylength = yrem.MyLength();
00211       for (int i=0; i<mylength; i++) 
00212   if (yrem[j][i]!=0.0) lnz++;
00213     }
00214   }
00215 
00216   if (numchunks>0) {
00217     Epetra_MultiVector x(domainMap, chunksize);
00218     Epetra_MultiVector y(rangeMap, chunksize);
00219     for (int ichunk = 0; ichunk<numchunks; ichunk++) {
00220       int startCol = ichunk*chunksize+rem;
00221       // Put 1's in slots
00222       for (int j=0; j<chunksize; j++) {
00223   int curGlobalCol = rootDomainMap.GID(startCol+j);
00224   if (domainMap.MyGID(curGlobalCol)) x[j][domainMap.LID(curGlobalCol)] = 1.0;
00225       }
00226       EPETRA_CHK_ERR(A.Apply(x, y));
00227       for (int j=0; j<chunksize; j++) {
00228   int mylength = y.MyLength();
00229   for (int i=0; i<mylength; i++) 
00230     if (y[j][i]!=0.0) lnz++;
00231       }
00232       // Put 0's in slots
00233       for (int j=0; j<chunksize; j++) {
00234   int curGlobalCol = rootDomainMap.GID(startCol+j);
00235   if (domainMap.MyGID(curGlobalCol)) x[j][domainMap.LID(curGlobalCol)] = 0.0;
00236       }
00237     }
00238   }
00239     
00240   // Sum up nonzero counts
00241   EPETRA_CHK_ERR(A.Comm().SumAll(&lnz, &nz, 1));
00242   return(0);
00243 }
00244 } // namespace EpetraExt

Generated on Tue Jul 13 09:23:06 2010 for EpetraExt by  doxygen 1.4.7