27 const Epetra_Comm &Comm,
30 Teuchos::ParameterList ParamList,
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 );
78 if ( transpose )
OUR_CHK_ERR( Abase->SetUseTranspose( transpose ) );
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,
160 Teuchos::ParameterList ParamList,
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;
int PartialFactorizationOneStep(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond, int Steps)
int PartialFactorization(const char *AmesosClass, const Epetra_Comm &Comm, bool transpose, bool verbose, Teuchos::ParameterList ParamList, Epetra_CrsMatrix *&Amat, double Rcond)