src/eigsolv/AbstractEigsolvOperators.cpp

Go to the documentation of this file.

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

#include <cassert>
#include "myrlog.h"
#include "Epetra_Comm.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Map.h"
#include "AbstractEigsolvOperators.h"

#if 1

Epetra_CrsMatrix* AbstractEigsolvOperators::build_Asigma(const Epetra_CrsMatrix& A,
                                                         const Epetra_CrsMatrix& M,
                                                         double sigma)
{
    // FIXME: This routine could probably be sped up by inserting/adding whole rows instead of
    //        individual elements.
    const Epetra_Map& map_nedelec = A.DomainMap();
    Epetra_CrsMatrix* Asigma = new Epetra_CrsMatrix(Copy, map_nedelec, 0);
    int sizeLocal = map_nedelec.NumMyElements();
    int i;
    for (i = 0; i < sizeLocal; ++ i) {
        int numEntries = 0;
        int *indices;
        double *values;
        int i_global = map_nedelec.GID(i);
        int info;
        info = A.ExtractMyRowView(i, numEntries, values, indices);
        assert(info >= 0);
        int jj;
        for (jj = 0; jj < numEntries; ++ jj) {
            int col = A.GCID(indices[jj]);
            assert(Asigma->InsertGlobalValues(i_global, 1, values + jj, &col) >=0);
        }
        info = M.ExtractMyRowView(i, numEntries, values, indices);
        assert(info >= 0);
        for (jj = 0; jj < numEntries; ++ jj) {
            int col = M.GCID(indices[jj]);
            double val = - sigma*values[jj];
            info = Asigma->SumIntoGlobalValues(i_global, 1, &val, &col);
            assert(info >= 0);
            if (info > 0) {
                assert(Asigma->InsertGlobalValues(i_global, 1, &val, &col) >= 0);
            }
        }
    }

    Asigma->FillComplete(map_nedelec, map_nedelec);
    log_matrix_stats("Matrix K", Asigma);
    return Asigma;
}
#else

#include "EpetraExt_MatrixMatrix.h"

Epetra_CrsMatrix* AbstractEigsolvOperators::build_Asigma(const Epetra_CrsMatrix& A,
                                                         const Epetra_CrsMatrix& M,
                                                         double sigma)
{
    Epetra_CrsMatrix* Asigma = new Epetra_CrsMatrix(M);
    assert(EpetraExt::MatrixMatrix::Add(A, false, 1.0, *Asigma, -sigma) == 0);
    return Asigma;
}

#endif

---