00001 /* @HEADER@ */ 00002 // 00003 /* @HEADER@ */ 00004 00005 #include "PlayaEpetraMatrix.hpp" 00006 #include "PlayaEpetraMatrixFactory.hpp" 00007 #include "PlayaEpetraVector.hpp" 00008 #include "PlayaVectorSpaceDecl.hpp" // changed from Impl 00009 #include "PlayaVectorDecl.hpp" 00010 #include "PlayaLinearOperatorDecl.hpp" // changed from Impl 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 /* Extract the underlying Epetra matrix. Type checking is done 00029 * ny rcp_dynamic_cast, so we need no error check here. */ 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 }