src/eigsolv/BlockDiagOperator.cpp

Go to the documentation of this file.

#include "BlockDiagOperator.h"

#include "Epetra_Comm.h"

//A is a symmetric matrix
//A[i,i] is nonzero

BlockDiagOperator::BlockDiagOperator(Epetra_CrsMatrix* AA)
    : A(AA),
      _useTranspose(false)
{
    _NumMyRows = A->NumMyRows();
    _diag = new int[_NumMyRows];
    _temp = new double[_NumMyRows];

    int* Indices = 0;
    double* Values = 0;
    int NumEntries;

    for (int i=0; i<_NumMyRows; i++) {
        assert(A->ExtractMyRowView(i, NumEntries, Values, Indices) >= 0);
        _diag[i] = -1;
        for (int j=0; j<NumEntries; j++) {
            if (i == Indices[j]) {
                _diag[i] = j;
                break;
            }
        }
    }
}

BlockDiagOperator::~BlockDiagOperator() {
    delete[] _diag;
    delete[] _temp;
}

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

int BlockDiagOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    int ierr = 1;
    cout << "Not implemented yet !!!!!!!!!!" << endl;
    return(ierr);
}

int BlockDiagOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    int ierr = 0;
    double* Xvalues = 0;
    double* Yvalues = 0;
    int i = 0;
    int j = 0;
    int* Indices = 0;
    double* Values = 0;
    int NumEntries = 0;
    double tsum;
    int start = 0;
    double diag = 0;

    for (int index=0; index<X.NumVectors(); index++) {

        Xvalues = X[index];
        Yvalues = Y[index];

        // (L+D) * temp = X
        for (i=0; i<_NumMyRows; i++) {
            assert(A->ExtractMyRowView(i, NumEntries, Values, Indices) >= 0);
            tsum = 0;
            for (j=0; j<_diag[i]; j++) {
                tsum += *Values++ * _temp[*Indices++];
            }
            _temp[i] = (Xvalues[i] - tsum) / *Values; 
        }

        // (L+D)^T * Y = D * temp
        for (i = _NumMyRows - 1; i >= 0; i --) {
            assert(A->ExtractMyRowView(i, NumEntries, Values, Indices) >= 0);
            tsum = 0;
            start = _diag[i];
            Values += start;
            diag = *Values++;
            Indices += ++start;
            for (j=start; j<NumEntries; j++) {
                tsum += *Values++ * Yvalues[*Indices++];
            }
            Yvalues[i] = _temp[i] - tsum / diag; 
        }
    }

    return(ierr);
}

double BlockDiagOperator::NormInf() const {
    cout << "Not implemented yet !!!!!!!!!!" << endl;
    return(1);
}

const char*  BlockDiagOperator::Label() const {
    return "Block Diagonal Operator";
}

bool BlockDiagOperator::UseTranspose() const {
    return(false);
}

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

const Epetra_Comm& BlockDiagOperator::Comm() const {
    return(A->Comm());
}

const Epetra_Map& BlockDiagOperator::OperatorDomainMap() const {
    return(A->DomainMap());
}

const Epetra_Map& BlockDiagOperator::OperatorRangeMap() const {
    return(A->RangeMap());
}

---