src/directsolveroperator.cpp

Go to the documentation of this file.

//
// C++ Implementation: directsolveroperator
//
// Description: Epetra_Operator wrapper for direct solvers (e.g. Amesos solvers)
//
//
// Author: Roman Geus <geus@maxwell>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#include <Epetra_Map.h>
#include <Epetra_Operator.h>
#include <Epetra_MultiVector.h>
#include <Epetra_LinearProblem.h>
#include <Teuchos_ParameterList.hpp>
#include <Amesos.h>

#include "myrlog.h"

#include "directsolveroperator.h"

DirectSolverOperator::DirectSolverOperator(Epetra_RowMatrix& A, Teuchos::ParameterList& params)
      : domain_map_(A.OperatorDomainMap()),
        range_map_(A.OperatorRangeMap()),
        comm_(A.Comm()),
        use_transpose_(false)
{
    char* solver_name = const_cast<char *>(params.get<string>("solver_type").c_str());
    linear_problem_ = new Epetra_LinearProblem();
    linear_problem_->SetOperator(&A);
    Amesos factory;
    if (!factory.Query(solver_name))
        rExit("Unknown direct solver: %s", solver_name);

    amesos_solver_ = factory.Create(solver_name, *linear_problem_);
    assert(amesos_solver_ != 0);
    int ierr = amesos_solver_->SetParameters(params);
    assert(ierr == 0);
    ierr = amesos_solver_->SymbolicFactorization();
    assert(ierr == 0);
    ierr = amesos_solver_->NumericFactorization();
    assert(ierr == 0);
}

DirectSolverOperator::~DirectSolverOperator() {
    delete amesos_solver_;
    delete linear_problem_;
}

int DirectSolverOperator::SetUseTranspose(bool use_transpose) {
    return amesos_solver_->SetUseTranspose(use_transpose);
}

int DirectSolverOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    return 1;
}

int DirectSolverOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    int ierr;

    if (!X.Map().SameAs(OperatorDomainMap()))
        EPETRA_CHK_ERR(-1);
    if (!Y.Map().SameAs(OperatorRangeMap()))
        EPETRA_CHK_ERR(-2);
    if (Y.NumVectors() != X.NumVectors())
        EPETRA_CHK_ERR(-3);
    if (!OperatorDomainMap().SameAs(OperatorRangeMap()))
        EPETRA_CHK_ERR(-4);

    linear_problem_->SetLHS(&Y);
    // FIXME: theoretically X could be modified (e.g. scaled) by the solver. Should we operate on a copy?
    linear_problem_->SetRHS(&const_cast<Epetra_MultiVector&>(X));
    EPETRA_CHK_ERR(ierr = amesos_solver_->Solve());

    return ierr;
}

double DirectSolverOperator::NormInf() const {
    return 0.0;
}

const char * DirectSolverOperator::Label() const {
    return "DirectSolverOperator";
}

bool DirectSolverOperator::UseTranspose() const {
    return amesos_solver_->UseTranspose();
}

bool DirectSolverOperator::HasNormInf() const {
    return false;
}

const Epetra_Comm& DirectSolverOperator::Comm() const {
    return comm_;
}

const Epetra_Map & DirectSolverOperator::OperatorDomainMap() const {
    return domain_map_;
}

const Epetra_Map & DirectSolverOperator::OperatorRangeMap() const {
    return range_map_;
}

---