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

Generated on Thu Sep 18 12:37:56 2008 for Epetra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1