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