Epetra Package Browser (Single Doxygen Collection) Development
bug1_main.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //               Epetra: 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 
00043 //  Epetra_LinearProblemRedistor Test routine
00044 
00045 //
00046 //  This file demonstrates a problem in some destructor somewhere - or
00047 //  perhaps the proble could be in this test code.
00048 //
00049 //  In any case, the symptom is a Segmentation fualt
00050 //
00051 
00052 
00053 #ifdef EPETRA_MPI
00054 #include "Epetra_MpiComm.h"
00055 #include <mpi.h>
00056 #endif
00057 #include "Epetra_CrsMatrix.h"
00058 #include "Epetra_VbrMatrix.h"
00059 #include "Epetra_Vector.h"
00060 #include "Epetra_MultiVector.h"
00061 #include "Epetra_LocalMap.h"
00062 #include "Epetra_IntVector.h"
00063 #include "Epetra_Map.h"
00064 
00065 #include "Epetra_SerialComm.h"
00066 #include "Epetra_Time.h"
00067 #include "Epetra_LinearProblem.h"
00068 #include "Epetra_LinearProblemRedistor.h"
00069 #include "Trilinos_Util.h"
00070 #ifdef HAVE_COMM_ASSERT_EQUAL
00071 #include "Comm_assert_equal.h"
00072 #endif
00073 
00074 int checkResults( bool trans,
00075       Epetra_LinearProblemRedistor * redistor,
00076       Epetra_LinearProblem * OrigProb,
00077       Epetra_LinearProblem * RedistProb,
00078       bool verbose) ;
00079 
00080 int main(int argc, char *argv[]) {
00081 
00082 {
00083   int i;
00084 
00085 #ifdef EPETRA_MPI
00086   // Initialize MPI
00087   MPI_Init(&argc,&argv);
00088   Epetra_MpiComm comm(MPI_COMM_WORLD);
00089 #else
00090   Epetra_SerialComm comm;
00091 #endif
00092 
00093   // Uncomment to debug in parallel int tmp; if (comm.MyPID()==0) cin >> tmp; comm.Barrier();
00094 
00095   bool verbose = false;
00096   bool veryVerbose = false;
00097 
00098   // Check if we should print results to standard out
00099   if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;
00100 
00101   if (!verbose) comm.SetTracebackMode(0); // This should shut down any error traceback reporting
00102 
00103   if (verbose) cout << comm << endl << flush;
00104 
00105   bool verbose1 = verbose;
00106   if (verbose) verbose = (comm.MyPID()==0);
00107 
00108   int nx = 4;
00109   int ny = comm.NumProc()*nx; // Scale y grid with number of processors
00110 
00111   // Create funky stencil to make sure the matrix is non-symmetric (transpose non-trivial):
00112 
00113   // (i-1,j-1) (i-1,j  )
00114   // (i  ,j-1) (i  ,j  ) (i  ,j+1)
00115   // (i+1,j-1) (i+1,j  )
00116 
00117   int npoints = 2;
00118 
00119   int xoff[] = {-1,  0,  1, -1,  0,  1,  0};
00120   int yoff[] = {-1, -1, -1,  0,  0,  0,  1};
00121 
00122   Epetra_Map * map;
00123   Epetra_CrsMatrix * A;
00124   Epetra_Vector * x, * b, * xexact;
00125   
00126   Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);
00127 
00128   if (verbose)
00129     cout << "npoints = " << npoints << " nx = " << nx << " ny = " << ny  << endl ;
00130 
00131   if (verbose && nx<6 ) {
00132     cout << *A << endl;
00133     cout << "B       = " << endl << *b << endl;
00134   }
00135   // Construct linear problem object
00136 
00137   Epetra_LinearProblem origProblem(A, x, b);
00138   Epetra_LinearProblem *redistProblem;
00139 
00140   Epetra_Time timer(comm);
00141 
00142   // Construct redistor object, use all processors and replicate full problem on each
00143 
00144   double start = timer.ElapsedTime();
00145   Epetra_LinearProblemRedistor redistor(&origProblem, comm.NumProc(), true);
00146   if (verbose) cout << "\nTime to construct redistor  = "
00147         << timer.ElapsedTime() - start << endl;
00148 
00149   bool ConstructTranspose = true;
00150   bool MakeDataContiguous = true;
00151 
00152   //
00153   //  Commenting out the following define avoids the segmentation fault
00154   //
00155 #define CREATEREDISTPROBLEM
00156 #ifdef CREATEREDISTPROBLEM
00157   start = timer.ElapsedTime();
00158   redistor.CreateRedistProblem(ConstructTranspose, MakeDataContiguous, redistProblem);
00159   if (verbose) cout << "\nTime to create redistributed problem = "
00160         << timer.ElapsedTime() - start << endl;
00161 
00162 #endif
00163 #if 0
00164   // Now test output of redistor by performing matvecs
00165 
00166   int ierr = 0;
00167   ierr += checkResults( ConstructTranspose, &redistor, &origProblem,
00168       redistProblem, verbose);
00169 #endif
00170   cout << " Here we are just before the return()  "  << endl ;
00171 }
00172  return(0) ;
00173 
00174 }
00175 
00176 //
00177 //  checkResults computes Ax = b1 and either Rx=b2 or R^T x=b2 (depending on trans) where
00178 //  A = origProblem and R = redistProblem.
00179 //  checkResults returns 0 (OK) if norm(b1-b2)/norm(b1) < size(b1) * 1e-15;
00180 //
00181 //  I guess that we need the redistor so that we can move x and b2 around.
00182 //
00183 const double error_tolerance = 1e-15 ;
00184 
00185 int checkResults( bool trans,
00186       Epetra_LinearProblemRedistor * redistor,
00187       Epetra_LinearProblem * A,
00188       Epetra_LinearProblem * R,
00189       bool verbose) {
00190 
00191   int m = A->GetRHS()->MyLength();
00192   int n = A->GetLHS()->MyLength();
00193   assert( m == n ) ;
00194 
00195   Epetra_MultiVector *x = A->GetLHS() ;
00196   Epetra_MultiVector x1( *x ) ;
00197   //  Epetra_MultiVector Difference( x1 ) ;
00198   Epetra_MultiVector *b = A->GetRHS();
00199   Epetra_RowMatrix *matrixA = A->GetMatrix();
00200   assert( matrixA != 0 ) ;
00201   int iam = matrixA->Comm().MyPID();
00202 
00203   //  Epetra_Time timer(A->Comm());
00204   //  double start = timer.ElapsedTime();
00205 
00206   matrixA->Multiply(trans, *b, x1) ;   // x = Ab
00207 
00208   int M,N,nz;
00209   int *ptr, *ind;
00210   double *val, *rhs, *lhs;
00211   int Nrhs, ldrhs, ldlhs;
00212 
00213   redistor->ExtractHbData( M, N, nz, ptr, ind,
00214         val, Nrhs, rhs, ldrhs,
00215         lhs, ldlhs);
00216 
00217   assert( M == N ) ;
00218   if ( verbose ) {
00219     cout << " iam = " << iam
00220    << " m = " << m  << " n = " << n  << " M = " << M << endl ;
00221 
00222     cout << " iam = " << iam << " ptr = " << ptr[0] << " "   << ptr[1] << " "  << ptr[2] << " "  << ptr[3] << " "  << ptr[4] << " " << ptr[5] << endl ;
00223 
00224     cout << " iam = " << iam << " ind = " << ind[0] << " "   << ind[1] << " "  << ind[2] << " "  << ind[3] << " "  << ind[4] << " " << ind[5] << endl ;
00225 
00226     cout << " iam = " << iam << " val = " << val[0] << " "   << val[1] << " "  << val[2] << " "  << val[3] << " "  << val[4] << " " << val[5] << endl ;
00227   }
00228   //  Create a serial map in case we end up needing it
00229   //  If it is created inside the else block below it would have to
00230   //  be with a call to new().
00231   int NumMyElements_ = 0 ;
00232   if (matrixA->Comm().MyPID()==0) NumMyElements_ = n;
00233   Epetra_Map SerialMap( n, NumMyElements_, 0, matrixA->Comm() );
00234 
00235   //  These are unnecessary and useless
00236   //  Epetra_Vector serial_A_rhs( SerialMap ) ;
00237   //  Epetra_Vector serial_A_lhs( SerialMap ) ;
00238 
00239   //  Epetra_Export exporter( matrixA->BlockRowMap(), SerialMap) ;
00240 
00241   //
00242   //  In each process, we will compute Rb putting the answer into LHS
00243   //
00244 
00245 
00246   for ( int k = 0 ; k < Nrhs; k ++ ) {
00247     for ( int i = 0 ; i < M ; i ++ ) {
00248       lhs[ i + k * ldlhs ] = 0.0;
00249     }
00250     for ( int i = 0 ; i < M ; i++ ) {
00251       for ( int l = ptr[i]; l < ptr[i+1]; l++ ) {
00252   int j = ind[l] ;
00253   if ( verbose && N < 40 ) {
00254     cout << " i = " << i << " j = " << j ;
00255     cout << " l = " << l << " val[l] = " << val[l] ;
00256     cout << " rhs = " << rhs[ j + k * ldrhs ] << endl ;
00257   }
00258   lhs[ i + k * ldrhs ] += val[l] * rhs[ j + k * ldrhs ] ;
00259       }
00260     }
00261 
00262     if ( verbose && N < 40 ) {
00263       cout << " lhs = " ;
00264       for ( int j = 0 ; j < N ; j++ ) cout << " " << lhs[j] ;
00265       cout << endl ;
00266       cout << " rhs = " ;
00267       for ( int j = 0 ; j < N ; j++ ) cout << " " << rhs[j] ;
00268       cout << endl ;
00269     }
00270 
00271     const Epetra_Comm &comm = matrixA->Comm() ;
00272 #ifdef HAVE_COMM_ASSERT_EQUAL
00273     //
00274     //  Here we double check to make sure that lhs and rhs are
00275     //  replicated.
00276     //
00277     for ( int j = 0 ; j < N ; j++ ) {
00278       assert( Comm_assert_equal( &comm, lhs[ j + k * ldrhs ] ) ) ;
00279       assert( Comm_assert_equal( &comm, rhs[ j + k * ldrhs ] ) ) ;
00280     }
00281 #endif
00282   }
00283 
00284   //
00285   //  Now we have to redistribue them back
00286   //
00287   redistor->UpdateOriginalLHS( A->GetLHS() ) ;
00288   //
00289   //  Now we want to compare x and x1 which have been computed as follows:
00290   //  x = Rb
00291   //  x1 = Ab
00292   //
00293 
00294   double Norm_x1, Norm_diff ;
00295   EPETRA_CHK_ERR( x1.Norm2( &Norm_x1 ) ) ;
00296 
00297   //  cout << " x1 = " << x1 << endl ;
00298   //  cout << " *x = " << *x << endl ;
00299 
00300   x1.Update( -1.0, *x, 1.0 ) ;
00301   EPETRA_CHK_ERR( x1.Norm2( &Norm_diff ) ) ;
00302 
00303   //  cout << " diff, i.e. updated x1 = " << x1 << endl ;
00304 
00305   int ierr = 0;
00306 
00307   if ( verbose ) {
00308     cout << " Norm_diff = "  << Norm_diff << endl ;
00309     cout << " Norm_x1 = "  << Norm_x1 << endl ;
00310   }
00311 
00312   if ( Norm_diff / Norm_x1 > n * error_tolerance ) ierr++ ;
00313 
00314   if (ierr!=0 && verbose) cerr << "Status: Test failed" << endl;
00315   else if (verbose) cerr << "Status: Test passed" << endl;
00316 
00317   return(ierr);
00318 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines