src/directsolveroperator.cpp

Go to the documentation of this file.
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 }

Generated on Fri Oct 26 13:35:11 2007 for FEMAXX (Finite Element Maxwell Eigensolver) by  doxygen 1.4.7