EpetraExt Package Browser (Single Doxygen Collection) Development
test/Transpose/cxx_main.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - 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 // Epetra_BlockMap Test routine
00043 
00044 #include "EpetraExt_Version.h"
00045 
00046 #include "Epetra_CrsMatrix.h"
00047 #include "Epetra_VbrMatrix.h"
00048 #include "Epetra_Vector.h"
00049 #include "Epetra_MultiVector.h"
00050 #include "Epetra_LocalMap.h"
00051 #include "Epetra_IntVector.h"
00052 #include "Epetra_Map.h"
00053 
00054 #ifdef EPETRA_MPI
00055 #include "Epetra_MpiComm.h"
00056 #include <mpi.h>
00057 #else
00058 #include "Epetra_SerialComm.h"
00059 #endif
00060 
00061 #include "Epetra_Time.h"
00062 #include "EpetraExt_Transpose_RowMatrix.h"
00063 #include "Trilinos_Util.h"
00064 
00065 int checkResults(Epetra_RowMatrix * A, Epetra_CrsMatrix * transA, 
00066                  Epetra_Vector * xexact, bool verbose);
00067 
00068 int main(int argc, char *argv[])
00069 {
00070   int i;
00071 
00072 #ifdef EPETRA_MPI
00073   // Initialize MPI
00074   MPI_Init(&argc,&argv);
00075   Epetra_MpiComm comm(MPI_COMM_WORLD);
00076 #else
00077   Epetra_SerialComm comm;
00078 #endif
00079 
00080   // Uncomment to debug in parallel int tmp; if (comm.MyPID()==0) cin >> tmp; comm.Barrier();
00081 
00082   bool verbose = false;
00083 
00084   // Check if we should print results to standard out
00085   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00086 
00087   if (!verbose) comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00088 
00089   if (verbose) std::cout << comm << std::endl << std::flush;
00090 
00091   if (verbose) verbose = (comm.MyPID()==0);
00092 
00093   if (verbose)
00094     std::cout << EpetraExt::EpetraExt_Version() << std::endl << std::endl;
00095 
00096   int nx = 128;
00097   int ny = comm.NumProc()*nx; // Scale y grid with number of processors
00098 
00099   // Create funky stencil to make sure the matrix is non-symmetric (transpose non-trivial):
00100 
00101   // (i-1,j-1) (i-1,j  )
00102   // (i  ,j-1) (i  ,j  ) (i  ,j+1)
00103   // (i+1,j-1) (i+1,j  )
00104 
00105   int npoints = 7;
00106 
00107   int xoff[] = {-1,  0,  1, -1,  0,  1,  0};
00108   int yoff[] = {-1, -1, -1,  0,  0,  0,  1};
00109 
00110   Epetra_Map * map;
00111   Epetra_CrsMatrix * A;
00112   Epetra_Vector * x, * b, * xexact;
00113   
00114   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
00115 
00116   if (nx<8)
00117   {
00118     std::cout << *A << std::endl;
00119     std::cout << "X exact = " << std::endl << *xexact << std::endl;
00120     std::cout << "B       = " << std::endl << *b << std::endl;
00121   }
00122 
00123   // Construct transposer 
00124   Epetra_Time timer(comm);
00125 
00126   double start = timer.ElapsedTime();
00127 
00128   //bool IgnoreNonLocalCols = false;
00129   EpetraExt::RowMatrix_Transpose transposer;
00130 
00131   if (verbose) std::cout << "\nTime to construct transposer  = " << timer.ElapsedTime() - start << std::endl;
00132   
00133   Epetra_CrsMatrix & transA = dynamic_cast<Epetra_CrsMatrix&>(transposer(*A));
00134 
00135   start = timer.ElapsedTime();
00136   if (verbose) std::cout << "\nTime to create transpose matrix  = " << timer.ElapsedTime() - start << std::endl;
00137   
00138   // Now test output of transposer by performing matvecs
00139   int ierr = 0;
00140   ierr += checkResults(A, &transA, xexact, verbose);
00141 
00142 
00143   // Now change values in original matrix and test update facility of transposer
00144   // Add 2 to the diagonal of each row
00145   double Value = 2.0;
00146   for (i=0; i< A->NumMyRows(); i++)
00147   A->SumIntoMyValues(i, 1, &Value, &i);
00148 
00149   start = timer.ElapsedTime();
00150   transposer.fwd();
00151 
00152   if (verbose) std::cout << "\nTime to update transpose matrix  = " << timer.ElapsedTime() - start << std::endl;
00153   
00154   ierr += checkResults(A, &transA, xexact, verbose);
00155 
00156   delete A;
00157   delete b;
00158   delete x;
00159   delete xexact;
00160   delete map;
00161 
00162   if (verbose) std::cout << std::endl << "Checking transposer for VbrMatrix objects" << std::endl<< std::endl;
00163 
00164   int nsizes = 4;
00165   int sizes[] = {4, 6, 5, 3};
00166 
00167   Epetra_VbrMatrix * Avbr;
00168   Epetra_BlockMap * bmap;
00169 
00170   Trilinos_Util_GenerateVbrProblem(nx, ny, npoints, xoff, yoff, nsizes, sizes,
00171                                    comm, bmap, Avbr, x, b, xexact);
00172 
00173   if (nx<8)
00174   {
00175     std::cout << *Avbr << std::endl;
00176     std::cout << "X exact = " << std::endl << *xexact << std::endl;
00177     std::cout << "B       = " << std::endl << *b << std::endl;
00178   }
00179 
00180   start = timer.ElapsedTime();
00181   EpetraExt::RowMatrix_Transpose transposer1;
00182 
00183   Epetra_CrsMatrix & transA1 = dynamic_cast<Epetra_CrsMatrix&>(transposer1(*Avbr));
00184   if (verbose) std::cout << "\nTime to create transpose matrix  = " << timer.ElapsedTime() - start << std::endl;
00185   
00186   // Now test output of transposer by performing matvecs
00187 ;
00188   ierr += checkResults(Avbr, &transA1, xexact, verbose);
00189 
00190   // Now change values in original matrix and test update facility of transposer
00191   // Scale matrix on the left by rowsums
00192 
00193   Epetra_Vector invRowSums(Avbr->RowMap());
00194 
00195   Avbr->InvRowSums(invRowSums);
00196   Avbr->LeftScale(invRowSums);
00197 
00198   start = timer.ElapsedTime();
00199   transposer1.fwd();
00200   if (verbose) std::cout << "\nTime to update transpose matrix  = " << timer.ElapsedTime() - start << std::endl;
00201   
00202   ierr += checkResults(Avbr, &transA1, xexact, verbose);
00203 
00204   delete Avbr;
00205   delete b;
00206   delete x;
00207   delete xexact;
00208   delete bmap;
00209 
00210 #ifdef EPETRA_MPI
00211   MPI_Finalize();
00212 #endif
00213 
00214   return ierr;
00215 }
00216 
00217 int checkResults(Epetra_RowMatrix * A, Epetra_CrsMatrix * transA, 
00218                  Epetra_Vector * xexact, bool verbose) {
00219 
00220   int n = A->NumGlobalRows();
00221 
00222   if (n<100) std::cout << "A transpose = " << std::endl << *transA << std::endl;
00223 
00224   Epetra_Vector x1(View,A->OperatorDomainMap(), &((*xexact)[0]));
00225   Epetra_Vector b1(A->OperatorRangeMap());
00226 
00227   A->SetUseTranspose(true);
00228 
00229   Epetra_Time timer(A->Comm());
00230   double start = timer.ElapsedTime();
00231   A->Apply(x1, b1);
00232   if (verbose) std::cout << "\nTime to compute b1: matvec with original matrix using transpose flag  = " << timer.ElapsedTime() - start << std::endl;
00233 
00234   if (n<100) std::cout << "b1 = " << std::endl << b1 << std::endl;
00235   Epetra_Vector x2(View,transA->OperatorRangeMap(), &((*xexact)[0]));
00236   Epetra_Vector b2(transA->OperatorDomainMap());
00237   start = timer.ElapsedTime();
00238   transA->Multiply(false, x2, b2);
00239   if (verbose) std::cout << "\nTime to compute b2: matvec with transpose matrix                      = " << timer.ElapsedTime() - start << std::endl;
00240 
00241   if (n<100) std::cout << "b1 = " << std::endl << b1 << std::endl;
00242 
00243   double residual;
00244   Epetra_Vector resid(A->OperatorDomainMap());
00245 
00246   resid.Update(1.0, b1, -1.0, b2, 0.0);
00247   resid.Norm2(&residual);
00248   if (verbose) std::cout << "Norm of b1 - b2 = " << residual << std::endl;
00249 
00250   int ierr = 0;
00251 
00252   if (residual > 1.0e-10) ierr++;
00253 
00254   if (ierr!=0 && verbose) std::cerr << "Status: Test failed" << std::endl;
00255   else if (verbose) std::cerr << "Status: Test passed" << std::endl;
00256 
00257   return(ierr);
00258 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines