test/LinearProblemRedistor/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_LinearProblemRedistor Test routine
00030 
00031 //
00032 //  Limitations
00033 //  -----------
00034 //
00035 //  We test only square matrices.  Indeed we perform this test only on 
00036 //  a single square (but non-symmetric) matrix.
00037 //
00038 //  We test only one type of redistor:  the one needed by 
00039 //  SuperLUdist.  i.e. ConstructTranspose = true -and- 
00040 //  MakeDataContiguous = true.  
00041 //
00042 //  We exercise only one constructor, again the constructor required by 
00043 //  SuperLUdist.  i.e. OrigProblem, NumProc, Replicate with 
00044 //  NumProc = comm.NumProc() and Replicate = true.  
00045 //
00046 //  Tests performed
00047 //  ---------------
00048 //
00049 //  We create three Epetra_LinearProblems and redistributions:
00050 //   1)  Create an Epetra_LinearProblem and redistribute it
00051 //   2)  Update the right hand side of an Epetra_LinearProblem and 
00052 //       update the redistributed version of the right hand side
00053 //   3)  Update the matrix of an Epetra_LinearProblem and 
00054 //       update the redistributed version of the matrix
00055 //  For each of these we compare the linear problem and its
00056 //  redistribution as follows:
00057 //
00058 //  Method
00059 //  ------
00060 //
00061 //  We compare a linear problem and the redistributed version of same
00062 //  by multiplying the matrix in the linear problem and the redistributed
00063 //  version of the linear problem by the right hand side to yield a 
00064 //  left hand side which we can then compare.
00065 //
00066 //  We perform a matrix vector multiply instead of a solve to avoid 
00067 //  dependence on any particular solver.  However, this results in 
00068 //  computing x = Ab instead of solving Ax =b.  This only makes sense 
00069 //  if A is square.
00070 //
00071 
00072 
00073 #ifdef EPETRA_MPI
00074 #include "Epetra_MpiComm.h"
00075 #include <mpi.h>
00076 #endif
00077 #include "Epetra_CrsMatrix.h"
00078 #include "Epetra_VbrMatrix.h"
00079 #include "Epetra_Vector.h"
00080 #include "Epetra_MultiVector.h"
00081 #include "Epetra_LocalMap.h"
00082 #include "Epetra_IntVector.h"
00083 #include "Epetra_Map.h"
00084 
00085 #include "Epetra_SerialComm.h"
00086 #include "Epetra_Time.h"
00087 #include "Epetra_LinearProblem.h"
00088 #include "Epetra_LinearProblemRedistor.h"
00089 #include "Trilinos_Util.h"
00090 #ifdef HAVE_COMM_ASSERT_EQUAL
00091 #include "Comm_assert_equal.h"
00092 #endif
00093 
00094 int checkResults( bool trans, 
00095       Epetra_LinearProblemRedistor * redistor, 
00096       Epetra_LinearProblem * OrigProb,
00097       Epetra_LinearProblem * RedistProb,
00098       bool verbose) ;
00099   
00100 int main(int argc, char *argv[]) {
00101 
00102   int i;
00103 
00104 #ifdef EPETRA_MPI
00105   // Initialize MPI
00106   MPI_Init(&argc,&argv);
00107   Epetra_MpiComm comm(MPI_COMM_WORLD);
00108 #else
00109   Epetra_SerialComm comm;
00110 #endif
00111 
00112   // Uncomment to debug in parallel int tmp; if (comm.MyPID()==0) cin >> tmp; comm.Barrier();
00113 
00114   bool verbose = false;
00115   bool veryVerbose = false;
00116 
00117   // Check if we should print results to standard out
00118   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00119 
00120   if (!verbose) comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00121 
00122   if (verbose) cout << comm << endl << flush;
00123 
00124   bool verbose1 = verbose;
00125   if (verbose) verbose = (comm.MyPID()==0);
00126 
00127   int nx = 4;
00128   int ny = comm.NumProc()*nx; // Scale y grid with number of processors
00129   
00130   // Create funky stencil to make sure the matrix is non-symmetric (transpose non-trivial):
00131 
00132   // (i-1,j-1) (i-1,j  )
00133   // (i  ,j-1) (i  ,j  ) (i  ,j+1)
00134   // (i+1,j-1) (i+1,j  )
00135   
00136   int npoints = 2;
00137 
00138   int xoff[] = {-1,  0,  1, -1,  0,  1,  0};
00139   int yoff[] = {-1, -1, -1,  0,  0,  0,  1};
00140 
00141   Epetra_Map * map;
00142   Epetra_CrsMatrix * A;
00143   Epetra_Vector * x, * b, * xexact;
00144   
00145   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
00146 
00147   if (verbose) 
00148     cout << "npoints = " << npoints << " nx = " << nx << " ny = " << ny  << endl ; 
00149 
00150   if (verbose && nx<6 ) {
00151     cout << *A << endl;
00152     cout << "B       = " << endl << *b << endl;
00153   }
00154   // Construct linear problem object
00155   
00156   Epetra_LinearProblem origProblem(A, x, b);
00157   Epetra_LinearProblem *redistProblem;
00158   
00159   Epetra_Time timer(comm);
00160   
00161   // Construct redistor object, use all processors and replicate full problem on each
00162 
00163   double start = timer.ElapsedTime();
00164   Epetra_LinearProblemRedistor redistor(&origProblem, comm.NumProc(), true);
00165   if (verbose) cout << "\nTime to construct redistor  = " 
00166         << timer.ElapsedTime() - start << endl;
00167 
00168   bool ConstructTranspose = true; 
00169   bool MakeDataContiguous = true;
00170   
00171   start = timer.ElapsedTime();
00172   redistor.CreateRedistProblem(ConstructTranspose, MakeDataContiguous, redistProblem);
00173   if (verbose) cout << "\nTime to create redistributed problem = " 
00174         << timer.ElapsedTime() - start << endl;
00175   
00176   
00177   // Now test output of redistor by performing matvecs
00178 
00179   int ierr = 0;
00180   ierr += checkResults( ConstructTranspose, &redistor, &origProblem, 
00181       redistProblem, verbose);
00182   
00183   
00184   // Now change values in original rhs and test update facility of redistor
00185   // Multiply b by 2
00186   
00187   double Value = 2.0;
00188   
00189   b->Scale(Value); // b = 2*b
00190   
00191   redistor.UpdateRedistRHS(b);
00192   if (verbose) cout << "\nTime to update redistributed RHS  = " 
00193         << timer.ElapsedTime() - start << endl;
00194   
00195   ierr += checkResults( ConstructTranspose, &redistor, 
00196       &origProblem, redistProblem, verbose);
00197  
00198   // Now change values in original matrix and test update facility of redistor
00199 
00200 #define  CREATE_CONST_MATRIX
00201 #ifdef CREATE_CONST_MATRIX
00202   //  The easiest way that I could find to change the matrix without EPETRA_CHK_ERRs
00203   A->PutScalar(13.0); 
00204 #else 
00205 
00206   //  This has no effect on matrices, such as when nx = 4, that have no 
00207   //  diagonal entries.  However, it does cause many EPETRA_CHK_ERR prints.
00208 
00209   // Add 2 to the diagonal of each row 
00210   for (i=0; i< A->NumMyRows(); i++)  {
00211   //  for (i=0; i < 1; i++)
00212     cout << " i = " << i ; 
00213     A->SumIntoMyValues(i, 1, &Value, &i);
00214   }
00215 #endif  
00216   
00217   
00218   start = timer.ElapsedTime();
00219   redistor.UpdateRedistProblemValues(&origProblem);
00220   if (verbose) cout << "\nTime to update redistributed problem  = " 
00221         << timer.ElapsedTime() - start << endl;
00222   
00223   ierr += checkResults(ConstructTranspose, &redistor, &origProblem, redistProblem, verbose);
00224   
00225   delete A;
00226   delete b;
00227   delete x;
00228   delete xexact;
00229   delete map;
00230   
00231   
00232 #ifdef EPETRA_MPI
00233   MPI_Finalize();
00234 #endif
00235 
00236   return ierr;
00237 }
00238 
00239 //
00240 //  checkResults computes Ax = b1 and either Rx=b2 or R^T x=b2 (depending on trans) where
00241 //  A = origProblem and R = redistProblem.  
00242 //  checkResults returns 0 (OK) if norm(b1-b2)/norm(b1) < size(b1) * 1e-15; 
00243 //
00244 //  I guess that we need the redistor so that we can move x and b2 around.  
00245 //
00246 const double error_tolerance = 1e-15 ; 
00247 
00248 int checkResults( bool trans, 
00249       Epetra_LinearProblemRedistor * redistor, 
00250       Epetra_LinearProblem * A,
00251       Epetra_LinearProblem * R,
00252       bool verbose) {
00253   
00254   int m = A->GetRHS()->MyLength();
00255   int n = A->GetLHS()->MyLength();
00256   assert( m == n ) ; 
00257 
00258   Epetra_MultiVector *x = A->GetLHS() ; 
00259   Epetra_MultiVector x1( *x ) ; 
00260   //  Epetra_MultiVector Difference( x1 ) ; 
00261   Epetra_MultiVector *b = A->GetRHS();
00262   Epetra_RowMatrix *matrixA = A->GetMatrix();
00263   assert( matrixA != 0 ) ; 
00264   int iam = matrixA->Comm().MyPID();
00265   
00266   //  Epetra_Time timer(A->Comm());
00267   //  double start = timer.ElapsedTime();
00268 
00269   matrixA->Multiply(trans, *b, x1) ;   // x = Ab 
00270 
00271   int M,N,nz;
00272   int *ptr, *ind;
00273   double *val, *rhs, *lhs;
00274   int Nrhs, ldrhs, ldlhs;
00275 
00276   redistor->ExtractHbData( M, N, nz, ptr, ind, 
00277         val, Nrhs, rhs, ldrhs, 
00278         lhs, ldlhs);
00279 
00280   assert( M == N ) ; 
00281   if ( verbose ) {
00282     cout << " iam = " << iam 
00283    << " m = " << m  << " n = " << n  << " M = " << M << endl ;  
00284     
00285     cout << " iam = " << iam << " ptr = " << ptr[0] << " "   << ptr[1] << " "  << ptr[2] << " "  << ptr[3] << " "  << ptr[4] << " " << ptr[5] << endl ;
00286     
00287     cout << " iam = " << iam << " ind = " << ind[0] << " "   << ind[1] << " "  << ind[2] << " "  << ind[3] << " "  << ind[4] << " " << ind[5] << endl ;
00288     
00289     cout << " iam = " << iam << " val = " << val[0] << " "   << val[1] << " "  << val[2] << " "  << val[3] << " "  << val[4] << " " << val[5] << endl ;
00290   }
00291   //  Create a serial map in case we end up needing it 
00292   //  If it is created inside the else block below it would have to
00293   //  be with a call to new().
00294   int NumMyElements_ = 0 ;
00295   if (matrixA->Comm().MyPID()==0) NumMyElements_ = n;
00296   Epetra_Map SerialMap( n, NumMyElements_, 0, matrixA->Comm() );
00297 
00298   //  These are unnecessary and useless
00299   //  Epetra_Vector serial_A_rhs( SerialMap ) ; 
00300   //  Epetra_Vector serial_A_lhs( SerialMap ) ;
00301 
00302   //  Epetra_Export exporter( matrixA->BlockRowMap(), SerialMap) ;
00303 
00304   //
00305   //  In each process, we will compute Rb putting the answer into LHS
00306   //
00307 
00308 
00309   for ( int k = 0 ; k < Nrhs; k ++ ) { 
00310     for ( int i = 0 ; i < M ; i ++ ) { 
00311       lhs[ i + k * ldlhs ] = 0.0; 
00312     }
00313     for ( int i = 0 ; i < M ; i++ ) { 
00314       for ( int l = ptr[i]; l < ptr[i+1]; l++ ) {
00315   int j = ind[l] ; 
00316   if ( verbose && N < 40 ) {
00317     cout << " i = " << i << " j = " << j ;
00318     cout << " l = " << l << " val[l] = " << val[l] ;
00319     cout << " rhs = " << rhs[ j + k * ldrhs ] << endl ;
00320   }
00321   lhs[ i + k * ldrhs ] += val[l] * rhs[ j + k * ldrhs ] ;
00322       }
00323     }
00324 
00325     if ( verbose && N < 40 ) { 
00326       cout << " lhs = " ; 
00327       for ( int j = 0 ; j < N ; j++ ) cout << " " << lhs[j] ; 
00328       cout << endl ; 
00329       cout << " rhs = " ; 
00330       for ( int j = 0 ; j < N ; j++ ) cout << " " << rhs[j] ; 
00331       cout << endl ; 
00332     }
00333 
00334     const Epetra_Comm &comm = matrixA->Comm() ; 
00335 #ifdef HAVE_COMM_ASSERT_EQUAL
00336     //
00337     //  Here we double check to make sure that lhs and rhs are 
00338     //  replicated.  
00339     //
00340     for ( int j = 0 ; j < N ; j++ ) { 
00341       assert( Comm_assert_equal( &comm, lhs[ j + k * ldrhs ] ) ) ; 
00342       assert( Comm_assert_equal( &comm, rhs[ j + k * ldrhs ] ) ) ; 
00343     } 
00344 #endif
00345   }
00346 
00347   //
00348   //  Now we have to redistribue them back 
00349   //
00350   redistor->UpdateOriginalLHS( A->GetLHS() ) ; 
00351   //
00352   //  Now we want to compare x and x1 which have been computed as follows:
00353   //  x = Rb  
00354   //  x1 = Ab 
00355   //
00356 
00357   double Norm_x1, Norm_diff ;
00358   EPETRA_CHK_ERR( x1.Norm2( &Norm_x1 ) ) ; 
00359 
00360   //  cout << " x1 = " << x1 << endl ; 
00361   //  cout << " *x = " << *x << endl ; 
00362 
00363   x1.Update( -1.0, *x, 1.0 ) ; 
00364   EPETRA_CHK_ERR( x1.Norm2( &Norm_diff ) ) ; 
00365 
00366   //  cout << " diff, i.e. updated x1 = " << x1 << endl ; 
00367 
00368   int ierr = 0;
00369 
00370   if ( verbose ) {
00371     cout << " Norm_diff = "  << Norm_diff << endl ; 
00372     cout << " Norm_x1 = "  << Norm_x1 << endl ; 
00373   }
00374 
00375   if ( Norm_diff / Norm_x1 > n * error_tolerance ) ierr++ ; 
00376 
00377   if (ierr!=0 && verbose) cerr << "Status: Test failed" << endl;
00378   else if (verbose) cerr << "Status: Test passed" << endl;
00379 
00380   return(ierr);
00381 }

Generated on Wed May 12 21:41:04 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7