Epetra Package Browser (Single Doxygen Collection) Development
test/ImportExport/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_Vector.h"
00047 #include "Epetra_IntVector.h"
00048 #include "Epetra_Import.h"
00049 #include "Epetra_Export.h"
00050 #include "Epetra_OffsetIndex.h"
00051 #ifdef EPETRA_MPI
00052 #include "Epetra_MpiComm.h"
00053 #include "mpi.h"
00054 #else
00055 #include "Epetra_SerialComm.h"
00056 #endif
00057 #ifndef __cplusplus
00058 #define __cplusplus
00059 #endif
00060 #include "../epetra_test_err.h"
00061 #include "Epetra_Version.h"
00062 
00063 int special_submap_import_test(Epetra_Comm& Comm);
00064 int combine_mode_test(Epetra_Comm& Comm);
00065 int alternate_import_constructor_test(Epetra_Comm& Comm);
00066 
00067 int main(int argc, char *argv[])
00068 {
00069   int ierr = 0, i, j, forierr = 0;
00070 
00071 #ifdef EPETRA_MPI
00072   // Initialize MPI
00073   MPI_Init(&argc,&argv);
00074   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00075 #else
00076   Epetra_SerialComm Comm;
00077 #endif
00078 
00079   bool verbose = false;
00080 
00081   // Check if we should print results to standard out
00082   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00083 
00084 
00085 
00086 
00087   //char tmp;
00088   //if (Comm.MyPID()==0) cout << "Press any key to continue..."<< endl;
00089   //if (Comm.MyPID()==0) cin >> tmp;
00090   //Comm.Barrier();
00091 
00092   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00093   int MyPID = Comm.MyPID();
00094   int NumProc = Comm.NumProc();
00095 
00096   if (verbose && MyPID==0)
00097     cout << Epetra_Version() << endl << endl;
00098 
00099   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
00100               << " is alive."<<endl;
00101 
00102   // Redefine verbose to only print on PE 0
00103   if (verbose && Comm.MyPID()!=0) verbose = false;
00104 
00105   int NumMyEquations = 20;
00106   int NumGlobalEquations = NumMyEquations*NumProc+EPETRA_MIN(NumProc,3);
00107   if (MyPID < 3) NumMyEquations++;
00108   // Construct a Source Map that puts approximately the same Number of equations on each processor in 
00109   // uniform global ordering
00110 
00111   Epetra_Map SourceMap(NumGlobalEquations, NumMyEquations, 0, Comm);
00112   
00113   // Get update list and number of local equations from newly created Map
00114   int NumMyElements = SourceMap.NumMyElements();
00115   int * SourceMyGlobalElements = new int[NumMyElements];
00116   SourceMap.MyGlobalElements(SourceMyGlobalElements);
00117 
00118   // Construct a Target Map that will contain:
00119   //  some unchanged elements (relative to the soure map),
00120   //  some permuted elements
00121   //  some off-processor elements
00122   Epetra_Vector RandVec(SourceMap);
00123   RandVec.Random(); // This creates a vector of random numbers between negative one and one.
00124 
00125   int *TargetMyGlobalElements = new int[NumMyElements];
00126 
00127   int MinGID = SourceMap.MinMyGID();
00128   for (i=0; i< NumMyEquations/2; i++) TargetMyGlobalElements[i] = i + MinGID; // Half will be the same...
00129   for (i=NumMyEquations/2; i<NumMyEquations; i++) {
00130     int index = abs((int)(((double) (NumGlobalEquations-1) ) * RandVec[i]));
00131     TargetMyGlobalElements[i] = EPETRA_MIN(NumGlobalEquations-1,EPETRA_MAX(0,index));
00132   }
00133 
00134   int NumSameIDs = 0;
00135   int NumPermutedIDs = 0;
00136   int NumRemoteIDs = 0;
00137   bool StillContiguous = true;
00138   for (i=0; i < NumMyEquations; i++) {
00139     if (SourceMyGlobalElements[i]==TargetMyGlobalElements[i] && StillContiguous)
00140       NumSameIDs++;
00141     else if (SourceMap.MyGID(TargetMyGlobalElements[i])) {
00142       StillContiguous = false;
00143       NumPermutedIDs++;
00144     }
00145     else {
00146       StillContiguous = false;
00147       NumRemoteIDs++;
00148     }
00149   }
00150   EPETRA_TEST_ERR(!(NumMyEquations==NumSameIDs+NumPermutedIDs+NumRemoteIDs),ierr);
00151 
00152   Epetra_Map TargetMap(-1, NumMyElements, TargetMyGlobalElements, 0, Comm);
00153 
00154   // Create a multivector whose elements are GlobalID * (column number +1)
00155 
00156   int NumVectors = 3;
00157   Epetra_MultiVector SourceMultiVector(SourceMap, NumVectors);
00158   for (j=0; j < NumVectors; j++)
00159     for (i=0; i < NumMyElements; i++)
00160       SourceMultiVector[j][i] = (double) SourceMyGlobalElements[i]*(j+1);
00161 
00162   // Create a target multivector that we will fill using an Import
00163 
00164   Epetra_MultiVector TargetMultiVector(TargetMap, NumVectors);
00165 
00166   Epetra_Import Importer(TargetMap, SourceMap);
00167 
00168   EPETRA_TEST_ERR(!(TargetMultiVector.Import(SourceMultiVector, Importer, Insert)==0),ierr);
00169 
00170   // Test Target against expected values
00171   forierr = 0;
00172   for (j=0; j < NumVectors; j++)
00173     for (i=0; i < NumMyElements; i++) {
00174       if (TargetMultiVector[j][i]!= (double) TargetMyGlobalElements[i]*(j+1))
00175   cout << "TargetMultiVector["<<i<<"]["<<j<<"] = " << TargetMultiVector[j][i] 
00176        <<  "  TargetMyGlobalElements[i]*(j+1) = " <<  TargetMyGlobalElements[i]*(j+1) << endl;
00177       forierr += !(TargetMultiVector[j][i]== (double) TargetMyGlobalElements[i]*(j+1));
00178     }
00179   EPETRA_TEST_ERR(forierr,ierr);
00180 
00181   if (verbose) cout << "MultiVector Import using Importer Check OK" << endl << endl;
00182 
00183 
00185 
00186   // Now use Importer to do an export
00187 
00188   Epetra_Vector TargetVector(SourceMap);
00189   Epetra_Vector ExpectedTarget(SourceMap);
00190   Epetra_Vector SourceVector(TargetMap);
00191 
00192   NumSameIDs = Importer.NumSameIDs();
00193   int NumPermuteIDs = Importer.NumPermuteIDs();
00194   int NumExportIDs = Importer.NumExportIDs();
00195   int *PermuteFromLIDs = Importer.PermuteFromLIDs();
00196   int *ExportLIDs = Importer.ExportLIDs();
00197   int *ExportPIDs = Importer.ExportPIDs();
00198 
00199   for (i=0; i < NumSameIDs; i++) ExpectedTarget[i] = (double) (MyPID+1);
00200   for (i=0; i < NumPermuteIDs; i++) ExpectedTarget[PermuteFromLIDs[i]] = 
00201               (double) (MyPID+1);
00202   for (i=0; i < NumExportIDs; i++) ExpectedTarget[ExportLIDs[i]] += 
00203              (double) (ExportPIDs[i]+1);
00204 
00205   for (i=0; i < NumMyElements; i++) SourceVector[i] =  (double) (MyPID+1);
00206 
00207   EPETRA_TEST_ERR(!(TargetVector.Export(SourceVector, Importer, Add)==0),ierr);
00208 
00209   forierr = 0;
00210   for (i=0; i < NumMyElements; i++) {
00211     if (TargetVector[i]!= ExpectedTarget[i])
00212       cout <<  "     TargetVector["<<i<<"] = " << TargetVector[i] 
00213      <<  "   ExpectedTarget["<<i<<"] = " <<  ExpectedTarget[i] << " on PE " << MyPID << endl;
00214     forierr += !(TargetVector[i]== ExpectedTarget[i]);
00215   }
00216   EPETRA_TEST_ERR(forierr,ierr);
00217 
00218   if (verbose) cout << "Vector Export using Importer Check OK" << endl << endl;
00219 
00221   // Now use Importer to create a reverse exporter
00222   TargetVector.PutScalar(0.0);
00223   Epetra_Export ReversedImport(Importer);
00224 
00225   EPETRA_TEST_ERR(!(TargetVector.Export(SourceVector, ReversedImport, Add)==0),ierr);
00226 
00227   forierr = 0;
00228   for (i=0; i < NumMyElements; i++) {
00229     if (TargetVector[i]!= ExpectedTarget[i])
00230       cout <<  "     TargetVector["<<i<<"] = " << TargetVector[i] 
00231      <<  "   ExpectedTarget["<<i<<"] = " <<  ExpectedTarget[i] << " on PE " << MyPID << endl;
00232     forierr += !(TargetVector[i]== ExpectedTarget[i]);
00233   }
00234   EPETRA_TEST_ERR(forierr,ierr);
00235 
00236   if (verbose) cout << "Vector Export using Reversed Importer Check OK" << endl << endl;
00237 
00239   // Now use Exporter to create a reverse importer
00240   TargetVector.PutScalar(0.0);
00241   Epetra_Import ReversedExport(ReversedImport);
00242 
00243   EPETRA_TEST_ERR(!(TargetVector.Export(SourceVector, ReversedExport, Add)==0),ierr);
00244 
00245   forierr = 0;
00246   for (i=0; i < NumMyElements; i++) {
00247     if (TargetVector[i]!= ExpectedTarget[i])
00248       cout <<  "     TargetVector["<<i<<"] = " << TargetVector[i] 
00249      <<  "   ExpectedTarget["<<i<<"] = " <<  ExpectedTarget[i] << " on PE " << MyPID << endl;
00250     forierr += !(TargetVector[i]== ExpectedTarget[i]);
00251   }
00252   EPETRA_TEST_ERR(forierr,ierr);
00253 
00254   if (verbose) cout << "Vector Export using Reversed Exporter Check OK" << endl << endl;
00255 
00257   //  Build a tridiagonal system two ways:
00258   //  1) From "standard" matrix view where equations are uniquely owned.
00259   //  2) From 1D PDE view where nodes (equations) between processors are shared and partial contributions are done
00260   //     in parallel, then merged together at the end of the construction process.
00261   //
00263 
00264 
00265 
00266   // Construct a Standard Map that puts approximately the same number of equations on each processor in 
00267   // uniform global ordering
00268 
00269   Epetra_Map StandardMap(NumGlobalEquations, NumMyEquations, 0, Comm);
00270   
00271   // Get update list and number of local equations from newly created Map
00272   NumMyElements = StandardMap.NumMyElements();
00273   int * StandardMyGlobalElements = new int[NumMyElements];
00274   StandardMap.MyGlobalElements(StandardMyGlobalElements);
00275 
00276 
00277   // Create a standard Epetra_CrsGraph
00278 
00279   Epetra_CrsGraph StandardGraph(Copy, StandardMap, 3);
00280   EPETRA_TEST_ERR(StandardGraph.IndicesAreGlobal(),ierr);
00281   EPETRA_TEST_ERR(StandardGraph.IndicesAreLocal(),ierr);
00282   
00283   // Add  rows one-at-a-time
00284   // Need some vectors to help
00285   // Off diagonal Values will always be -1
00286 
00287 
00288   int *Indices = new int[2];
00289   int NumEntries;
00290   
00291   forierr = 0;
00292   for (i=0; i<NumMyEquations; i++)
00293     {
00294     if (StandardMyGlobalElements[i]==0)
00295       {
00296   Indices[0] = 1;
00297   NumEntries = 1;
00298       }
00299     else if (StandardMyGlobalElements[i] == NumGlobalEquations-1)
00300       {
00301   Indices[0] = NumGlobalEquations-2;
00302   NumEntries = 1;
00303       }
00304     else
00305       {
00306   Indices[0] = StandardMyGlobalElements[i]-1;
00307   Indices[1] = StandardMyGlobalElements[i]+1;
00308   NumEntries = 2;
00309       }
00310     forierr += !(StandardGraph.InsertGlobalIndices(StandardMyGlobalElements[i], NumEntries, Indices)==0);
00311     forierr += !(StandardGraph.InsertGlobalIndices(StandardMyGlobalElements[i], 1, StandardMyGlobalElements+i)==0); // Put in the diagonal entry
00312     }
00313   EPETRA_TEST_ERR(forierr,ierr);
00314 
00315   // Finish up
00316   EPETRA_TEST_ERR(!(StandardGraph.IndicesAreGlobal()),ierr);
00317   EPETRA_TEST_ERR(!(StandardGraph.FillComplete()==0),ierr);
00318   EPETRA_TEST_ERR(!(StandardGraph.IndicesAreLocal()),ierr);
00319   EPETRA_TEST_ERR(StandardGraph.StorageOptimized(),ierr);
00320   StandardGraph.OptimizeStorage();
00321   EPETRA_TEST_ERR(!(StandardGraph.StorageOptimized()),ierr);
00322   EPETRA_TEST_ERR(StandardGraph.UpperTriangular(),ierr);
00323   EPETRA_TEST_ERR(StandardGraph.LowerTriangular(),ierr);
00324 
00325   // Create Epetra_CrsMatrix using the just-built graph
00326 
00327   Epetra_CrsMatrix StandardMatrix(Copy, StandardGraph);
00328   EPETRA_TEST_ERR(StandardMatrix.IndicesAreGlobal(),ierr);
00329   EPETRA_TEST_ERR(!(StandardMatrix.IndicesAreLocal()),ierr);
00330   
00331   // Add  rows one-at-a-time
00332   // Need some vectors to help
00333   // Off diagonal Values will always be -1
00334 
00335 
00336   double *Values = new double[2];
00337   Values[0] = -1.0; Values[1] = -1.0;
00338   double two = 2.0;
00339   
00340   forierr = 0;
00341   for (i=0; i<NumMyEquations; i++)
00342     {
00343     if (StandardMyGlobalElements[i]==0)
00344       {
00345   Indices[0] = 1;
00346   NumEntries = 1;
00347       }
00348     else if (StandardMyGlobalElements[i] == NumGlobalEquations-1)
00349       {
00350   Indices[0] = NumGlobalEquations-2;
00351   NumEntries = 1;
00352       }
00353     else
00354       {
00355   Indices[0] = StandardMyGlobalElements[i]-1;
00356   Indices[1] = StandardMyGlobalElements[i]+1;
00357   NumEntries = 2;
00358       }
00359     forierr += !(StandardMatrix.ReplaceGlobalValues(StandardMyGlobalElements[i], NumEntries, Values, Indices)==0);
00360     // Put in the diagonal entry
00361     forierr += !(StandardMatrix.ReplaceGlobalValues(StandardMyGlobalElements[i], 1, &two, StandardMyGlobalElements+i)==0); 
00362     }
00363   EPETRA_TEST_ERR(forierr,ierr);
00364 
00365   // Finish up
00366   EPETRA_TEST_ERR(!(StandardMatrix.IndicesAreLocal()),ierr);
00367   EPETRA_TEST_ERR(!(StandardMatrix.FillComplete()==0),ierr);
00368   EPETRA_TEST_ERR(!(StandardMatrix.IndicesAreLocal()),ierr);
00369   //  EPETRA_TEST_ERR((StandardMatrix.StorageOptimized()),ierr);
00370   EPETRA_TEST_ERR((StandardMatrix.OptimizeStorage()),ierr);
00371   EPETRA_TEST_ERR(!(StandardMatrix.StorageOptimized()),ierr);
00372   EPETRA_TEST_ERR(StandardMatrix.UpperTriangular(),ierr);
00373   EPETRA_TEST_ERR(StandardMatrix.LowerTriangular(),ierr);
00374 
00375   // Construct an Overlapped Map of StandardMap that include the endpoints from two neighboring processors.
00376 
00377   int OverlapNumMyElements;
00378   int OverlapMinMyGID;
00379 
00380   OverlapNumMyElements = NumMyElements + 1;
00381   if (MyPID==0) OverlapNumMyElements--;
00382 
00383   if (MyPID==0) OverlapMinMyGID = StandardMap.MinMyGID();
00384   else OverlapMinMyGID = StandardMap.MinMyGID()-1;
00385 
00386   int * OverlapMyGlobalElements = new int[OverlapNumMyElements];
00387 
00388   for (i=0; i< OverlapNumMyElements; i++) OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
00389 
00390   Epetra_Map OverlapMap(-1, OverlapNumMyElements, OverlapMyGlobalElements, 0, Comm);
00391 
00392   // Create the Overlap Epetra_Matrix
00393 
00394   Epetra_CrsMatrix OverlapMatrix(Copy, OverlapMap, 4);
00395   EPETRA_TEST_ERR(OverlapMatrix.IndicesAreGlobal(),ierr);
00396   EPETRA_TEST_ERR(OverlapMatrix.IndicesAreLocal(),ierr);
00397   
00398   // Add  matrix element one cell at a time.
00399   // Each cell does an incoming and outgoing flux calculation
00400 
00401 
00402   double pos_one = 1.0;
00403   double neg_one = -1.0;
00404 
00405   forierr = 0;
00406   for (i=0; i<OverlapNumMyElements; i++)
00407     {
00408       int node_left = OverlapMyGlobalElements[i]-1;
00409       int node_center = node_left + 1;
00410       int node_right = node_left + 2;
00411       if (i>0) {
00412   if (node_left>-1)
00413     forierr += !(OverlapMatrix.InsertGlobalValues(node_center, 1, &neg_one, &node_left)==0);
00414   forierr += !(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
00415       }
00416       if (i<OverlapNumMyElements-1) {
00417   forierr += !(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
00418   if (node_right<NumGlobalEquations) 
00419     forierr += !(OverlapMatrix.InsertGlobalValues(node_center, 1, &neg_one, &node_right)==0);
00420       }
00421     }
00422   EPETRA_TEST_ERR(forierr,ierr);
00423 
00424   // Handle endpoints
00425   if (MyPID==0) {
00426     int node_center = 0;
00427     EPETRA_TEST_ERR(!(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0),ierr);
00428   }
00429   if (MyPID==NumProc-1) {
00430     int node_center = OverlapMyGlobalElements[OverlapNumMyElements-1];
00431     EPETRA_TEST_ERR(!(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0),ierr);
00432   }
00433     
00434   EPETRA_TEST_ERR(!(OverlapMatrix.FillComplete()==0),ierr);
00435 
00436   // Make a gathered matrix from OverlapMatrix.  It should be identical to StandardMatrix
00437 
00438   Epetra_CrsMatrix GatheredMatrix(Copy, StandardGraph);
00439   Epetra_Export Exporter(OverlapMap, StandardMap);
00440   EPETRA_TEST_ERR(!(GatheredMatrix.Export(OverlapMatrix, Exporter, Add)==0),ierr);
00441   EPETRA_TEST_ERR(!(GatheredMatrix.FillComplete()==0),ierr);
00442 
00443   // Check if entries of StandardMatrix and GatheredMatrix are identical
00444 
00445   int StandardNumEntries, GatheredNumEntries;
00446   int * StandardIndices, * GatheredIndices;
00447   double * StandardValues, * GatheredValues;
00448 
00449   int StandardNumMyNonzeros = StandardMatrix.NumMyNonzeros();
00450   int GatheredNumMyNonzeros = GatheredMatrix.NumMyNonzeros();
00451   EPETRA_TEST_ERR(!(StandardNumMyNonzeros==GatheredNumMyNonzeros),ierr);
00452 
00453   int StandardNumMyRows = StandardMatrix.NumMyRows();
00454   int GatheredNumMyRows = GatheredMatrix.NumMyRows();
00455   EPETRA_TEST_ERR(!(StandardNumMyRows==GatheredNumMyRows),ierr);
00456 
00457   forierr = 0;
00458   for (i=0; i< StandardNumMyRows; i++)
00459     {
00460       forierr += !(StandardMatrix.ExtractMyRowView(i, StandardNumEntries, StandardValues, StandardIndices)==0);
00461       forierr += !(GatheredMatrix.ExtractMyRowView(i, GatheredNumEntries, GatheredValues, GatheredIndices)==0);
00462       forierr += !(StandardNumEntries==GatheredNumEntries);
00463       for (j=0; j < StandardNumEntries; j++) {
00464   //if (StandardIndices[j]!=GatheredIndices[j])
00465   // cout << "MyPID = " << MyPID << " i = " << i << "   StandardIndices[" << j << "] = " << StandardIndices[j] 
00466   //      << "   GatheredIndices[" << j << "] = " << GatheredIndices[j] << endl;
00467   //if (StandardValues[j]!=GatheredValues[j])
00468   //cout << "MyPID = " << MyPID << " i = " << i << "    StandardValues[" << j << "] = " <<  StandardValues[j] 
00469   //     << "    GatheredValues[" << j << "] = " <<  GatheredValues[j] << endl;
00470   forierr += !(StandardIndices[j]==GatheredIndices[j]);
00471   forierr += !(StandardValues[j]==GatheredValues[j]);
00472       }
00473     }
00474   EPETRA_TEST_ERR(forierr,ierr);
00475 
00476   if (verbose) cout << "Matrix Export Check OK" << endl << endl;
00477 
00478   //Do Again with use of Epetra_OffsetIndex object for speed
00479   Epetra_OffsetIndex OffsetIndex( OverlapMatrix.Graph(), GatheredMatrix.Graph(), Exporter );
00480   EPETRA_TEST_ERR(!(GatheredMatrix.Export(OverlapMatrix, Exporter, Add)==0),ierr);
00481 
00482   if (verbose) cout << "Optimized Matrix Export Check OK" << endl << endl;
00483 
00484   bool passed;
00485   Epetra_IntVector v1(StandardMap); v1.PutValue(2);
00486   Epetra_IntVector v2(StandardMap); v2.PutValue(3);
00487 
00488   Epetra_Export identExporter(StandardMap,StandardMap); // Identity exporter
00489   EPETRA_TEST_ERR(!(v2.Export(v1, identExporter, Insert)==0),ierr);
00490   passed = (v2.MinValue()==2);
00491   EPETRA_TEST_ERR(!passed,ierr);
00492 
00493   v1.PutValue(1);
00494   Epetra_Import identImporter(StandardMap,StandardMap); // Identity importer
00495   EPETRA_TEST_ERR(!(v2.Import(v1, identExporter, Insert)==0),ierr);
00496   passed = passed && (v2.MaxValue()==1);
00497   EPETRA_TEST_ERR(!passed,ierr);
00498 
00499   if (verbose) {
00500     if (passed) cout << "Identity Import/Export Check OK" << endl << endl;
00501     else cout << "Identity Import/Export Check Failed" << endl << endl;
00502   }
00503 
00504   int NumSubMapElements = StandardMap.NumMyElements()/2;
00505   int SubStart = Comm.MyPID();
00506   NumSubMapElements = EPETRA_MIN(NumSubMapElements,StandardMap.NumMyElements()-SubStart);
00507   Epetra_Map SubMap(-1, NumSubMapElements, StandardMyGlobalElements+SubStart, 0, Comm);
00508 
00509   Epetra_IntVector v3(View, SubMap, SubMap.MyGlobalElements()); // Fill v3 with GID values for variety
00510   Epetra_Export subExporter(SubMap, StandardMap); // Export to a subset of indices of standard map
00511   EPETRA_TEST_ERR(!(v2.Export(v3,subExporter,Insert)==0),ierr);
00512 
00513   forierr = 0;
00514   for (i=0; i<SubMap.NumMyElements(); i++) {
00515     int i1 = StandardMap.LID(SubMap.GID(i));
00516     forierr += !(v3[i]==v2[i1]);
00517   }
00518   EPETRA_TEST_ERR(forierr,ierr);
00519 
00520   Epetra_Import subImporter(StandardMap, SubMap); // Import to a subset of indices of standard map
00521   EPETRA_TEST_ERR(!(v1.Import(v3,subImporter,Insert)==0),ierr);
00522 
00523   for (i=0; i<SubMap.NumMyElements(); i++) {
00524     int i1 = StandardMap.LID(SubMap.GID(i));
00525     forierr += !(v3[i]==v1[i1]);
00526   }
00527   EPETRA_TEST_ERR(forierr,ierr);
00528 
00529   if (verbose) {
00530     if (forierr==0) cout << "SubMap Import/Export Check OK" << endl << endl;
00531     else cout << "SubMap Import/Export Check Failed" << endl << endl;
00532   }
00533 
00534 
00535 #ifdef DOESNT_WORK_IN_PARALLEL
00536   forierr = special_submap_import_test(Comm);
00537   EPETRA_TEST_ERR(forierr, ierr);
00538 
00539   if (verbose) {
00540     if (forierr==0) cout << "Special SubMap Import Check OK" << endl << endl;
00541     else cout << "Special SubMap Import Check Failed" << endl << endl;
00542   }
00543 #endif
00544 
00545 
00546   forierr =  alternate_import_constructor_test(Comm);
00547   EPETRA_TEST_ERR(forierr, ierr);
00548 
00549   if (verbose) {
00550     if (forierr==0) cout << "Alternative Import Constructor Check OK" << endl << endl;
00551     else cout << "Alternative Import Constructor Check Failed" << endl << endl;
00552   }
00553 
00554 
00555   // Now let's test Importer replacement:  Coalesce to 1 proc...
00556   Epetra_CrsMatrix FunMatrix(StandardMatrix);
00557   if(Comm.NumProc()!=1) {
00558     forierr=0;
00559     long long num_global_elements1 = FunMatrix.DomainMap().NumGlobalElements64();
00560     long long num_global_elements2 = FunMatrix.DomainMap().MaxAllGID64()- FunMatrix.DomainMap().MinAllGID64()+1;
00561     if(num_global_elements1 == num_global_elements2) {
00562       // The original domain map is linear.  Let's have fun
00563       int NumMyElements = Comm.MyPID()==0 ? num_global_elements1 : 0;
00564       if(FunMatrix.DomainMap().GlobalIndicesLongLong()) {
00565   Epetra_Map NewMap((long long)-1,NumMyElements,(long long)0,Comm);
00566   Epetra_Import NewImport(FunMatrix.ColMap(),NewMap);
00567   FunMatrix.ReplaceDomainMapAndImporter(NewMap,&NewImport);
00568       }
00569       else {
00570   Epetra_Map NewMap(-1,NumMyElements,0,Comm);
00571   Epetra_Import NewImport(FunMatrix.ColMap(),NewMap);
00572   FunMatrix.ReplaceDomainMapAndImporter(NewMap,&NewImport);       
00573       }
00574 
00575       // Now let's test the new importer...
00576 
00577       // Fill a random vector on the original map
00578       Epetra_Vector OriginalVec(StandardMatrix.DomainMap());
00579       Epetra_Vector OriginalY(FunMatrix.RangeMap(),true);
00580       OriginalVec.SetSeed(24601);
00581       OriginalVec.Random();
00582       
00583       // Move said random vector to a single proc
00584       Epetra_Vector NewVec(FunMatrix.DomainMap(),true);
00585       Epetra_Vector NewY(FunMatrix.RangeMap(),true);
00586       Epetra_Import ImportOld2New(FunMatrix.DomainMap(),StandardMatrix.DomainMap());
00587       NewVec.Import(OriginalVec,ImportOld2New,Add);
00588 
00589       // Test the Importer Copy Constructor
00590       Epetra_Vector ColVec1(FunMatrix.ColMap(),true);
00591       Epetra_Import ColImport(FunMatrix.ColMap(),FunMatrix.DomainMap());
00592       ColVec1.Import(NewVec,ColImport,Add);
00593 
00594       Epetra_Vector ColVec2(FunMatrix.ColMap(),true);
00595       Epetra_Import ColImport2(ColImport);
00596       ColVec2.Import(NewVec,ColImport2,Add);
00597 
00598       double norm;
00599       ColVec1.Update(-1.0,ColVec2,1.0);
00600       NewY.Norm2(&norm);
00601       if(norm > 1e-12) forierr=-1;
00602       if (verbose) {
00603   if (forierr==0) cout << "Import Copy Constructor Check OK" << endl << endl;
00604   else cout << "Import Copy Constructor Check Failed" << endl << endl;
00605       }
00606 
00607       // Test replaceDomainMapAndImporter
00608       // Now do two multiplies and compare
00609       StandardMatrix.Apply(OriginalVec,OriginalY);
00610       FunMatrix.Apply(NewVec,NewY);
00611       NewY.Update(-1.0,OriginalY,1.0);
00612       NewY.Norm2(&norm);
00613       if(norm > 1e-12) forierr=-1;
00614       if (verbose) {
00615   if (forierr==0) cout << "ReplaceDomainMapAndImporter Check OK" << endl << endl;
00616   else cout << "ReplaceDomainMapAndImporter Check Failed" << endl << endl;
00617       }
00618     }
00619   }
00620 
00621 
00622 
00623   // Release all objects
00624 
00625   delete [] SourceMyGlobalElements;
00626   delete [] TargetMyGlobalElements;
00627   delete [] OverlapMyGlobalElements;
00628   delete [] StandardMyGlobalElements;
00629 
00630   delete [] Values;
00631   delete [] Indices;
00632 
00633 #ifdef EPETRA_MPI
00634   MPI_Finalize() ;
00635 #endif
00636 
00637 /* end main
00638 */
00639 return ierr ;
00640 }
00641 
00642 int special_submap_import_test(Epetra_Comm& Comm)
00643 {
00644   int localProc = Comm.MyPID();
00645 
00646   //set up ids_source and ids_target such that ids_source are only
00647   //a subset of ids_target, and furthermore that ids_target are ordered
00648   //such that the LIDs don't match up. In other words, even if gid 2 does
00649   //exist in both ids_source and ids_target, it will correspond to different
00650   //LIDs on at least 1 proc.
00651   //
00652   //This is to test a certain bug-fix in Epetra_Import where the 'RemoteLIDs'
00653   //array wasn't being calculated correctly on all procs.
00654 
00655   int ids_source[1];
00656   ids_source[0] = localProc*2+2;
00657 
00658   int ids_target[3];
00659   ids_target[0] = localProc*2+2;
00660   ids_target[1] = localProc*2+1;
00661   ids_target[2] = localProc*2+0;
00662 
00663   Epetra_Map map_source(-1, 1, &ids_source[0], 0, Comm);
00664   Epetra_Map map_target(-1, 3, &ids_target[0], 0, Comm);
00665 
00666   Epetra_Import importer(map_target, map_source);
00667 
00668   Epetra_IntVector vec_source(map_source);
00669   Epetra_IntVector vec_target(map_target);
00670 
00671   vec_target.PutValue(0);
00672 
00673   //set vec_source's contents so that entry[i] == GID[i].
00674   int* GIDs = map_source.MyGlobalElements();
00675   for(int i=0; i<map_source.NumMyElements(); ++i) {
00676     vec_source[i] = GIDs[i];
00677   }
00678 
00679   //Import vec_source into vec_target. This should result in the contents
00680   //of vec_target remaining 0 for the entries that don't exist in vec_source,
00681   //and other entries should be equal to the corresponding GID in the map.
00682 
00683   vec_target.Import(vec_source, importer, Insert);
00684 
00685   GIDs = map_target.MyGlobalElements();
00686   int test_failed = 0;
00687 
00688   //the test passes if the i-th entry in vec_target equals either 0 or
00689   //GIDs[i].
00690   for(int i=0; i<vec_target.MyLength(); ++i) {
00691     if (vec_target[i] != GIDs[i] && vec_target[i] != 0) test_failed = 1;
00692   }
00693 
00694   int global_result;
00695   Comm.MaxAll(&test_failed, &global_result, 1);
00696 
00697   //If test didn't fail on any procs, global_result should be 0.
00698   //If test failed on any proc, global_result should be 1.
00699   return global_result;
00700 }
00701 
00702 int combine_mode_test(Epetra_Comm& Comm)
00703 {
00704   int localProc = Comm.MyPID();
00705 
00706 
00707   int ids_source[1];
00708   ids_source[0] = localProc*2+2;
00709 
00710   int ids_target[3];
00711   ids_target[0] = localProc*2+2;
00712   ids_target[1] = localProc*2+1;
00713   ids_target[2] = localProc*2+0;
00714 
00715   Epetra_Map map_source(-1, 1, &ids_source[0], 0, Comm);
00716   Epetra_Map map_target(-1, 3, &ids_target[0], 0, Comm);
00717 
00718   Epetra_Import importer(map_target, map_source);
00719 
00720   Epetra_IntVector vec_source(map_source);
00721   Epetra_IntVector vec_target(map_target);
00722 
00723   vec_target.PutValue(0);
00724 
00725   //set vec_source's contents so that entry[i] == GID[i].
00726   int* GIDs = map_source.MyGlobalElements();
00727   for(int i=0; i<map_source.NumMyElements(); ++i) {
00728     vec_source[i] = GIDs[i];
00729   }
00730 
00731   //Import vec_source into vec_target. This should result in the contents
00732   //of vec_target remaining 0 for the entries that don't exist in vec_source,
00733   //and other entries should be equal to the corresponding GID in the map.
00734 
00735   vec_target.Import(vec_source, importer, Insert);
00736 
00737   GIDs = map_target.MyGlobalElements();
00738   int test_failed = 0;
00739 
00740   //the test passes if the i-th entry in vec_target equals either 0 or
00741   //GIDs[i].
00742   for(int i=0; i<vec_target.MyLength(); ++i) {
00743     if (vec_target[i] != GIDs[i] && vec_target[i] != 0) test_failed = 1;
00744   }
00745 
00746   int global_result;
00747   Comm.MaxAll(&test_failed, &global_result, 1);
00748 
00749   //If test didn't fail on any procs, global_result should be 0.
00750   //If test failed on any proc, global_result should be 1.
00751   return global_result;
00752 }
00753 
00754 
00755 int test_import_gid(const char * name,Epetra_IntVector & Source, Epetra_IntVector & Target, const Epetra_Import & Import){
00756   int i;
00757   bool test_passed=true;
00758 
00759   // Setup
00760   for(i=0; i<Source.MyLength(); i++)
00761     Source[i] = Source.Map().GID(i);
00762   Target.PutValue(0);
00763 
00764   // Import
00765   Target.Import(Source,Import,Add);
00766 
00767   // Test
00768   for(i=0; i<Target.MyLength(); i++){
00769     if(Target[i] != Target.Map().GID(i)) test_passed=false;
00770   }
00771 
00772   if(!test_passed){
00773     printf("[%d] test_import_gid %s failed: ",Source.Map().Comm().MyPID(),name);
00774     for(i=0; i<Target.MyLength(); i++)
00775       printf("%2d(%2d) ",Target[i],Target.Map().GID(i));
00776     printf("\n");
00777     fflush(stdout);
00778   }
00779 
00780   return !test_passed;
00781 }
00782 
00783 
00784 
00785 int alternate_import_constructor_test(Epetra_Comm& Comm) {
00786   int rv=0;
00787   int nodes_per_proc=10;
00788   int numprocs = Comm.NumProc();
00789   int mypid    = Comm.MyPID();
00790 
00791   // Only run if we have multiple procs & MPI
00792   if(numprocs==0) return 0;
00793 #ifndef HAVE_MPI
00794   return 0;
00795 #endif
00796 
00797   // Build Map 1 - linear
00798   Epetra_Map Map1(-1,nodes_per_proc,0,Comm);
00799   
00800   // Build Map 2 - mod striped
00801   std::vector<int> MyGIDs(nodes_per_proc);
00802   for(int i=0; i<nodes_per_proc; i++)
00803     MyGIDs[i] = (mypid*nodes_per_proc + i) % numprocs;
00804   Epetra_Map Map2(-1,nodes_per_proc,&MyGIDs[0],0,Comm);
00805 
00806   // For testing
00807   Epetra_IntVector Source(Map1), Target(Map2);
00808 
00809 
00810   // Build Import 1 - normal
00811   Epetra_Import Import1(Map2,Map1);
00812   rv = rv|| test_import_gid("Alt test: 2 map constructor",Source,Target, Import1);
00813 
00814   // Build Import 2 - no-comm constructor
00815   int Nremote=Import1.NumRemoteIDs();
00816   const int * RemoteLIDs = Import1.RemoteLIDs();
00817   std::vector<int> RemotePIDs(Nremote+1); // I hate you, stl vector....
00818   std::vector<int> AllPIDs;  
00819   Epetra_Util::GetPids(Import1,AllPIDs,true);
00820 
00821   for(int i=0; i<Nremote; i++) {
00822     RemotePIDs[i]=AllPIDs[RemoteLIDs[i]];
00823   }
00824   Epetra_Import Import2(Import1.TargetMap(),Import1.SourceMap(),Nremote,&RemotePIDs[0],Import1.NumExportIDs(),Import1.ExportLIDs(),Import1.ExportPIDs());
00825 
00826   rv = rv || test_import_gid("Alt test: no comm constructor",Source,Target,Import2);
00827 
00828 
00829   // Build Import 3 - Remotes only
00830   Epetra_Import Import3(Import1.TargetMap(),Import1.SourceMap(),Nremote,&RemotePIDs[0]);
00831   rv = rv || test_import_gid("Alt test: remote only constructor",Source,Target, Import3);
00832 
00833 
00834   return rv;
00835 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines