00001
00002
00003
00004
00005 #include "PlayaEpetraMatrix.hpp"
00006 #include "PlayaEpetraMatrixMatrixProduct.hpp"
00007 #include "PlayaEpetraVector.hpp"
00008 #include "PlayaExceptions.hpp"
00009 #include "EpetraExt_MatrixMatrix.h"
00010
00011 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00012 #include "PlayaVectorImpl.hpp"
00013 #include "PlayaLinearOperatorImpl.hpp"
00014 #endif
00015
00016
00017
00018 namespace Playa
00019 {
00020 using namespace Teuchos;
00021
00022
00023 LinearOperator<double> epetraLeftScale(
00024 const Vector<double>& d,
00025 const LinearOperator<double>& A)
00026 {
00027
00028
00029 RCP<const Epetra_CrsMatrix> A_crs = EpetraMatrix::getConcretePtr(A);
00030
00031
00032 RCP<Epetra_CrsMatrix> mtxCopy = rcp(new Epetra_CrsMatrix(*A_crs));
00033
00034
00035
00036 const Epetra_Vector& epv = EpetraVector::getConcrete(d);
00037
00038
00039 mtxCopy->LeftScale(epv);
00040
00041 RCP<LinearOperatorBase<double> > rtn
00042 = rcp(new EpetraMatrix(mtxCopy, A.domain(), A.range()));
00043 return rtn;
00044
00045 }
00046
00047 LinearOperator<double> epetraRightScale(
00048 const LinearOperator<double>& A,
00049 const Vector<double>& d)
00050 {
00051
00052
00053 RCP<const Epetra_CrsMatrix> A_crs = EpetraMatrix::getConcretePtr(A);
00054
00055
00056 RCP<Epetra_CrsMatrix> mtxCopy = rcp(new Epetra_CrsMatrix(*A_crs));
00057
00058
00059
00060 const Epetra_Vector& epv = EpetraVector::getConcrete(d);
00061
00062
00063 mtxCopy->RightScale(epv);
00064
00065
00066 RCP<LinearOperatorBase<double> > rtn
00067 = rcp(new EpetraMatrix(mtxCopy, A.domain(), A.range()));
00068 return rtn;
00069
00070 }
00071
00072
00073 LinearOperator<double> epetraMatrixMatrixProduct(
00074 const LinearOperator<double>& A,
00075 const LinearOperator<double>& B)
00076 {
00077
00078
00079 RCP<const Epetra_CrsMatrix> A_crs = EpetraMatrix::getConcretePtr(A);
00080
00081
00082
00083 RCP<const Epetra_CrsMatrix> B_crs = EpetraMatrix::getConcretePtr(B);
00084
00085 bool transA = false;
00086 bool transB = false;
00087
00088
00089
00090 const Epetra_Map* rowmap
00091 = transA ? &(A_crs->DomainMap()) : &(A_crs->RowMap());
00092
00093
00094 RCP<Epetra_CrsMatrix> C = rcp(new Epetra_CrsMatrix(Copy, *rowmap, 1));
00095
00096
00097 int ierr
00098 = EpetraExt::MatrixMatrix::Multiply(*A_crs, transA, *B_crs, transB, *C);
00099 TEUCHOS_TEST_FOR_EXCEPTION(ierr != 0, RuntimeError,
00100 "EpetraExt Matrix-matrix multiply failed with error code ierr=" << ierr);
00101
00102
00103 RCP<LinearOperatorBase<double> > rtn
00104 = rcp(new EpetraMatrix(C, B.domain(), A.range()));
00105 return rtn;
00106
00107 }
00108
00109 }