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