00001 // 00002 // C++ Implementation: directsolveroperator 00003 // 00004 // Description: Epetra_Operator wrapper for direct solvers (e.g. Amesos solvers) 00005 // 00006 // 00007 // Author: Roman Geus <geus@maxwell>, (C) 2004 00008 // 00009 // Copyright: See COPYING file that comes with this distribution 00010 // 00011 // 00012 00013 #include <Epetra_Map.h> 00014 #include <Epetra_Operator.h> 00015 #include <Epetra_MultiVector.h> 00016 #include <Epetra_LinearProblem.h> 00017 #include <Teuchos_ParameterList.hpp> 00018 #include <Amesos.h> 00019 00020 #include "myrlog.h" 00021 00022 #include "directsolveroperator.h" 00023 00024 DirectSolverOperator::DirectSolverOperator(Epetra_RowMatrix& A, Teuchos::ParameterList& params) 00025 : domain_map_(A.OperatorDomainMap()), 00026 range_map_(A.OperatorRangeMap()), 00027 comm_(A.Comm()), 00028 use_transpose_(false) 00029 { 00030 char* solver_name = const_cast<char *>(params.get<string>("solver_type").c_str()); 00031 linear_problem_ = new Epetra_LinearProblem(); 00032 linear_problem_->SetOperator(&A); 00033 Amesos factory; 00034 if (!factory.Query(solver_name)) 00035 rExit("Unknown direct solver: %s", solver_name); 00036 00037 amesos_solver_ = factory.Create(solver_name, *linear_problem_); 00038 assert(amesos_solver_ != 0); 00039 int ierr = amesos_solver_->SetParameters(params); 00040 assert(ierr == 0); 00041 ierr = amesos_solver_->SymbolicFactorization(); 00042 assert(ierr == 0); 00043 ierr = amesos_solver_->NumericFactorization(); 00044 assert(ierr == 0); 00045 } 00046 00047 DirectSolverOperator::~DirectSolverOperator() { 00048 delete amesos_solver_; 00049 delete linear_problem_; 00050 } 00051 00052 int DirectSolverOperator::SetUseTranspose(bool use_transpose) { 00053 return amesos_solver_->SetUseTranspose(use_transpose); 00054 } 00055 00056 int DirectSolverOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const { 00057 return 1; 00058 } 00059 00060 int DirectSolverOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const { 00061 int ierr; 00062 00063 if (!X.Map().SameAs(OperatorDomainMap())) 00064 EPETRA_CHK_ERR(-1); 00065 if (!Y.Map().SameAs(OperatorRangeMap())) 00066 EPETRA_CHK_ERR(-2); 00067 if (Y.NumVectors() != X.NumVectors()) 00068 EPETRA_CHK_ERR(-3); 00069 if (!OperatorDomainMap().SameAs(OperatorRangeMap())) 00070 EPETRA_CHK_ERR(-4); 00071 00072 linear_problem_->SetLHS(&Y); 00073 // FIXME: theoretically X could be modified (e.g. scaled) by the solver. Should we operate on a copy? 00074 linear_problem_->SetRHS(&const_cast<Epetra_MultiVector&>(X)); 00075 EPETRA_CHK_ERR(ierr = amesos_solver_->Solve()); 00076 00077 return ierr; 00078 } 00079 00080 double DirectSolverOperator::NormInf() const { 00081 return 0.0; 00082 } 00083 00084 const char * DirectSolverOperator::Label() const { 00085 return "DirectSolverOperator"; 00086 } 00087 00088 bool DirectSolverOperator::UseTranspose() const { 00089 return amesos_solver_->UseTranspose(); 00090 } 00091 00092 bool DirectSolverOperator::HasNormInf() const { 00093 return false; 00094 } 00095 00096 const Epetra_Comm& DirectSolverOperator::Comm() const { 00097 return comm_; 00098 } 00099 00100 const Epetra_Map & DirectSolverOperator::OperatorDomainMap() const { 00101 return domain_map_; 00102 } 00103 00104 const Epetra_Map & DirectSolverOperator::OperatorRangeMap() const { 00105 return range_map_; 00106 }