Epetra Package Browser (Single Doxygen Collection) Development
test/FusedImportExport/cxx_main.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //               Epetra: 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 
00042 
00043 #include "Epetra_Map.h"
00044 #include "Epetra_Time.h"
00045 #include "Epetra_CrsMatrix.h"
00046 #include "Epetra_Import.h"
00047 #include "Epetra_Export.h"
00048 #include "Epetra_Vector.h"
00049 #include "Epetra_Flops.h"
00050 
00051 #ifdef EPETRA_MPI
00052 
00053 #include "Epetra_MpiComm.h"
00054 #include "mpi.h"
00055 #include "../epetra_test_err.h"
00056 #include "Epetra_Version.h"
00057 
00058 // prototypes
00059 
00060 int check(Epetra_CrsMatrix& A, int NumMyRows1, int NumGlobalRows1, int NumMyNonzeros1,
00061     int NumGlobalNonzeros1, int * MyGlobalElements, bool verbose);
00062 
00063 int power_method(bool TransA, Epetra_CrsMatrix& A, 
00064      Epetra_Vector& q,
00065      Epetra_Vector& z, 
00066      Epetra_Vector& resid, 
00067      double * lambda, int niters, double tolerance,
00068      bool verbose);
00069 
00070 int check_graph_sharing(Epetra_Comm& Comm);
00071 
00072 
00073 double test_with_matvec(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B){
00074   const Epetra_Map & Xamap  = A.DomainMap();
00075   const Epetra_Map & Yamap  = A.RangeMap();
00076   const Epetra_Map & Xbmap  = B.DomainMap();
00077   const Epetra_Map & Ybmap  = B.RangeMap();
00078 
00079   Epetra_Vector Xa(Xamap), Xb(Xbmap), Ya(Yamap), Yb(Ybmap), Diff(Yamap);
00080 
00081   Xa.SetSeed(24601);
00082   Xa.Random();
00083 
00084   // Handle domain map change
00085   if(!Xamap.SameAs(Xbmap)) {
00086     Epetra_Import Ximport(Xbmap,Xamap);
00087     Xb.Import(Xa,Ximport,Insert);
00088   }
00089   else {
00090     Xb=Xa;
00091   }
00092   
00093   // Do the multiplies
00094   A.Apply(Xa,Ya);
00095   B.Apply(Xb,Yb);
00096 
00097   // Handle Rangemap change
00098   if(!Yamap.SameAs(Ybmap)) {
00099     Epetra_Import Yimport(Yamap,Ybmap);
00100     Diff.Import(Yb,Yimport,Insert);
00101   }
00102   else {
00103     Diff=Yb;
00104   }
00105 
00106   // Check solution
00107   Diff.Update(-1.0,Ya,1.0);
00108   double norm;
00109   Diff.Norm2(&norm);
00110 
00111   return norm;
00112 }
00113 
00114 
00115 // B here is the "reduced" matrix.  Square matrices w/ Row=Domain=Range only.
00116 double test_with_matvec_reduced_maps(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B, const Epetra_Map & Bfullmap){
00117   const Epetra_Map & Amap  = A.DomainMap();
00118   Epetra_Vector Xa(Amap), Ya(Amap), Diff(Amap);
00119   const Epetra_Map *Bmap  = Bfullmap.NumMyElements() > 0 ? &B.DomainMap() : 0;
00120   Epetra_Vector *Xb = Bmap ? new Epetra_Vector(*Bmap) : 0;
00121   Epetra_Vector *Yb = Bmap ? new Epetra_Vector(*Bmap) : 0;
00122   
00123   Epetra_Vector Xb_alias(View,Bfullmap, Bmap ? Xb->Values(): 0);
00124   Epetra_Vector Yb_alias(View,Bfullmap, Bmap ? Yb->Values(): 0);
00125 
00126   Epetra_Import Ximport(Bfullmap,Amap);
00127 
00128   // Set the input vector
00129   Xa.SetSeed(24601);
00130   Xa.Random();
00131   Xb_alias.Import(Xa,Ximport,Insert);
00132   
00133   // Do the multiplies
00134   A.Apply(Xa,Ya);
00135   if(Bmap) B.Apply(*Xb,*Yb);
00136   
00137   // Check solution
00138   Epetra_Import Yimport(Amap,Bfullmap);
00139   Diff.Import(Yb_alias,Yimport,Insert);
00140 
00141 
00142   Diff.Update(-1.0,Ya,1.0);
00143   double norm;
00144   Diff.Norm2(&norm);
00145 
00146   delete Xb; delete Yb;
00147   return norm;
00148 }
00149 
00150 
00151 
00152 int build_matrix_unfused(const Epetra_CrsMatrix & SourceMatrix, Epetra_Import & RowImporter, Epetra_CrsMatrix *&A){
00153   int rv=0;
00154   rv=A->Import(SourceMatrix, RowImporter, Insert);
00155   if(rv) {cerr<<"build_matrix_unfused: Import failed"<<endl; return rv;}
00156 
00157   rv=A->FillComplete(SourceMatrix.DomainMap(), SourceMatrix.RangeMap()); 
00158   return rv;
00159 }
00160 
00161 int build_matrix_unfused(const Epetra_CrsMatrix & SourceMatrix, Epetra_Export & RowExporter, Epetra_CrsMatrix *&A){
00162   int rv=0;
00163   rv=A->Export(SourceMatrix, RowExporter, Insert);
00164   if(rv) {cerr<<"build_matrix_unfused: Export failed"<<endl; return rv;}
00165 
00166   rv=A->FillComplete(SourceMatrix.DomainMap(), SourceMatrix.RangeMap()); 
00167   return rv;
00168 }
00169 
00170 
00171 
00172 int build_test_matrix(Epetra_MpiComm & Comm, int test_number, Epetra_CrsMatrix *&A){
00173   int NumProc = Comm.NumProc();
00174   int MyPID   = Comm.MyPID();
00175 
00176   if(test_number==1){
00177     // Case 1: Tridiagonal
00178     int NumMyEquations = 100;
00179 
00180     int NumGlobalEquations = (NumMyEquations * NumProc) + EPETRA_MIN(NumProc,3);
00181     if(MyPID < 3)  NumMyEquations++;
00182     
00183     // Construct a Map that puts approximately the same Number of equations on each processor
00184     Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
00185     
00186     // Get update list and number of local equations from newly created Map
00187     int* MyGlobalElements = new int[Map.NumMyElements()];
00188     Map.MyGlobalElements(MyGlobalElements);
00189     
00190     // Create an integer vector NumNz that is used to build the Petra Matrix.
00191     // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation on this processor
00192     
00193     int* NumNz = new int[NumMyEquations];
00194     
00195     // We are building a tridiagonal matrix where each row has (-1 2 -1)
00196     // So we need 2 off-diagonal terms (except for the first and last equation)
00197     
00198     for (int i = 0; i < NumMyEquations; i++)
00199       if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
00200   NumNz[i] = 1;
00201       else
00202   NumNz[i] = 2;
00203     
00204     // Create a Epetra_Matrix     
00205     A=new Epetra_CrsMatrix(Copy, Map, NumNz);
00206     
00207     // Add  rows one-at-a-time
00208     // Need some vectors to help
00209     // Off diagonal Values will always be -1
00210     
00211     double* Values = new double[2];
00212     Values[0] = -1.0; 
00213     Values[1] = -1.0;
00214     int* Indices = new int[2];
00215     double two = 2.0;
00216     int NumEntries;
00217     
00218     for (int i = 0; i < NumMyEquations; i++) {
00219       if(MyGlobalElements[i] == 0) {
00220   Indices[0] = 1;
00221   NumEntries = 1;
00222       }
00223       else if (MyGlobalElements[i] == NumGlobalEquations-1) {
00224   Indices[0] = NumGlobalEquations-2;
00225   NumEntries = 1;
00226       }
00227       else {
00228   Indices[0] = MyGlobalElements[i]-1;
00229   Indices[1] = MyGlobalElements[i]+1;
00230   NumEntries = 2;
00231       }
00232       A->InsertGlobalValues(MyGlobalElements[i], NumEntries, Values, Indices);
00233       A->InsertGlobalValues(MyGlobalElements[i], 1, &two, MyGlobalElements+i);
00234     }
00235     
00236     A->FillComplete();
00237     
00238     // Cleanup
00239     delete [] MyGlobalElements;
00240     delete [] NumNz;
00241     delete [] Values;
00242     delete [] Indices;
00243     
00244   }
00245   return 0; 
00246 }
00247 
00248 
00249 void build_test_map(const Epetra_Map & oldMap, Epetra_Map *& newMap) {
00250   int NumProc = oldMap.Comm().NumProc();
00251   int MyPID   = oldMap.Comm().MyPID();
00252 
00253   int num_global = oldMap.NumGlobalElements();
00254   if(NumProc<3) {
00255     // Dump everything onto -proc 0
00256     int num_local = MyPID==0 ? num_global : 0;
00257     newMap = new Epetra_Map(num_global,num_local,0,oldMap.Comm());
00258   }
00259   else {
00260     // Split everything between procs 0 and 2 (leave proc 1 empty)
00261     int num_local=0;
00262     if(MyPID==0) num_local = num_global/2;
00263     else if(MyPID==2) num_local =  num_global - ((int)num_global/2);
00264     newMap = new Epetra_Map(num_global,num_local,0,oldMap.Comm());
00265   } 
00266 }
00267 
00268 
00269 int main(int argc, char *argv[])
00270 {
00271   int ierr = 0, total_err=0;
00272 
00273   // Initialize MPI
00274 
00275   MPI_Init(&argc,&argv);
00276   int rank; // My process ID
00277 
00278   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00279   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00280 
00281   bool verbose = false;
00282 
00283   // Check if we should print results to standard out
00284   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00285 
00286   int verbose_int = verbose ? 1 : 0;
00287   Comm.Broadcast(&verbose_int, 1, 0);
00288   verbose = verbose_int==1 ? true : false;
00289 
00290   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00291   int MyPID = Comm.MyPID();
00292   int NumProc = Comm.NumProc();
00293 
00294   if(verbose && MyPID==0)
00295     cout << Epetra_Version() << std::endl << std::endl;
00296 
00297   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
00298         << " is alive."<<endl;
00299 
00300   // Redefine verbose to only print on PE 0
00301   if(verbose && rank!=0) verbose = false;
00302 
00303   // Matrix & Map pointers
00304   Epetra_CrsMatrix *A, *B, *C;
00305   Epetra_Map* Map1;
00306   Epetra_Import* Import1;
00307   Epetra_Export* Export1;
00308   double diff_tol=1e-12;
00309 
00310 #define ENABLE_TEST_1
00311 #define ENABLE_TEST_2
00312 #define ENABLE_TEST_3
00313 #define ENABLE_TEST_4
00314 #define ENABLE_TEST_5
00315 #define ENABLE_TEST_6
00316 
00318   // Test #1: Tridiagonal Matrix; Migrate to Proc 0
00320 #ifdef ENABLE_TEST_1
00321   {
00322     double diff;
00323     ierr=build_test_matrix(Comm,1,A); 
00324     int num_global = A->RowMap().NumGlobalElements();
00325     
00326     // New map with all on Proc1
00327     if(MyPID==0) Map1=new Epetra_Map(num_global,num_global,0,Comm);
00328     else         Map1=new Epetra_Map(num_global,0,0,Comm);
00329 
00330     // Execute fused import constructor
00331     Import1 = new Epetra_Import(*Map1,A->RowMap());
00332     B=new Epetra_CrsMatrix(*A,*Import1,0,&A->RangeMap());
00333 
00334     diff=test_with_matvec(*A,*B);
00335     if(diff > diff_tol){
00336       if(MyPID==0) cout<<"FusedImport: Test #1 FAILED with norm diff = "<<diff<<"."<<endl;
00337       total_err--;
00338     }
00339     
00340     // Execute fused export constructor
00341     delete B;
00342     Export1 = new Epetra_Export(A->RowMap(),*Map1);
00343     B=new Epetra_CrsMatrix(*A,*Export1,0,&A->RangeMap());
00344 
00345     diff=test_with_matvec(*A,*B);
00346     if(diff > diff_tol){
00347       if(MyPID==0) cout<<"FusedExport: Test #1 FAILED with norm diff = "<<diff<<"."<<endl;
00348       total_err--;
00349     }
00350 
00351     delete A; delete B; delete Map1; delete Import1; delete Export1;
00352   }
00353 #endif
00354 
00355 
00357   // Test #2: Tridiagonal Matrix; Locally Reversed Map
00359 #ifdef ENABLE_TEST_2
00360   {
00361     double diff;
00362     ierr=build_test_matrix(Comm,1,A); 
00363     int num_local = A->RowMap().NumMyElements();
00364     
00365     std::vector<int> MyGIDS(num_local);
00366     for(int i=0; i<num_local; i++)
00367       MyGIDS[i] = A->RowMap().GID(num_local-i-1);
00368 
00369     // New map with all on Proc1
00370     Map1=new Epetra_Map(-1,num_local,&MyGIDS[0],0,Comm);
00371 
00372     // Execute fused import constructor
00373     Import1 = new Epetra_Import(*Map1,A->RowMap());
00374     B=new Epetra_CrsMatrix(*A,*Import1,0,&A->RangeMap());
00375 
00376     diff=test_with_matvec(*A,*B);
00377     if(diff > diff_tol){
00378       if(MyPID==0) cout<<"FusedImport: Test #2 FAILED with norm diff = "<<diff<<"."<<endl;
00379       total_err--;
00380     }
00381 
00382     // Execute fused export constructor
00383     delete B;
00384     Export1 = new Epetra_Export(A->RowMap(),*Map1);
00385     B=new Epetra_CrsMatrix(*A,*Export1,0,&A->RangeMap());
00386 
00387     diff=test_with_matvec(*A,*B);
00388     if(diff > diff_tol){
00389       if(MyPID==0) cout<<"FusedExport: Test #2 FAILED with norm diff = "<<diff<<"."<<endl;
00390       total_err--;
00391     }
00392 
00393     delete A; delete B; delete Map1; delete Import1; delete Export1;
00394   }
00395 #endif
00396 
00398   // Test #3: Tridiagonal Matrix; Globally Reversed Map
00400 #ifdef ENABLE_TEST_3
00401   {
00402     double diff;
00403     ierr=build_test_matrix(Comm,1,A); 
00404     int num_local  = A->RowMap().NumMyElements();
00405     int num_global = A->RowMap().NumGlobalElements();
00406     int num_scansum = 0;
00407 
00408     Comm.ScanSum(&num_local,&num_scansum,1);
00409 
00410     // New Map
00411     std::vector<int> MyGIDS(num_local);
00412     for(int i=0; i<num_local; i++)
00413       MyGIDS[i] = num_global - num_scansum + num_local - i - 1;
00414     Map1=new Epetra_Map(-1,num_local,&MyGIDS[0],0,Comm);
00415 
00416 
00417     // Execute fused import constructor
00418     Import1 = new Epetra_Import(*Map1,A->RowMap());   
00419     B=new Epetra_CrsMatrix(*A,*Import1,0,&A->RangeMap());
00420     
00421     diff=test_with_matvec(*A,*B);
00422     if(diff > diff_tol){
00423       if(MyPID==0) cout<<"FusedImport: Test #3 FAILED with norm diff = "<<diff<<"."<<endl;
00424       total_err--;
00425     }
00426 
00427     // Execute fused export constructor
00428     delete B;
00429     Export1 = new Epetra_Export(A->RowMap(),*Map1);
00430     B=new Epetra_CrsMatrix(*A,*Export1,0,&A->RangeMap());
00431 
00432     diff=test_with_matvec(*A,*B);
00433     if(diff > diff_tol){
00434       if(MyPID==0) cout<<"FusedExport: Test #3 FAILED with norm diff = "<<diff<<"."<<endl;
00435       total_err--;
00436     }
00437 
00438     delete A; delete B; delete Map1; delete Import1; delete Export1;
00439   }
00440 #endif
00441 
00442 
00444   // Test #4: Tridiagonal Matrix; MMM style halo import
00446 #ifdef ENABLE_TEST_4
00447   {
00448     double diff;
00449     ierr=build_test_matrix(Comm,1,A); 
00450 
00451     // Assume we always own the diagonal
00452     int num_local = A->NumMyCols()-A->NumMyRows();
00453     std::vector<int> MyGIDS(num_local);
00454 
00455     for(int i=0, idx=0; i<A->NumMyCols(); i++)
00456       if(A->LRID(A->GCID(i)) == -1){
00457   MyGIDS[idx] = A->GCID(i);
00458   idx++;
00459       }
00460 
00461     // New map
00462     const int * MyGIDS_ptr = MyGIDS.size() ? &MyGIDS[0] : 0;
00463     Map1=new Epetra_Map(-1,num_local,MyGIDS_ptr,0,Comm);
00464 
00465     
00466     // Execute fused import constructor
00467     Import1 = new Epetra_Import(*Map1,A->RowMap());
00468     B=new Epetra_CrsMatrix(*A,*Import1,0,&A->RangeMap());
00469 
00470     // Build unfused matrix to compare
00471     C=new Epetra_CrsMatrix(Copy,*Map1,0);
00472     build_matrix_unfused(*A,*Import1,C);
00473 
00474     diff=test_with_matvec(*B,*C);
00475     if(diff > diff_tol){
00476       if(MyPID==0) cout<<"FusedImport: Test #4 FAILED with norm diff = "<<diff<<"."<<endl;
00477       total_err--;
00478     }
00479 
00480     // Execute fused export constructor
00481     delete B;
00482     Export1 = new Epetra_Export(A->RowMap(),*Map1);
00483     B=new Epetra_CrsMatrix(*A,*Export1,0,&A->RangeMap());
00484 
00485     diff=test_with_matvec(*B,*C);
00486     if(diff > diff_tol){
00487       if(MyPID==0) cout<<"FusedExport: Test #4 FAILED with norm diff = "<<diff<<"."<<endl;
00488       total_err--;
00489     }
00490 
00491     delete A; delete B; delete C; delete Map1; delete Import1; delete Export1;
00492   }
00493 #endif
00494 
00495 
00497   // Test 5: Tridiagonal Matrix; Migrate to Proc 0, Replace Maps
00499 #ifdef ENABLE_TEST_5
00500   {
00501     double diff;
00502     ierr=build_test_matrix(Comm,1,A); 
00503 
00504     // New map with all on Procs 0 and 2
00505     build_test_map(A->RowMap(),Map1);    
00506 
00507     // Execute fused import constructor
00508     Import1 = new Epetra_Import(*Map1,A->RowMap());
00509     B=new Epetra_CrsMatrix(*A,*Import1,Map1,Map1);
00510 
00511     diff=test_with_matvec(*A,*B);
00512     if(diff > diff_tol){
00513       if(MyPID==0) cout<<"FusedImport: Test #5 FAILED with norm diff = "<<diff<<"."<<endl;
00514       total_err--;
00515     }
00516     
00517     // Execute fused export constructor
00518     delete B;
00519     Export1 = new Epetra_Export(A->RowMap(),*Map1);
00520     B=new Epetra_CrsMatrix(*A,*Export1,Map1,Map1);
00521     
00522     diff=test_with_matvec(*A,*B);
00523     if(diff > diff_tol){
00524       if(MyPID==0) cout<<"FusedExport: Test #5 FAILED with norm diff = "<<diff<<"."<<endl;
00525       total_err--;
00526     }
00527 
00528     delete A; delete B; delete Map1; delete Import1; delete Export1;
00529   }
00530 #endif
00531 
00532 
00534   // Test 6: Tridiagonal Matrix; Migrate to Proc 0, Replace Comm
00536 #ifdef ENABLE_TEST_6
00537   {
00538     double diff;
00539     ierr=build_test_matrix(Comm,1,A); 
00540     
00541     // New map with all on Procs 0 and 2
00542     build_test_map(A->RowMap(),Map1);
00543 
00544     // Execute fused import constructor
00545     Import1 = new Epetra_Import(*Map1,A->RowMap());
00546     B=new Epetra_CrsMatrix(*A,*Import1,Map1,Map1,true);   
00547 
00548     diff=test_with_matvec_reduced_maps(*A,*B,*Map1);
00549     if(diff > diff_tol){
00550       if(MyPID==0) cout<<"FusedImport: Test #6 FAILED with norm diff = "<<diff<<"."<<endl;
00551       total_err--;
00552     }
00553     
00554     // Execute fused export constructor
00555     delete B;
00556     Export1 = new Epetra_Export(A->RowMap(),*Map1);
00557     B=new Epetra_CrsMatrix(*A,*Export1,Map1,Map1,true);
00558     
00559     diff=test_with_matvec_reduced_maps(*A,*B,*Map1);
00560     if(diff > diff_tol){
00561       if(MyPID==0) cout<<"FusedExport: Test #6 FAILED with norm diff = "<<diff<<"."<<endl;
00562       total_err--;
00563     }
00564 
00565     delete A; delete B; delete Map1; delete Import1; delete Export1;
00566   }
00567 #endif
00568 
00569 
00570   // Final output for OK
00571   if(MyPID==0 && total_err==0)
00572     cout<<"FusedImportExport: All tests PASSED."<<endl;
00573 
00574   // Cleanup  
00575   MPI_Finalize();
00576 
00577   return total_err ;
00578 }
00579 
00580 
00581 
00582 #else
00583 int main(){
00584   
00585   return 0;
00586 }
00587 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines