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