Epetra Package Browser (Single Doxygen Collection) Development
test/RowMatrixTransposer/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_LocalMap.h"
00045 #include "Epetra_Time.h"
00046 #include "Epetra_CrsMatrix.h"
00047 #include "Epetra_VbrMatrix.h"
00048 #include "Epetra_Vector.h"
00049 #include "Epetra_Flops.h"
00050 #ifdef EPETRA_MPI
00051 #include "Epetra_MpiComm.h"
00052 #include "mpi.h"
00053 #else
00054 #include "Epetra_SerialComm.h"
00055 #endif
00056 #include "../epetra_test_err.h"
00057 #include "Epetra_IntVector.h"
00058 #include "Epetra_Version.h"
00059 #include "Epetra_RowMatrixTransposer.h"
00060 #include "Epetra_Time.h"
00061 
00062 int checkResults(Epetra_RowMatrix * A,
00063                  Epetra_CrsMatrix * transA,
00064                  Epetra_MultiVector * xexact,
00065                  bool verbose);
00066 
00067 void GenerateCrsProblem(int nx, int ny, int npoints,
00068                         int * xoff, int * yoff, int nrhs,
00069                         const Epetra_Comm  &comm,
00070                         Epetra_Map *& map,
00071                         Epetra_CrsMatrix *& A,
00072                         Epetra_MultiVector *& x,
00073                         Epetra_MultiVector *& b,
00074                         Epetra_MultiVector *&xexact);
00075 
00076 void GenerateVbrProblem(int nx, int ny, int npoints, int * xoff, int * yoff,
00077                         int nsizes, int * sizes, int nrhs,
00078                         const Epetra_Comm  &comm,
00079                         Epetra_BlockMap *& map,
00080                         Epetra_VbrMatrix *& A,
00081                         Epetra_MultiVector *& x,
00082                         Epetra_MultiVector *& b,
00083                         Epetra_MultiVector *&xexact);
00084 
00085 int main(int argc, char *argv[]) {
00086 
00087   int ierr = 0, i;
00088 
00089 #ifdef EPETRA_MPI
00090 
00091   // Initialize MPI
00092 
00093   MPI_Init(&argc,&argv);
00094 
00095   Epetra_MpiComm Comm( MPI_COMM_WORLD );
00096 
00097 #else
00098 
00099   Epetra_SerialComm Comm;
00100 
00101 #endif
00102 
00103   bool verbose = false;
00104 
00105   // Check if we should print results to standard out
00106   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00107 
00108   int verbose_int = verbose ? 1 : 0;
00109   Comm.Broadcast(&verbose_int, 1, 0);
00110   verbose = verbose_int==1 ? true : false;
00111 
00112 
00113   //  char tmp;
00114   //  if (Comm.MyPID()==0) cout << "Press any key to continue..."<< endl;
00115   //  if (Comm.MyPID()==0) cin >> tmp;
00116   //  Comm.Barrier();
00117 
00118   Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00119   int MyPID = Comm.MyPID();
00120   int NumProc = Comm.NumProc();
00121 
00122   if(verbose && MyPID==0)
00123     cout << Epetra_Version() << endl << endl;
00124 
00125   int nx = 128;
00126   int ny = NumProc*nx; // Scale y grid with number of processors
00127 
00128   // Create funky stencil to make sure the matrix is non-symmetric (transpose non-trivial):
00129 
00130   // (i-1,j-1) (i-1,j  )
00131   // (i  ,j-1) (i  ,j  ) (i  ,j+1)
00132   // (i+1,j-1) (i+1,j  )
00133 
00134   int npoints = 7;
00135 
00136   int xoff[] = {-1,  0,  1, -1,  0,  1,  0};
00137   int yoff[] = {-1, -1, -1,  0,  0,  0,  1};
00138 
00139   Epetra_Map * map;
00140   Epetra_CrsMatrix * A;
00141   Epetra_MultiVector * x, * b, * xexact;
00142 
00143   GenerateCrsProblem(nx, ny, npoints, xoff, yoff, 1, Comm, map, A, x, b, xexact);
00144 
00145   if (nx<8) {
00146     cout << *A << endl;
00147     cout << "X exact = " << endl << *xexact << endl;
00148     cout << "B       = " << endl << *b << endl;
00149   }
00150   // Construct transposer object
00151 
00152   Epetra_Time timer(Comm);
00153 
00154   double start = timer.ElapsedTime();
00155   Epetra_RowMatrixTransposer transposer(A);
00156   if (verbose) cout << "\nTime to construct transposer  = "
00157        << timer.ElapsedTime() - start << endl;
00158 
00159   bool MakeDataContiguous = true;
00160   Epetra_CrsMatrix * transA;
00161 
00162   start = timer.ElapsedTime();
00163   transposer.CreateTranspose(MakeDataContiguous, transA);
00164   if (verbose) cout << "\nTime to create transpose matrix  = "
00165        << timer.ElapsedTime() - start << endl;
00166 
00167 
00168   // Now test output of transposer by performing matvecs
00169   ierr += checkResults(A, transA, xexact, verbose);
00170 
00171 
00172   // Now change values in original matrix and test update facility of transposer
00173   // Add 2 to the diagonal of each row
00174 
00175   double Value = 2.0;
00176   for (i=0; i< A->NumMyRows(); i++)
00177     A->SumIntoMyValues(i, 1, &Value, &i);
00178 
00179 
00180   start = timer.ElapsedTime();
00181   transposer.UpdateTransposeValues(A);
00182   if (verbose) cout << "\nTime to update transpose matrix  = "
00183      << timer.ElapsedTime() - start << endl;
00184 
00185   ierr += checkResults(A, transA, xexact, verbose);
00186 
00187   delete A;
00188   delete b;
00189   delete x;
00190   delete xexact;
00191   delete map;
00192 
00193 
00194   if (verbose) cout << endl << "Checking transposer for VbrMatrix objects" << endl<< endl;
00195 
00196   int nsizes = 4;
00197   int sizes[] = {4, 6, 5, 3};
00198 
00199   Epetra_VbrMatrix * Avbr;
00200   Epetra_BlockMap * bmap;
00201 
00202   GenerateVbrProblem(nx, ny, npoints, xoff, yoff, nsizes, sizes,
00203                      1,  Comm, bmap, Avbr, x, b, xexact);
00204 
00205   if (nx<8) {
00206     cout << *Avbr << endl;
00207     cout << "X exact = " << endl << *xexact << endl;
00208     cout << "B       = " << endl << *b << endl;
00209   }
00210 
00211   start = timer.ElapsedTime();
00212   Epetra_RowMatrixTransposer transposer1(Avbr);
00213   if (verbose) cout << "\nTime to construct transposer  = "
00214        << timer.ElapsedTime() - start << endl;
00215 
00216   start = timer.ElapsedTime();
00217   transposer1.CreateTranspose(MakeDataContiguous, transA);
00218   if (verbose) cout << "\nTime to create transpose matrix  = "
00219        << timer.ElapsedTime() - start << endl;
00220 
00221 
00222   // Now test output of transposer by performing matvecs
00223   ierr += checkResults(Avbr, transA, xexact, verbose);
00224 
00225   // Now change values in original matrix and test update facility of transposer
00226   // Scale matrix on the left by rowsums
00227 
00228   Epetra_Vector invRowSums(Avbr->RowMap());
00229 
00230   Avbr->InvRowSums(invRowSums);
00231   Avbr->LeftScale(invRowSums);
00232 
00233   start = timer.ElapsedTime();
00234   transposer1.UpdateTransposeValues(Avbr);
00235   if (verbose) cout << "\nTime to update transpose matrix  = "
00236     << timer.ElapsedTime() - start << endl;
00237 
00238   ierr += checkResults(Avbr, transA, xexact, verbose);
00239 
00240   delete Avbr;
00241   delete b;
00242   delete x;
00243   delete xexact;
00244   delete bmap;
00245 
00246 #ifdef EPETRA_MPI
00247   MPI_Finalize() ;
00248 #endif
00249 
00250 /* end main */
00251   return ierr ;
00252 }
00253 
00254 int checkResults(Epetra_RowMatrix * A, Epetra_CrsMatrix * transA,
00255                  Epetra_MultiVector * xexact, bool verbose)
00256 {
00257   int n = A->NumGlobalRows();
00258 
00259   if (n<100) cout << "A transpose = " << endl << *transA << endl;
00260 
00261   Epetra_MultiVector x1(View, A->OperatorRangeMap(), &((*xexact)[0]), 1);
00262   Epetra_MultiVector b1(A->OperatorDomainMap(), 1);
00263 
00264   A->SetUseTranspose(true);
00265 
00266   Epetra_Time timer(A->Comm());
00267   double start = timer.ElapsedTime();
00268   A->Apply(x1, b1);
00269   if (verbose) cout << "\nTime to compute b1: matvec with original matrix using transpose flag  = " << timer.ElapsedTime() - start << endl;
00270 
00271   if (n<100) cout << "b1 = " << endl << b1 << endl;
00272   Epetra_MultiVector x2(View, transA->OperatorDomainMap(), &((*xexact)[0]), 1);
00273   Epetra_MultiVector b2(transA->OperatorRangeMap(), 1);
00274   start = timer.ElapsedTime();
00275   transA->Multiply(false, x2, b2);
00276   if (verbose) cout << "\nTime to compute b2: matvec with transpose matrix                      = " << timer.ElapsedTime() - start << endl;
00277 
00278   if (n<100) cout << "b1 = " << endl << b1 << endl;
00279 
00280   double residual;
00281   Epetra_MultiVector resid(A->OperatorRangeMap(), 1);
00282 
00283   resid.Update(1.0, b1, -1.0, b2, 0.0);
00284   int ierr0 = resid.Norm2(&residual);
00285   assert(ierr0==0);
00286   if (verbose) cout << "Norm of b1 - b2 = " << residual << endl;
00287 
00288   int ierr = 0;
00289 
00290   if (residual > 1.0e-10) ierr++;
00291 
00292   if (ierr!=0 && verbose) {
00293     cerr << "Status: Test failed" << endl;
00294   }
00295   else if (verbose) cerr << "Status: Test passed" << endl;
00296 
00297   return(ierr);
00298 }
00299 
00300 void GenerateCrsProblem(int nx, int ny, int npoints,
00301                         int * xoff, int * yoff, int nrhs,
00302                         const Epetra_Comm  &comm,
00303                         Epetra_Map *& map,
00304                         Epetra_CrsMatrix *& A,
00305                         Epetra_MultiVector *& x,
00306                         Epetra_MultiVector *& b,
00307                         Epetra_MultiVector *&xexact)
00308 {
00309   // Number of global equations is nx*ny.  These will be distributed in a linear fashion
00310   int numGlobalEquations = nx*ny;
00311   map = new Epetra_Map(numGlobalEquations, 0, comm); // Create map with equal distribution of equations.
00312 
00313   int numMyEquations = map->NumMyElements();
00314 
00315   A = new Epetra_CrsMatrix(Copy, *map, 0); // Construct matrix
00316 
00317   int * indices = new int[npoints];
00318   double * values = new double[npoints];
00319 
00320   double dnpoints = (double) npoints;
00321 
00322   for (int i=0; i<numMyEquations; i++) {
00323 
00324     int rowID = map->GID(i);
00325     int numIndices = 0;
00326 
00327     for (int j=0; j<npoints; j++) {
00328       int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
00329       if (colID>-1 && colID<numGlobalEquations) {
00330         indices[numIndices] = colID;
00331         double value = - ((double) rand())/ ((double) RAND_MAX);
00332         if (colID==rowID)
00333           values[numIndices++] = dnpoints - value; // Make diagonal dominant
00334         else
00335           values[numIndices++] = -value;
00336       }
00337     }
00338 
00339     A->InsertGlobalValues(rowID, numIndices, values, indices);
00340   }
00341 
00342   delete [] indices;
00343   delete [] values;
00344 
00345   A->FillComplete();
00346 
00347   if (nrhs<=1) {
00348     x = new Epetra_Vector(*map);
00349     b = new Epetra_Vector(*map);
00350     xexact = new Epetra_Vector(*map);
00351   }
00352   else {
00353     x = new Epetra_MultiVector(*map, nrhs);
00354     b = new Epetra_MultiVector(*map, nrhs);
00355     xexact = new Epetra_MultiVector(*map, nrhs);
00356   }
00357 
00358   xexact->Random(); // Fill xexact with random values
00359 
00360   A->Multiply(false, *xexact, *b);
00361 
00362   return;
00363 }
00364 
00365 void GenerateVbrProblem(int nx, int ny, int npoints, int * xoff, int * yoff,
00366                         int nsizes, int * sizes, int nrhs,
00367                         const Epetra_Comm  &comm,
00368                         Epetra_BlockMap *& map,
00369                         Epetra_VbrMatrix *& A,
00370                         Epetra_MultiVector *& x,
00371                         Epetra_MultiVector *& b,
00372                         Epetra_MultiVector *&xexact)
00373 {
00374   int i, j;
00375 
00376   // Number of global equations is nx*ny.  These will be distributed in a linear fashion
00377   int numGlobalEquations = nx*ny;
00378   Epetra_Map ptMap(numGlobalEquations, 0, comm); // Create map with equal distribution of equations.
00379 
00380   int numMyElements = ptMap.NumMyElements();
00381 
00382   Epetra_IntVector elementSizes(ptMap); // This vector will have the list of element sizes
00383   for (i=0; i<numMyElements; i++)
00384     elementSizes[i] = sizes[ptMap.GID(i)%nsizes]; // cycle through sizes array
00385 
00386   map = new Epetra_BlockMap(-1, numMyElements, ptMap.MyGlobalElements(), elementSizes.Values(), ptMap.IndexBase(), ptMap.Comm());
00387 
00388 
00389   A = new Epetra_VbrMatrix(Copy, *map, 0); // Construct matrix
00390 
00391   int * indices = new int[npoints];
00392 
00393   // This section of code creates a vector of random values that will be used to create
00394   // light-weight dense matrices to pass into the VbrMatrix construction process.
00395 
00396   int maxElementSize = 0;
00397   for (i=0; i< nsizes; i++) maxElementSize = EPETRA_MAX(maxElementSize, sizes[i]);
00398 
00399   Epetra_LocalMap lmap(maxElementSize*maxElementSize, ptMap.IndexBase(), ptMap.Comm());
00400   Epetra_Vector randvec(lmap);
00401   randvec.Random();
00402   randvec.Scale(-1.0); // Make value negative
00403 
00404 
00405   for (i=0; i<numMyElements; i++) {
00406     int rowID = map->GID(i);
00407     int numIndices = 0;
00408     int rowDim = sizes[rowID%nsizes];
00409     for (j=0; j<npoints; j++) {
00410       int colID = rowID + xoff[j] + nx*yoff[j]; // Compute column ID based on stencil offsets
00411       if (colID>-1 && colID<numGlobalEquations)
00412               indices[numIndices++] = colID;
00413     }
00414 
00415     A->BeginInsertGlobalValues(rowID, numIndices, indices);
00416 
00417     for (j=0; j < numIndices; j++) {
00418       int colDim = sizes[indices[j]%nsizes];
00419       A->SubmitBlockEntry(&(randvec[0]), rowDim, rowDim, colDim);
00420     }
00421     A->EndSubmitEntries();
00422   }
00423 
00424   delete [] indices;
00425 
00426   A->FillComplete();
00427 
00428   // Compute the InvRowSums of the matrix rows
00429   Epetra_Vector invRowSums(A->RowMap());
00430   Epetra_Vector rowSums(A->RowMap());
00431   A->InvRowSums(invRowSums);
00432   rowSums.Reciprocal(invRowSums);
00433 
00434   // Jam the row sum values into the diagonal of the Vbr matrix (to make it diag dominant)
00435   int numBlockDiagonalEntries;
00436   int * rowColDims;
00437   int * diagoffsets = map->FirstPointInElementList();
00438   A->BeginExtractBlockDiagonalView(numBlockDiagonalEntries, rowColDims);
00439   for (i=0; i< numBlockDiagonalEntries; i++) {
00440     double * diagVals;
00441     int diagLDA;
00442     A->ExtractBlockDiagonalEntryView(diagVals, diagLDA);
00443     int rowDim = map->ElementSize(i);
00444     for (j=0; j<rowDim; j++) diagVals[j+j*diagLDA] = rowSums[diagoffsets[i]+j];
00445   }
00446 
00447   if (nrhs<=1) {
00448     x = new Epetra_Vector(*map);
00449     b = new Epetra_Vector(*map);
00450     xexact = new Epetra_Vector(*map);
00451   }
00452   else {
00453     x = new Epetra_MultiVector(*map, nrhs);
00454     b = new Epetra_MultiVector(*map, nrhs);
00455     xexact = new Epetra_MultiVector(*map, nrhs);
00456   }
00457 
00458   xexact->Random(); // Fill xexact with random values
00459 
00460   A->Multiply(false, *xexact, *b);
00461 
00462   return;
00463 }
00464 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines