00001
00002
00003
00004
00005 #include "PlayaIfpackOperator.hpp"
00006 #include "Teuchos_Array.hpp"
00007 #include "PlayaMPIComm.hpp"
00008 #include "PlayaEpetraVector.hpp"
00009
00010
00011 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00012 #include "PlayaVectorImpl.hpp"
00013 #include "PlayaLinearOperatorImpl.hpp"
00014 #endif
00015
00016 namespace Playa
00017 {
00018
00019 using namespace Teuchos;
00020
00021 IfpackOperator::IfpackOperator(const EpetraMatrix* A,
00022 int fillLevels,
00023 int overlapFill,
00024 double relaxationValue,
00025 double relativeThreshold,
00026 double absoluteThreshold)
00027 : LinearOpWithSpaces<double>(A->domain(), A->range()),
00028 precondGraph_(),
00029 precond_()
00030 {
00031 const Epetra_CrsMatrix* matrix = A->crsMatrix();
00032
00033 const Epetra_CrsGraph& matrixGraph = matrix->Graph();
00034
00035 Ifpack_IlukGraph* precondGraph
00036 = new Ifpack_IlukGraph(matrixGraph, fillLevels, overlapFill);
00037 precondGraph_ = rcp(precondGraph);
00038
00039 int ierr = precondGraph->ConstructFilledGraph();
00040
00041 TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error,
00042 "IfpackOperator ctor: "
00043 "precondGraph->ConstructFilledGraph() failed with ierr="
00044 << ierr);
00045
00046 Ifpack_CrsRiluk* precond = new Ifpack_CrsRiluk(*precondGraph);
00047 precond_ = rcp(precond);
00048
00049 precond->SetRelaxValue(relaxationValue);
00050 precond->SetRelativeThreshold(relativeThreshold);
00051 precond->SetAbsoluteThreshold(absoluteThreshold);
00052
00053 ierr = precond->InitValues(*matrix);
00054
00055 TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error,
00056 "IfpackOperator ctor: "
00057 "precond->InitValues() failed with ierr="
00058 << ierr);
00059
00060 ierr = precond->Factor();
00061
00062 TEUCHOS_TEST_FOR_EXCEPTION(ierr < 0, std::runtime_error,
00063 "IfpackOperator ctor: "
00064 "precond->Factor() failed with ierr="
00065 << ierr);
00066 }
00067
00068 void IfpackOperator::apply(Teuchos::ETransp transApplyType,
00069 const Vector<double>& in,
00070 Vector<double> out) const
00071 {
00072
00073 const Epetra_Vector& epIn = EpetraVector::getConcrete(in);
00074 Epetra_Vector& epOut = EpetraVector::getConcrete(out);
00075
00076
00077
00078 Ifpack_CrsRiluk* p = const_cast<Ifpack_CrsRiluk*>(precond_.get());
00079
00080 int ierr;
00081
00082
00083 if (transApplyType==NO_TRANS)
00084 {
00085 ierr = p->Solve(false, epIn, epOut);
00086 }
00087 else
00088 {
00089 ierr = p->Solve(true, epIn, epOut);
00090 }
00091 }
00092
00093 }