Epetra Package Browser (Single Doxygen Collection) Development
B.cc
Go to the documentation of this file.
00001 /*
00002 //@HEADER
00003 // ************************************************************************
00004 //
00005 //               Epetra: Linear Algebra Services Package
00006 //                 Copyright 2011 Sandia Corporation
00007 //
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00039 //
00040 // ************************************************************************
00041 //@HEADER
00042 */
00043 
00044 #include <stdio.h>
00045 #include <stdlib.h>
00046 #include <ctype.h>
00047 #include <assert.h>
00048 #include <string.h>
00049 #include <math.h>
00050 #ifdef PETRA_MPI
00051 #include "mpi.h"
00052 #endif
00053 #ifndef __cplusplus
00054 #define __cplusplus
00055 #endif
00056 #include "Petra_Comm.h"
00057 #include "Petra_Map.h"
00058 #include "Petra_Time.h"
00059 #include "Petra_RDP_CRS_Matrix.h"
00060 
00061 int main(int argc, char *argv[])
00062 {
00063   int ierr = 0, i, j;
00064   bool debug = false;
00065 
00066 #ifdef PETRA_MPI
00067 
00068   // Initialize MPI
00069 
00070   MPI_Init(&argc,&argv);
00071   int size, rank; // Number of MPI processes, My process ID
00072 
00073   MPI_Comm_size(MPI_COMM_WORLD, &size);
00074   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
00075 
00076 #else
00077 
00078   int size = 1; // Serial case (not using MPI)
00079   int rank = 0;
00080 
00081 #endif
00082 
00083   bool verbose = false;
00084 
00085   // Check if we should print results to standard out
00086   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00087 
00088 
00089 
00090 #ifdef PETRA_MPI
00091   Petra_Comm & Comm = *new Petra_Comm( MPI_COMM_WORLD );
00092 #else
00093   Petra_Comm & Comm = *new Petra_Comm();
00094 #endif
00095 
00096 
00097   //char tmp;
00098   //if (rank==0) cout << "Press any key to continue..."<< endl;
00099   //if (rank==0) cin >> tmp;
00100   //Comm.Barrier();
00101 
00102   int MyPID = Comm.MyPID();
00103   int NumProc = Comm.NumProc();
00104   if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
00105               << " is alive."<<endl;
00106 
00107   bool verbose1 = verbose;
00108 
00109   // Redefine verbose to only print on PE 0
00110   if (verbose && rank!=0) verbose = false;
00111 
00112   int NumMyPoints = 10000;
00113   int NumGlobalPoints = NumMyPoints*NumProc+minfn(NumProc,3);
00114   if (MyPID < 3) NumMyPoints++;
00115   int IndexBase = 0;
00116   bool DistributedGlobal = (NumGlobalPoints>NumMyPoints);
00117 
00118   // Construct a Source Map that puts approximately the same Number of equations on each processor in
00119   // uniform global ordering
00120 
00121   Petra_Map& SourceMap = *new Petra_Map(NumGlobalPoints, NumMyPoints, 0, Comm);
00122 
00123   // Get update list and number of local equations from newly created Map
00124   int NumMyElements = SourceMap.NumMyElements();
00125   int * SourceMyGlobalElements = new int[NumMyElements];
00126   SourceMap.MyGlobalElements(SourceMyGlobalElements);
00127 
00128 
00129   // Construct a Target Map that will contain:
00130   //  some unchanged elements (relative to the soure map),
00131   //  some permuted elements
00132   //  some off-processor elements
00133   Petra_RDP_Vector & RandVec = *new Petra_RDP_Vector(SourceMap);
00134   RandVec.Random(); // This creates a vector of random numbers between negative one and one.
00135 
00136   int *TargetMyGlobalElements = new int[NumMyElements];
00137 
00138   for (i=0; i< NumMyPoints/2; i++) TargetMyGlobalElements[i] = i; // Half will be the same...
00139   for (i=NumMyPoints/2; i<NumMyPoints; i++) {
00140     int index = abs((int)(((double) (NumGlobalPoints-1) ) * RandVec[i]));
00141     TargetMyGlobalElements[i] = minfn(NumGlobalPoints-1,maxfn(0,index));
00142   }
00143 
00144   int NumSameIDs = 0;
00145   int NumPermutedIDs = 0;
00146   int NumRemoteIDs = 0;
00147   bool StillContiguous = true;
00148   for (i=0; i < NumMyPoints; i++) {
00149     if (SourceMyGlobalElements[i]==TargetMyGlobalElements[i] && StillContiguous)
00150       NumSameIDs++;
00151     else if (SourceMap.MyGID(TargetMyGlobalElements[i])) {
00152       StillContiguous = false;
00153       NumPermutedIDs++;
00154     }
00155     else {
00156       StillContiguous = false;
00157       NumRemoteIDs++;
00158     }
00159   }
00160   assert(NumMyPoints==NumSameIDs+NumPermutedIDs+NumRemoteIDs);
00161 
00162   Petra_Map & TargetMap = *new Petra_Map(-1, NumMyElements, TargetMyGlobalElements, 0, Comm);
00163 
00164   // Create a multivector whose elements are GlobalID * (column number +1)
00165 
00166   int NumVectors = 3;
00167   Petra_RDP_MultiVector & SourceMultiVector = *new Petra_RDP_MultiVector(SourceMap, NumVectors);
00168   for (j=0; j < NumVectors; j++)
00169     for (i=0; i < NumMyElements; i++)
00170       SourceMultiVector[j][i] = (double) SourceMyGlobalElements[i]*(j+1);
00171 
00172   // Create a target multivector that we will fill using an Import
00173 
00174   Petra_RDP_MultiVector & TargetMultiVector = *new Petra_RDP_MultiVector(TargetMap, NumVectors);
00175 
00176   Petra_Import & Importer = *new Petra_Import(TargetMap, SourceMap);
00177 
00178   assert(TargetMultiVector.Import(SourceMultiVector, Importer, Insert)==0);
00179 
00180   // Test Target against expected values
00181 
00182   for (j=0; j < NumVectors; j++)
00183     for (i=0; i < NumMyElements; i++) {
00184       if (TargetMultiVector[j][i]!= (double) TargetMyGlobalElements[i]*(j+1))
00185   cout << "TargetMultiVector["<<i<<"]["<<j<<"] = " << TargetMultiVector[j][i]
00186        <<  "  TargetMyGlobalElements[i]*(j+1) = " <<  TargetMyGlobalElements[i]*(j+1) << endl;
00187       assert(TargetMultiVector[j][i]== (double) TargetMyGlobalElements[i]*(j+1));
00188     }
00189 
00190   if (verbose) cout << "MultiVector Import using Importer Check OK" << endl << endl;
00191 
00192 
00194 
00195   // Now use Importer to do an export
00196 
00197   Petra_RDP_Vector & TargetVector = *new  Petra_RDP_Vector(SourceMap);
00198   Petra_RDP_Vector & ExpectedTarget = *new  Petra_RDP_Vector(SourceMap);
00199   Petra_RDP_Vector & SourceVector = *new  Petra_RDP_Vector(TargetMap);
00200 
00201   NumSameIDs = Importer.NumSameIDs();
00202   int NumPermuteIDs = Importer.NumPermuteIDs();
00203   int NumExportIDs = Importer.NumExportIDs();
00204   int *PermuteToLIDs = Importer.PermuteToLIDs();
00205   int *PermuteFromLIDs = Importer.PermuteFromLIDs();
00206   int *ExportLIDs = Importer.ExportLIDs();
00207   int *ExportPIDs = Importer.ExportPIDs();
00208 
00209   for (i=0; i < NumSameIDs; i++) ExpectedTarget[i] = (double) (MyPID+1);
00210   for (i=0; i < NumPermuteIDs; i++) ExpectedTarget[PermuteFromLIDs[i]] =
00211               (double) (MyPID+1);
00212   for (i=0; i < NumExportIDs; i++) ExpectedTarget[ExportLIDs[i]] +=
00213              (double) (ExportPIDs[i]+1);
00214 
00215   for (i=0; i < NumMyElements; i++) SourceVector[i] =  (double) (MyPID+1);
00216 
00217   assert(TargetVector.Export(SourceVector, Importer, Add)==0);
00218 
00219     for (i=0; i < NumMyElements; i++) {
00220       if (TargetVector[i]!= ExpectedTarget[i])
00221   cout <<  "     TargetVector["<<i<<"] = " << TargetVector[i]
00222        <<  "   ExpectedTarget["<<i<<"] = " <<  ExpectedTarget[i] << " on PE " << MyPID << endl;
00223       assert(TargetVector[i]== ExpectedTarget[i]);
00224     }
00225 
00226   if (verbose) cout << "Vector Export using Importer Check OK" << endl << endl;
00227 
00228 
00229 
00231   //  Build a tridiagonal system two ways:
00232   //  1) From "standard" matrix view where equations are uniquely owned.
00233   //  2) From 1D PDE view where nodes (equations) between processors are shared and partial contributions are done
00234   //     in parallel, then merged together at the end of the construction process.
00235   //
00237 
00238 
00239 
00240   // Construct a Standard Map that puts approximately the same number of equations on each processor in
00241   // uniform global ordering
00242 
00243   Petra_Map& StandardMap = *new Petra_Map(NumGlobalPoints, NumMyPoints, 0, Comm);
00244 
00245   // Get update list and number of local equations from newly created Map
00246   NumMyElements = StandardMap.NumMyElements();
00247   int * StandardMyGlobalElements = new int[NumMyElements];
00248   StandardMap.MyGlobalElements(StandardMyGlobalElements);
00249 
00250 
00251   // Create a standard Petra_CRS_Graph
00252 
00253   Petra_CRS_Graph& StandardGraph = *new Petra_CRS_Graph(Copy, StandardMap, 3);
00254   assert(!StandardGraph.IndicesAreGlobal());
00255   assert(!StandardGraph.IndicesAreLocal());
00256 
00257   // Add  rows one-at-a-time
00258   // Need some vectors to help
00259   // Off diagonal Values will always be -1
00260 
00261 
00262   int *Indices = new int[2];
00263   int NumEntries;
00264 
00265   for (i=0; i<NumMyPoints; i++)
00266     {
00267     if (StandardMyGlobalElements[i]==0)
00268       {
00269   Indices[0] = 1;
00270   NumEntries = 1;
00271       }
00272     else if (StandardMyGlobalElements[i] == NumGlobalPoints-1)
00273       {
00274   Indices[0] = NumGlobalPoints-2;
00275   NumEntries = 1;
00276       }
00277     else
00278       {
00279   Indices[0] = StandardMyGlobalElements[i]-1;
00280   Indices[1] = StandardMyGlobalElements[i]+1;
00281   NumEntries = 2;
00282       }
00283      assert(StandardGraph.InsertGlobalIndices(StandardMyGlobalElements[i], NumEntries, Indices)==0);
00284      assert(StandardGraph.InsertGlobalIndices(StandardMyGlobalElements[i], 1, StandardMyGlobalElements+i)==0); // Put in the diagonal entry
00285     }
00286 
00287   // Finish up
00288   assert(StandardGraph.IndicesAreGlobal());
00289   assert(StandardGraph.FillComplete()==0);
00290   assert(StandardGraph.IndicesAreLocal());
00291   assert(!StandardGraph.StorageOptimized());
00292   StandardGraph.OptimizeStorage();
00293   assert(StandardGraph.StorageOptimized());
00294   assert(!StandardGraph.UpperTriangular());
00295   assert(!StandardGraph.LowerTriangular());
00296 
00297 
00298   // Create Petra_RDP_CRS_Matrix using the just-built graph
00299 
00300   Petra_RDP_CRS_Matrix& StandardMatrix = *new Petra_RDP_CRS_Matrix(Copy, StandardGraph);
00301   assert(!StandardMatrix.IndicesAreGlobal());
00302   assert(StandardMatrix.IndicesAreLocal());
00303 
00304   // Add  rows one-at-a-time
00305   // Need some vectors to help
00306   // Off diagonal Values will always be -1
00307 
00308 
00309   double *Values = new double[2];
00310   Values[0] = -1.0; Values[1] = -1.0;
00311   double two = 2.0;
00312 
00313   for (i=0; i<NumMyPoints; i++)
00314     {
00315     if (StandardMyGlobalElements[i]==0)
00316       {
00317   Indices[0] = 1;
00318   NumEntries = 1;
00319       }
00320     else if (StandardMyGlobalElements[i] == NumGlobalPoints-1)
00321       {
00322   Indices[0] = NumGlobalPoints-2;
00323   NumEntries = 1;
00324       }
00325     else
00326       {
00327   Indices[0] = StandardMyGlobalElements[i]-1;
00328   Indices[1] = StandardMyGlobalElements[i]+1;
00329   NumEntries = 2;
00330       }
00331      assert(StandardMatrix.ReplaceGlobalValues(StandardMyGlobalElements[i], NumEntries, Values, Indices)==0);
00332      assert(StandardMatrix.ReplaceGlobalValues(StandardMyGlobalElements[i], 1, &two, StandardMyGlobalElements+i)==0); // Put in the diagonal entry
00333     }
00334 
00335   // Finish up
00336   assert(StandardMatrix.IndicesAreLocal());
00337   assert(StandardMatrix.FillComplete()==0);
00338   assert(StandardMatrix.IndicesAreLocal());
00339   assert(StandardMatrix.StorageOptimized());
00340   assert(!StandardMatrix.UpperTriangular());
00341   assert(!StandardMatrix.LowerTriangular());
00342 
00343   // Construct an Overlapped Map of StandardMap that include the endpoints from two neighboring processors.
00344 
00345   int OverlapNumMyElements;
00346   int OverlapMinMyGID;
00347 
00348   OverlapNumMyElements = NumMyElements + 1;
00349   if (MyPID==0) OverlapNumMyElements--;
00350 
00351   if (MyPID==0) OverlapMinMyGID = StandardMap.MinMyGID();
00352   else OverlapMinMyGID = StandardMap.MinMyGID()-1;
00353 
00354   int * OverlapMyGlobalElements = new int[OverlapNumMyElements];
00355 
00356   for (i=0; i< OverlapNumMyElements; i++) OverlapMyGlobalElements[i] = OverlapMinMyGID + i;
00357 
00358   Petra_Map& OverlapMap = *new Petra_Map(-1, OverlapNumMyElements, OverlapMyGlobalElements, 0, Comm);
00359 
00360   // Create the Overlap Petra_Matrix
00361 
00362   Petra_RDP_CRS_Matrix& OverlapMatrix = *new Petra_RDP_CRS_Matrix(Copy, OverlapMap, 4);
00363   assert(!OverlapMatrix.IndicesAreGlobal());
00364   assert(!OverlapMatrix.IndicesAreLocal());
00365 
00366   // Add  matrix element one cell at a time.
00367   // Each cell does an incoming and outgoing flux calculation
00368 
00369 
00370   double pos_one = 1.0;
00371   double neg_one = -1.0;
00372 
00373   for (i=0; i<OverlapNumMyElements; i++)
00374     {
00375       int node_left = OverlapMyGlobalElements[i]-1;
00376       int node_center = node_left + 1;
00377       int node_right = node_left + 2;
00378       if (i>0) {
00379   if (node_left>-1)
00380     assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &neg_one, &node_left)==0);
00381   assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
00382       }
00383       if (i<OverlapNumMyElements-1) {
00384   assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
00385   if (node_right<NumGlobalPoints)
00386     assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &neg_one, &node_right)==0);
00387       }
00388     }
00389 
00390   // Handle endpoints
00391   if (MyPID==0) {
00392     int node_center = 0;
00393     assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
00394   }
00395   if (MyPID==NumProc-1) {
00396     int node_center = OverlapMyGlobalElements[OverlapNumMyElements-1];
00397     assert(OverlapMatrix.InsertGlobalValues(node_center, 1, &pos_one, &node_center)==0);
00398   }
00399 
00400   assert(OverlapMatrix.FillComplete()==0);
00401 
00402   // Make a gathered matrix from OverlapMatrix.  It should be identical to StandardMatrix
00403 
00404   Petra_RDP_CRS_Matrix& GatheredMatrix = *new Petra_RDP_CRS_Matrix(Copy, StandardGraph);
00405   Petra_Export & Exporter = *new Petra_Export(OverlapMap, StandardMap);
00406   assert(GatheredMatrix.Export(OverlapMatrix, Exporter, Add)==0);
00407   assert(GatheredMatrix.FillComplete()==0);
00408 
00409   // Check if entries of StandardMatrix and GatheredMatrix are identical
00410 
00411   int StandardNumEntries, GatheredNumEntries;
00412   int * StandardIndices, * GatheredIndices;
00413   double * StandardValues, * GatheredValues;
00414 
00415   int StandardNumMyNonzeros = StandardMatrix.NumMyNonzeros();
00416   int GatheredNumMyNonzeros = GatheredMatrix.NumMyNonzeros();
00417   assert(StandardNumMyNonzeros==GatheredNumMyNonzeros);
00418 
00419   int StandardNumMyRows = StandardMatrix.NumMyRows();
00420   int GatheredNumMyRows = GatheredMatrix.NumMyRows();
00421   assert(StandardNumMyRows==GatheredNumMyRows);
00422 
00423   for (i=0; i< StandardNumMyRows; i++)
00424     {
00425       assert(StandardMatrix.ExtractMyRowView(i, StandardNumEntries, StandardValues, StandardIndices)==0);
00426       assert(GatheredMatrix.ExtractMyRowView(i, GatheredNumEntries, GatheredValues, GatheredIndices)==0);
00427       assert(StandardNumEntries==GatheredNumEntries);
00428       for (j=0; j < StandardNumEntries; j++) {
00429   //if (StandardIndices[j]!=GatheredIndices[j])
00430   // cout << "MyPID = " << MyPID << " i = " << i << "   StandardIndices[" << j << "] = " << StandardIndices[j]
00431   //      << "   GatheredIndices[" << j << "] = " << GatheredIndices[j] << endl;
00432   //if (StandardValues[j]!=GatheredValues[j])
00433   //cout << "MyPID = " << MyPID << " i = " << i << "    StandardValues[" << j << "] = " <<  StandardValues[j]
00434   //     << "    GatheredValues[" << j << "] = " <<  GatheredValues[j] << endl;
00435   assert(StandardIndices[j]==GatheredIndices[j]);
00436   assert(StandardValues[j]==GatheredValues[j]);
00437       }
00438     }
00439 
00440   if (verbose) cout << "Matrix Export Check OK" << endl;
00441   // Release all objects
00442 
00443   delete &SourceVector;
00444   delete &TargetVector;
00445   delete &ExpectedTarget;
00446 
00447 
00448   delete &Importer;
00449   delete &SourceMap;
00450   delete &TargetMap;
00451 
00452   delete [] SourceMyGlobalElements;
00453   delete [] TargetMyGlobalElements;
00454 
00455   delete &SourceMultiVector;
00456   delete &TargetMultiVector;
00457   delete &RandVec;
00458 
00459   delete &Exporter;
00460   delete &GatheredMatrix;
00461   delete &OverlapMatrix;
00462   delete &OverlapMap;
00463   delete [] OverlapMyGlobalElements;
00464 
00465   delete &StandardMatrix;
00466   delete &StandardGraph;
00467   delete &StandardMap;
00468   delete [] StandardMyGlobalElements;
00469 
00470   delete [] Values;
00471   delete [] Indices;
00472   delete &Comm;
00473 
00474 #ifdef PETRA_MPI
00475   MPI_Finalize() ;
00476 #endif
00477 
00478 /* end main
00479 */
00480 return ierr ;
00481 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines