4 #define OUR_CHK_ERR(a) { { int epetra_err = a; \ 5 if (epetra_err != 0) { std::cerr << "Amesos ERROR " << epetra_err << ", " \ 6 << __FILE__ << ", line " << __LINE__ << std::endl; \ 7 relresidual=1e15; return(1);} }\ 11 #include "Epetra_Comm.h" 14 #include "Epetra_CrsMatrix.h" 15 #include "Epetra_Map.h" 16 #include "Epetra_Vector.h" 17 #include "Epetra_LinearProblem.h" 27 const Epetra_Comm &Comm,
31 Epetra_CrsMatrix *& Amat,
38 int iam = Comm.MyPID() ;
41 const Epetra_Map *RangeMap =
42 transpose?&Amat->OperatorDomainMap():&Amat->OperatorRangeMap() ;
43 const Epetra_Map *DomainMap =
44 transpose?&Amat->OperatorRangeMap():&Amat->OperatorDomainMap() ;
46 Epetra_Vector xexact(*DomainMap);
47 Epetra_Vector x(*DomainMap);
49 Epetra_Vector b(*RangeMap);
50 Epetra_Vector bcheck(*RangeMap);
52 Epetra_Vector difference(*DomainMap);
54 Epetra_LinearProblem Problem;
58 Abase = Afactory.
Create( AmesosClass, Problem ) ;
60 std::string AC = AmesosClass ;
61 if ( AC ==
"Amesos_Mumps" ) {
62 ParamList.
set(
"NoDestroy",
true );
66 double relresidual = 0 ;
72 Problem.SetOperator( Amat );
86 xexact.PutScalar(1.0);
91 Amat->Multiply( transpose, xexact, b ) ;
94 std::cout << __FILE__ <<
"::" << __LINE__ <<
"b = " << std::endl ;
95 b.Print( std::cout ) ;
96 std::cout << __FILE__ <<
"::" << __LINE__ <<
"xexact = " << std::endl ;
97 xexact.Print( std::cout ) ;
98 std::cout << __FILE__ <<
"::" << __LINE__ <<
"x = " << std::endl ;
99 x.Print( std::cout ) ;
107 Problem.SetLHS( &x );
108 Problem.SetRHS( &b );
122 Amat->Multiply( transpose, x, bcheck ) ;
127 difference.Update( 1.0, x, -1.0, xexact, 0.0 ) ;
128 difference.Norm2( &norm_diff ) ;
129 x.Norm2( &norm_one ) ;
131 relresidual = norm_diff / norm_one ;
134 if ( relresidual * Rcond > 1e-16 ) {
135 if (
verbose) std::cout << __FILE__ <<
"::"<< __LINE__
136 <<
" norm( x - xexact ) / norm(x) = " 137 << norm_diff /norm_one << std::endl ;
157 const Epetra_Comm &Comm,
161 Epetra_CrsMatrix *& Amat,
165 std::string AC = AmesosClass ;
167 if (epetra_err != 0) {
168 std::cerr <<
"Amesos ERROR " << epetra_err <<
", " << __FILE__ <<
", line " << __LINE__ << std::endl;
virtual int Solve()=0
Solves A X = B (or AT x = B)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
virtual int NumericFactorization()=0
Performs NumericFactorization on the matrix A.
virtual int SymbolicFactorization()=0
Performs SymbolicFactorization on the matrix A.
virtual int SetParameters(Teuchos::ParameterList &ParameterList)=0
Updates internal variables.
int PartialFactorization(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond)
Factory for binding a third party direct solver to an Epetra_LinearProblem.
virtual int SetUseTranspose(bool UseTranspose)=0
If set true, X will be set to the solution of AT X = B (not A X = B)
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Amesos Create method.
int PartialFactorizationOneStep(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond, int Steps)
Amesos_BaseSolver: A pure virtual class for direct solution of real-valued double-precision operators...