test/RowMatrixTransposer/cxx_main.cpp

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

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