src/eigsolv/Projector.cpp

Go to the documentation of this file.

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

#include <pbedefs.h>

#include "Projector.h"

Projector::Projector(Epetra_Operator* Y, const Epetra_Operator* H, Epetra_Operator* C)
    : _Y(Y),
      _H(H),
      _C(C),
      _useTranspose(false),
      _vector1(new Epetra_Vector(H->OperatorRangeMap())),
      _vector2(new Epetra_Vector(H->OperatorRangeMap())),
      _vector3(new Epetra_Vector(Y->OperatorRangeMap())),
      _Comm(Y->Comm()) {

}

Projector::~Projector() {
    delete _vector1;
    delete _vector2;
    delete _vector3;
}

/*
int Projector::LinearSolverInit() {

    AMESOS::Parameter::List ParamList;
    Amesos_Factory Afactory;
    _Abase = Afactory.Create( AMESOS_MUMPS, *_LinProblem, ParamList ); 
    if (_Abase == 0) {
        cout << " AMESOS_MUMPS not implemented " << endl ; 
        exit(13);
    }  
    _LinProblem->SetOperator(_H);
    EPETRA_CHK_ERR(_Abase->SymbolicFactorization()); 
    EPETRA_CHK_ERR(_Abase->NumericFactorization()); 
} 
*/

int Projector::SetUseTranspose(bool UseTranspose) {
    int ierr = 0;
    _useTranspose = UseTranspose;
    return(ierr);
}

int Projector::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& B) const {

    int ierr = 0;
    int i = 0;
    int dim = X.NumVectors();
    const Epetra_Vector* x = 0;

    pbe_start(84, "Nullspace projector application");

    if (_useTranspose) {
        _Y->SetUseTranspose(true);
        _C->SetUseTranspose(false);
        for (i=0; i<dim; i++) {
            x = (Epetra_Vector*) X(i);
            _Y->Apply(*x, *_vector1);
            _H->ApplyInverse(*_vector1, *_vector2);         
            _C->Apply(*_vector2, *_vector3);        
            B(i)->Update(1.0, *x, -1.0, *_vector3, 0.0);
        }       
    }
    else {
        _Y->SetUseTranspose(false);
        _C->SetUseTranspose(true);      
        for (i=0; i<dim; i++) {
            x = (Epetra_Vector*) X(i);
            _C->Apply(*x, *_vector1);
            _H->ApplyInverse(*_vector1, *_vector2);
            _Y->Apply(*_vector2, *_vector3);
            B(i)->Update(1.0, *x, -1.0, *_vector3, 0.0);
        }
    }
    pbe_stop(84);

    return(ierr);
}

int  Projector::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    return(1);
}

double Projector::NormInf() const {
    return(1);
}

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

bool Projector::UseTranspose() const {
    return(_useTranspose);
}

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

const Epetra_Comm& Projector::Comm() const {
    return(_Comm);
}

// ???
const Epetra_Map& Projector::OperatorDomainMap() const {
    return(_Y->OperatorRangeMap());
}

// ???
const Epetra_Map& Projector::OperatorRangeMap() const {
    return(_Y->OperatorRangeMap());
}

---