src/operatortest.cpp

Go to the documentation of this file.

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

#include "Epetra_MultiVector.h"
#include "Epetra_Vector.h"
#include "Epetra_Operator.h"
#include "Epetra_Map.h"
#include "Epetra_Comm.h"

#include "operatortest.h"

#define EPETRA_CHK(EXPR) { \
    int ierr = EXPR; \
    if (op_->Comm().MyPID() == 0 && ierr != 0) { \
        std::cerr << "EPETRA error caught in " << __FILE__ << ", line "<< __LINE__ << " : " \
                  << #EXPR << " == " << ierr << std::endl; \
    } \
}

using std::endl;

const double OperatorTest::symmetric_tol_ = 1e-13;
const double OperatorTest::linear_tol_ = 1e-11;

OperatorTest::OperatorTest(const Epetra_Operator* op)
    : op_(op),
      use_inverse_(false)
{
    verbose_ = 1;
    nrows_ = op->OperatorRangeMap().NumGlobalElements();
    ncols_ = op->OperatorDomainMap().NumGlobalElements();
}

void OperatorTest::set_operator(const Epetra_Operator* op, bool use_inverse) {
    op_ = op;
    use_inverse_ = use_inverse;
    nrows_ = op->OperatorRangeMap().NumGlobalElements();
    ncols_ = op->OperatorDomainMap().NumGlobalElements();
}

void OperatorTest::dump_stats(std::ostream& ostr, std::string operator_name) {
    bool symmetric = is_symmetric();
    bool linear = is_linear();
    bool non_destructive = is_non_destructive();

    if (op_->Comm().MyPID() == 0) {
        ostr << "Statistics for operator \"" << operator_name << "\"" << endl;
        if (symmetric) {
            ostr << "    Appears to be symmetric" << endl;
        } else {
            ostr << "    Is non-symmetric" << endl;
        }
        if (linear) {
            ostr << "    Appears to be linear" << endl;
        } else {
            ostr << "    Is non-linear" << endl;
        }
        if (non_destructive) {
            ostr << "    Appears to be behaving nicely" << endl;
        } else {
            ostr << "    Overwrites the input vector!!!" << endl;
        }
        ostr << endl;
    }
}


bool OperatorTest::is_symmetric() {
    if (ncols_ != nrows_) {
        return false;
    } else {
        Epetra_Vector x(op_->OperatorDomainMap());
        Epetra_Vector y(op_->OperatorDomainMap());
        Epetra_Vector t(op_->OperatorRangeMap());

        x.Random();
        y.Random();

        // s1 = x^T*OP*y
        double s1;
        EPETRA_CHK(apply(y, t));
        EPETRA_CHK(t.Dot(x, &s1));

        // s2 = y^T*OP*x
        double s2;
        EPETRA_CHK(apply(x, t));
        EPETRA_CHK(t.Dot(y, &s2));

        double err = fabs(s1 - s2)/fabs(s1);

        if (verbose_ > 0 && op_->Comm().MyPID() == 0)
            std::cout << "OperatorTest::is_symmetric: error=" << err << endl;
        return fabs(s1 - s2)/s1 < symmetric_tol_;
    }
}

bool OperatorTest::is_linear() {
    Epetra_Vector x(op_->OperatorDomainMap());
    Epetra_Vector y(op_->OperatorDomainMap());
    Epetra_Vector t1(op_->OperatorDomainMap());
    Epetra_Vector z1(op_->OperatorRangeMap());
    Epetra_Vector z2(op_->OperatorRangeMap());
    Epetra_Vector t2(op_->OperatorRangeMap());

    x.Random();
    y.Random();

    // z1 = OP * (x + y)
    t1 = x;
    EPETRA_CHK(t1.Update(1.0, y, 1.0));
    EPETRA_CHK(apply(t1, z1));

    // z2 = OP*x + OP*y
    EPETRA_CHK(apply(x, z2));
    EPETRA_CHK(apply(y, t2));
    EPETRA_CHK(z2.Update(1.0, t2, 1.0));

    // test z1 == z2
    EPETRA_CHK(z1.Update(-1.0, z2, 1.0));
    double z1_nrm2;
    EPETRA_CHK(z1.Norm2(&z1_nrm2));

    if (verbose_ > 0 && op_->Comm().MyPID() == 0)
        std::cout << "OperatorTest::is_linear: error=" << z1_nrm2 << endl;
    return z1_nrm2 < linear_tol_;
}

bool OperatorTest::is_non_destructive() {
    Epetra_Vector x(op_->OperatorDomainMap());
    Epetra_Vector y(op_->OperatorRangeMap());
    Epetra_Vector t(x);

    EPETRA_CHK(apply(x, y));

    // test if x has not changed
    EPETRA_CHK(t.Update(-1.0, x, 1.0));
    double t_nrm2;
    EPETRA_CHK(t.Norm2(&t_nrm2));
    if (verbose_ > 0 && op_->Comm().MyPID() == 0)
        std::cout << "OperatorTest::is_non_destructive: error=" << t_nrm2 << endl;
    return t_nrm2 == 0.0;
}


int OperatorTest::apply(Epetra_MultiVector& X, Epetra_MultiVector& Y) {
    if (use_inverse_)
        return op_->ApplyInverse(X, Y);
    else
        return op_->Apply(X, Y);
}

---