PlayaEpetraMatrixOps.cpp
00001
00002
00003
00004
00005 #include "PlayaEpetraMatrix.hpp"
00006 #include "PlayaEpetraMatrixFactory.hpp"
00007 #include "PlayaEpetraVector.hpp"
00008 #include "PlayaVectorSpaceDecl.hpp"
00009 #include "PlayaVectorDecl.hpp"
00010 #include "PlayaLinearOperatorDecl.hpp"
00011
00012 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00013 #include "PlayaVectorImpl.hpp"
00014 #include "PlayaLinearOperatorImpl.hpp"
00015 #endif
00016
00017 using namespace Playa;
00018 using namespace Teuchos;
00019
00020
00021
00022
00023 namespace Playa
00024 {
00025
00026 Vector<double> getEpetraDiagonal(const LinearOperator<double>& A)
00027 {
00028
00029
00030 RCP<const Epetra_CrsMatrix> A_crs = EpetraMatrix::getConcretePtr(A);
00031
00032 VectorSpace<double> rowSpace = A.domain();
00033 Vector<double> rtn = rowSpace.createMember();
00034
00035 Epetra_Vector* xPtr = EpetraVector::getConcretePtr(rtn);
00036 A_crs->ExtractDiagonalCopy(*xPtr);
00037
00038 return rtn;
00039 }
00040
00041
00042 LinearOperator<double> makeEpetraDiagonalMatrix(const Vector<double>& d)
00043 {
00044 VectorSpace<double> space = d.space();
00045 RCP<const EpetraVectorSpace> eps
00046 = rcp_dynamic_cast<const EpetraVectorSpace>(space.ptr());
00047
00048 EpetraMatrixFactory mf(eps, eps);
00049
00050 int nLocal = space.numLocalElements();
00051 int offset = space.baseGlobalNaturalIndex();
00052 for (int i=0; i<nLocal; i++)
00053 {
00054 int rowIndex = offset + i;
00055 mf.initializeNonzerosInRow(rowIndex, 1, &rowIndex);
00056 }
00057
00058 mf.finalize();
00059 LinearOperator<double> rtn = mf.createMatrix();
00060
00061 RCP<EpetraMatrix> epm = rcp_dynamic_cast<EpetraMatrix>(rtn.ptr());
00062 epm->zero();
00063
00064 for (int i=0; i<nLocal; i++)
00065 {
00066 int rowIndex = offset + i;
00067 double val = d[i];
00068 epm->addToRow(rowIndex, 1, &rowIndex, &val);
00069 }
00070
00071 return rtn;
00072 }
00073
00074 }