src/diagonalpreconditioner.h

Go to the documentation of this file.

//
// C++ Interface: diagonalpreconditioner
//
// Description: Diagonal preconditioner
//
// Author: Roman Geus <roman.geus@psi.ch>, (C) 2004
//
// Copyright: See COPYING file that comes with this distribution
//
//

#ifndef DIAGONAL_PRECONDITIONER_H
#define DIAGONAL_PRECONDITIONER_H

#include "Epetra_Operator.h"
#include "Epetra_Vector.h"
#include "Epetra_CrsMatrix.h"

class Epetra_MultiVector;
class Epetra_Map;
class Epetra_Comm;

class DiagonalPreconditioner : public Epetra_Operator {
public:
    DiagonalPreconditioner(const Epetra_CrsMatrix& mat)
        : inv_diag_(new Epetra_Vector(mat.RowMap(), false)),
          domain_map_(mat.OperatorDomainMap()),
          range_map_(mat.OperatorRangeMap()),
          comm_(mat.Comm()),
          use_transpose_(false)
    {     
        mat.ExtractDiagonalCopy(*inv_diag_);
        for (int i = 0; i < inv_diag_->Map().NumMyElements(); i ++)
            (*inv_diag_)[i] = 1.0 / (*inv_diag_)[i];
    }

    ~DiagonalPreconditioner() {
        delete inv_diag_;
    }

    int SetUseTranspose(bool use_transpose) {
        use_transpose_ = use_transpose;
        return 0;
    }

    int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
        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 (!domain_map_.SameAs(range_map_))
            EPETRA_CHK_ERR(-4);

        return Y.Multiply(1.0, *inv_diag_, X, 0.0);
    }

    int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
        return Apply(X, Y);
    }

    double NormInf() const {
        return 0.0;
    }

    const char* Label() const {
        return "Identity operator";
    }

    bool UseTranspose() const {
        return use_transpose_;
    }

    bool HasNormInf() const {
        return false;
    }

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

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

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

private:
    Epetra_Vector* inv_diag_;
    const Epetra_Map& domain_map_;
    const Epetra_Map& range_map_;
    const Epetra_Comm& comm_;
    bool use_transpose_;
};

#endif

---