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