00001
00002
00003
00004
00005 #include "PlayaEpetraMatrix.hpp"
00006 #include "PlayaEpetraVector.hpp"
00007 #include "PlayaVectorSpaceDecl.hpp"
00008 #include "PlayaVectorDecl.hpp"
00009 #include "PlayaLinearOperatorDecl.hpp"
00010 #include "Teuchos_Array.hpp"
00011 #include "PlayaMPIComm.hpp"
00012 #include "PlayaIfpackOperator.hpp"
00013 #include "PlayaGenericLeftPreconditioner.hpp"
00014 #include "PlayaGenericRightPreconditioner.hpp"
00015 #include "Teuchos_dyn_cast.hpp"
00016 #include "Teuchos_getConst.hpp"
00017 #include "EpetraPlayaOperator.hpp"
00018 #include "Epetra_Comm.h"
00019 #include "Epetra_CrsMatrix.h"
00020
00021
00022 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00023 #include "PlayaVectorImpl.hpp"
00024 #include "PlayaLinearCombinationImpl.hpp"
00025 #include "PlayaLinearOperatorImpl.hpp"
00026 #include "PlayaLinearSolverImpl.hpp"
00027 #endif
00028
00029 using namespace Playa;
00030 using namespace Teuchos;
00031
00032 namespace Epetra
00033 {
00034
00035 Epetra_PlayaOperator::Epetra_PlayaOperator(const LinearOperator<double>& A,
00036 const LinearSolver<double>& solver)
00037 : A_(A), solver_(solver), useTranspose_(false), comm_(), domain_(), range_(),
00038 isNativeEpetra_(false), isCompoundEpetra_(false), label_(A.description())
00039 {
00040 const EpetraMatrix* em = dynamic_cast<const EpetraMatrix*>(A.ptr().get());
00041 const EpetraVectorSpace* ed = dynamic_cast<const EpetraVectorSpace*>(A.domain().ptr().get());
00042 const EpetraVectorSpace* er = dynamic_cast<const EpetraVectorSpace*>(A.range().ptr().get());
00043
00044 if (em)
00045 {
00046 isNativeEpetra_ = true;
00047 const Epetra_CrsMatrix* crs = em->crsMatrix();
00048 domain_ = rcp(new Epetra_Map(crs->OperatorDomainMap()));
00049 range_ = rcp(new Epetra_Map(crs->OperatorRangeMap()));
00050 useTranspose_ = crs->UseTranspose();
00051 comm_ = rcp(crs->Comm().Clone());
00052 }
00053 else if (er != 0 && ed != 0)
00054 {
00055 domain_ = ed->epetraMap();
00056 range_ = er->epetraMap();
00057 comm_ = rcp(domain_->Comm().Clone());
00058 isCompoundEpetra_ = true;
00059 }
00060 else
00061 {
00062 TEUCHOS_TEST_FOR_EXCEPT(true);
00063 }
00064 }
00065
00066
00067 int Epetra_PlayaOperator::Apply(const Epetra_MultiVector& in, Epetra_MultiVector& out) const
00068 {
00069 if (isNativeEpetra_)
00070 {
00071 const EpetraMatrix* em = dynamic_cast<const EpetraMatrix*>(A_.ptr().get());
00072 return em->crsMatrix()->Multiply(useTranspose_, in, out);
00073 }
00074 else if (isCompoundEpetra_)
00075 {
00076 const Epetra_Vector* cevIn = dynamic_cast<const Epetra_Vector*>(&in);
00077 Epetra_Vector* evIn = const_cast<Epetra_Vector*>(cevIn);
00078 Epetra_Vector* evOut = dynamic_cast<Epetra_Vector*>(&out);
00079 TEUCHOS_TEST_FOR_EXCEPTION(evIn==0, std::runtime_error, "Epetra_PlayaOperator::Apply "
00080 "cannot deal with multivectors");
00081 TEUCHOS_TEST_FOR_EXCEPTION(evOut==0, std::runtime_error, "Epetra_PlayaOperator::Apply "
00082 "cannot deal with multivectors");
00083
00084
00085 RCP<VectorBase<double> > vpIn
00086 = rcp(new EpetraVector(A_.domain(), rcp(evIn, false)));
00087 RCP<VectorBase<double> > vpOut
00088 = rcp(new EpetraVector(A_.range(), rcp(evOut, false)));
00089 Vector<double> vIn = vpIn;
00090 Vector<double> vOut = vpOut;
00091
00092 A_.apply(vIn, vOut);
00093 out = EpetraVector::getConcrete(vOut);
00094 return 0;
00095 }
00096 else
00097 {
00098 TEUCHOS_TEST_FOR_EXCEPT(true);
00099 return -1;
00100 }
00101 }
00102
00103 int Epetra_PlayaOperator::ApplyInverse(const Epetra_MultiVector& in, Epetra_MultiVector& out) const
00104 {
00105
00106 TEUCHOS_TEST_FOR_EXCEPTION(solver_.ptr().get()==0, std::runtime_error,
00107 "no solver provided for Epetra_PlayaOperator::ApplyInverse");
00108 TEUCHOS_TEST_FOR_EXCEPTION(!isNativeEpetra_ && !isCompoundEpetra_, std::runtime_error,
00109 "Epetra_PlayaOperator::ApplyInverse expects either "
00110 "a native epetra operator or a compound operator with "
00111 "Epetra domain and range spaces");
00112 const Epetra_Vector* cevIn = dynamic_cast<const Epetra_Vector*>(&in);
00113 Epetra_Vector* evIn = const_cast<Epetra_Vector*>(cevIn);
00114 Epetra_Vector* evOut = dynamic_cast<Epetra_Vector*>(&out);
00115
00116 TEUCHOS_TEST_FOR_EXCEPTION(evIn==0, std::runtime_error, "Epetra_PlayaOperator::Apply "
00117 "cannot deal with multivectors");
00118 TEUCHOS_TEST_FOR_EXCEPTION(evOut==0, std::runtime_error, "Epetra_PlayaOperator::Apply "
00119 "cannot deal with multivectors");
00120
00121 RCP<VectorBase<double> > vpIn
00122 = rcp(new EpetraVector(A_.domain(),
00123 rcp(evIn, false)));
00124 RCP<VectorBase<double> > vpOut
00125 = rcp(new EpetraVector(A_.range(),
00126 rcp(evOut, false)));
00127 Vector<double> vIn = vpIn;
00128 Vector<double> vOut = vpOut;
00129
00130 SolverState<double> state = solver_.solve(A_, vIn, vOut);
00131
00132 if (state.finalState() == SolveCrashed) return -1;
00133 else if (state.finalState() == SolveFailedToConverge) return -2;
00134 else out = EpetraVector::getConcrete(vOut);
00135
00136 return 0;
00137 }
00138
00139
00140
00141
00142 double Epetra_PlayaOperator::NormInf() const
00143 {
00144 TEUCHOS_TEST_FOR_EXCEPT(true);
00145 return -1;
00146 }
00147
00148 const char* Epetra_PlayaOperator::Label() const
00149 {
00150 return label_.c_str();
00151 }
00152
00153 }