B.cc

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

Generated on Thu Sep 18 12:37:56 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1