00001
00002
00003
00004
00005 #ifndef PLAYA_INVERSEOPERATOR_IMPL_HPP
00006 #define PLAYA_INVERSEOPERATOR_IMPL_HPP
00007
00008 #include "PlayaDefs.hpp"
00009 #include "PlayaTabs.hpp"
00010 #include "PlayaSolverState.hpp"
00011 #include "PlayaInverseOperatorDecl.hpp"
00012 #include "PlayaSimpleZeroOpDecl.hpp"
00013
00014 #include "Teuchos_RefCountPtr.hpp"
00015
00016 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00017 #include "PlayaSimpleTransposedOpImpl.hpp"
00018 #include "PlayaSimpleZeroOpDecl.hpp"
00019 #endif
00020
00021
00022 namespace Playa
00023 {
00024 using Teuchos::RCP;
00025
00026
00027
00028
00029 template <class Scalar> inline
00030 InverseOperator<Scalar>::InverseOperator(const LinearOperator<Scalar>& op,
00031 const LinearSolver<Scalar>& solver)
00032 : LinearOpWithSpaces<Scalar>(op.domain(), op.range()),
00033 op_(op), solver_(solver) {;}
00034
00035
00036 template <class Scalar> inline
00037 void InverseOperator<Scalar>::apply(
00038 Teuchos::ETransp applyType,
00039 const Vector<Scalar>& in,
00040 Vector<Scalar> out) const
00041 {
00042
00043 Tabs tab(0);
00044 PLAYA_MSG2(this->verb(), tab << "InverseOperator::generalApply()");
00045
00046 TEUCHOS_TEST_FOR_EXCEPTION(dynamic_cast<SimpleZeroOp<Scalar>* >(op_.ptr().get()) != 0, std::runtime_error,
00047 "InverseOperator<Scalar>::apply() called on a zero operator.");
00048
00049 TEUCHOS_TEST_FOR_EXCEPTION(op_.domain().dim() != op_.range().dim(), std::runtime_error,
00050 "InverseOperator<Scalar>::apply() called on a non-square operator.");
00051
00052 Vector<Scalar> result = out.space().createMember();
00053
00054 SolverState<Scalar> haveSoln;
00055 if (applyType==NO_TRANS)
00056 {
00057 haveSoln = solver_.solve(op_, in, result);
00058 }
00059 else
00060 {
00061 haveSoln = solver_.solve(op_.transpose(), in, result);
00062 }
00063
00064
00065 TEUCHOS_TEST_FOR_EXCEPTION(haveSoln.finalState() != SolveConverged,
00066 std::runtime_error,
00067 "InverseOperator<Scalar>::apply() "
00068 << haveSoln.stateDescription());
00069
00070 out.acceptCopyOf(result);
00071
00072 PLAYA_MSG2(this->verb(), tab << "done InverseOperator::generalApply()");
00073 }
00074
00075
00076
00077 template <class Scalar>
00078 void InverseOperator<Scalar>::print(std::ostream& os) const
00079 {
00080 Tabs tab(0);
00081 os << tab << "InverseOperator[" << std::endl;
00082 Tabs tab1;
00083 os << tab1 << "op=" << op_ << std::endl;
00084 os << tab << "]" << std::endl;
00085 }
00086
00087
00088 template <class Scalar>
00089 LinearOperator<Scalar>
00090 inverse(const LinearOperator<Scalar>& op,
00091 const LinearSolver<Scalar>& solver)
00092 {
00093
00094 const InverseOperator<Scalar>* invOp
00095 = dynamic_cast<const InverseOperator<Scalar>*>(op.ptr().get());
00096 if (invOp) return invOp->op();
00097
00098
00099 RCP<LinearOperatorBase<Scalar> > rtn
00100 = rcp(new InverseOperator<Scalar>(op, solver));
00101 return rtn;
00102 }
00103
00104 }
00105
00106 #endif