Amesos Package Browser (Single Doxygen Collection) Development
PartialFactorization.cpp
Go to the documentation of this file.
00001 //
00002 //  OUR_CHK_ERR always returns 1 on error.
00003 //
00004 #define OUR_CHK_ERR(a) { { int epetra_err = a; \
00005                       if (epetra_err != 0) { std::cerr << "Amesos ERROR " << epetra_err << ", " \
00006                            << __FILE__ << ", line " << __LINE__ << std::endl; \
00007 relerror = 1.3e15; relresidual=1e15; return(1);}  }\
00008                    }
00009 
00010 
00011 #include "Epetra_Comm.h"
00012 #include "Teuchos_ParameterList.hpp"
00013 #include "Amesos.h"
00014 #include "Epetra_CrsMatrix.h"
00015 #include "Epetra_Map.h"
00016 #include "Epetra_Vector.h"
00017 #include "Epetra_LinearProblem.h"
00018 #include "PerformOneSolveAndTest.h"
00019 #include "CreateTridi.h"
00020 //
00021 //  PartialFactorization checks to make sure that we clean up after our mess
00022 //  before everything is done.
00023 //  
00024 
00025 const int MaxNumSteps = 7 ; 
00026 
00027 int PartialFactorizationOneStep( const char* AmesosClass,
00028          const Epetra_Comm &Comm, 
00029          bool transpose, 
00030          bool verbose, 
00031          Teuchos::ParameterList ParamList, 
00032          Epetra_CrsMatrix *& Amat, 
00033          double Rcond, 
00034          int Steps ) 
00035 {
00036   
00037   assert( Steps >= 0 && Steps < MaxNumSteps ) ; 
00038 
00039   int iam = Comm.MyPID() ; 
00040   int errors = 0 ; 
00041 
00042   const Epetra_Map *RangeMap = 
00043     transpose?&Amat->OperatorDomainMap():&Amat->OperatorRangeMap() ; 
00044   const Epetra_Map *DomainMap = 
00045     transpose?&Amat->OperatorRangeMap():&Amat->OperatorDomainMap() ; 
00046 
00047   Epetra_Vector xexact(*DomainMap);
00048   Epetra_Vector x(*DomainMap);
00049 
00050   Epetra_Vector b(*RangeMap);
00051   Epetra_Vector bcheck(*RangeMap);
00052 
00053   Epetra_Vector difference(*DomainMap);
00054 
00055   Epetra_LinearProblem Problem;
00056   Amesos_BaseSolver* Abase ; 
00057   Amesos Afactory;
00058 
00059   Abase = Afactory.Create( AmesosClass, Problem ) ; 
00060 
00061   std::string AC = AmesosClass ;
00062   if ( AC == "Amesos_Mumps" ) { 
00063     ParamList.set( "NoDestroy", true );
00064    Abase->SetParameters( ParamList ) ; 
00065   }
00066 
00067   double relerror = 0 ; 
00068   double relresidual = 0 ; 
00069   
00070   if ( Steps > 0 ) {
00071     //
00072     //  Phase 1:  Compute b = A' A' A xexact
00073     //
00074     Problem.SetOperator( Amat );
00075    
00076     //
00077     //  We only set transpose if we have to - this allows valgrind to check
00078     //  that transpose is set to a default value before it is used.
00079     //
00080     if ( transpose ) OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) ); 
00081     //    if (verbose) ParamList.set( "DebugLevel", 1 );
00082     //    if (verbose) ParamList.set( "OutputLevel", 1 );
00083     if ( Steps > 1 ) {
00084       OUR_CHK_ERR( Abase->SetParameters( ParamList ) ); 
00085       if ( Steps > 2 ) {
00086     
00087   int ind[1];
00088   ind[0] = 0;
00089   xexact.Random();
00090   xexact.PutScalar(1.0);
00091   
00092   //
00093   //  Compute cAx = A' xexact
00094   //
00095   Amat->Multiply( transpose, xexact, b ) ;  //  b = A x2 = A A' A'' xexact
00096 
00097 #if 0 
00098   std::cout << __FILE__ << "::"  << __LINE__ << "b = " << std::endl ; 
00099   b.Print( std::cout ) ; 
00100   std::cout << __FILE__ << "::"  << __LINE__ << "xexact = " << std::endl ; 
00101   xexact.Print( std::cout ) ; 
00102   std::cout << __FILE__ << "::"  << __LINE__ << "x = " << std::endl ; 
00103   x.Print( std::cout ) ; 
00104 #endif
00105   //
00106   //  Phase 2:  Solve A' A' A x = b 
00107   //
00108   //
00109   //  Solve A sAAx = b 
00110   //
00111   Problem.SetLHS( &x );
00112   Problem.SetRHS( &b );
00113   OUR_CHK_ERR( Abase->SymbolicFactorization(  ) ); 
00114   if ( Steps > 2 ) {
00115     OUR_CHK_ERR( Abase->SymbolicFactorization(  ) ); 
00116     if ( Steps > 3 ) {
00117       OUR_CHK_ERR( Abase->NumericFactorization(  ) ); 
00118       if ( Steps > 4 ) {
00119         OUR_CHK_ERR( Abase->NumericFactorization(  ) ); 
00120         if ( Steps > 5 ) {
00121     OUR_CHK_ERR( Abase->Solve(  ) ); 
00122     if ( Steps > 6 ) {
00123       OUR_CHK_ERR( Abase->Solve(  ) ); 
00124 
00125 
00126       Amat->Multiply( transpose, x, bcheck ) ; //  temp = A" x2
00127       
00128       double norm_diff ;
00129       double norm_one ;
00130       
00131       difference.Update( 1.0, x, -1.0, xexact, 0.0 ) ;
00132       difference.Norm2( &norm_diff ) ; 
00133       x.Norm2( &norm_one ) ; 
00134       
00135       relerror = norm_diff / norm_one ; 
00136       
00137       relresidual = norm_diff / norm_one ; 
00138       
00139       if (iam == 0 ) {
00140         if ( relresidual * Rcond > 1e-16 ) {
00141           if (verbose) std::cout << __FILE__ << "::"<< __LINE__ 
00142           << " norm( x - xexact ) / norm(x) = " 
00143           << norm_diff /norm_one << std::endl ; 
00144           errors += 1 ; 
00145         }
00146       }
00147     }
00148         }
00149       }
00150     }
00151   }
00152       }
00153     }
00154 }
00155  delete Abase;
00156  
00157  return errors;
00158  
00159 }
00160 
00161 
00162 int PartialFactorization( const char* AmesosClass,
00163         const Epetra_Comm &Comm, 
00164         bool transpose, 
00165         bool verbose, 
00166         Teuchos::ParameterList ParamList, 
00167         Epetra_CrsMatrix *& Amat, 
00168         double Rcond ) {
00169 
00170   double relerror = 0 ; 
00171   double relresidual = 0 ; 
00172 
00173   for( int i =0 ; i < MaxNumSteps ; i ++ ) {
00174     std::string AC = AmesosClass ; 
00175 
00176     //
00177     //  I have turned this test back on because it no longer seems to fail.  Although Amesos_Dscpack
00178     //  fails other aspects of TestOptions.
00179     //
00180     //  Amesos_Dscpack dies if SymbolicFactorization() is called twice in a row - Bug #1237 
00181     //  Hence we do not test Amesos_Dscpack with 3 or more steps here.
00182     //    if ( AC != "Amesos_Dscpack" || i < 3 ) { 
00183     {
00184       OUR_CHK_ERR( PartialFactorizationOneStep( AmesosClass,
00185                Comm, 
00186                transpose, 
00187                verbose, 
00188                ParamList, 
00189                Amat, 
00190                Rcond,
00191                i ) ) ;
00192     }
00193   }
00194   return(0);
00195 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines