PlayaSimpleAddedOpImpl.hpp
00001
00002
00003
00004
00005 #ifndef PLAYA_SIMPLE_ADDED_OP_IMPL_HPP
00006 #define PLAYA_SIMPLE_ADDED_OP_IMPL_HPP
00007
00008
00009
00010 #include "PlayaSimpleAddedOpDecl.hpp"
00011 #include "PlayaSimpleZeroOpDecl.hpp"
00012 #include "PlayaLinearCombinationImpl.hpp"
00013 #include "PlayaOut.hpp"
00014 #include "PlayaTabs.hpp"
00015 #include "Teuchos_Array.hpp"
00016
00017 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00018 #include "PlayaSimpleZeroOpImpl.hpp"
00019 #include "PlayaLinearOperatorImpl.hpp"
00020 #include "PlayaSimpleTransposedOpImpl.hpp"
00021 #endif
00022
00023
00024 namespace Playa
00025 {
00026 using namespace Teuchos;
00027
00028
00029
00030
00031
00032
00033 template <class Scalar> inline
00034 SimpleAddedOp<Scalar>::SimpleAddedOp(
00035 const Array<LinearOperator<Scalar> >& ops)
00036 : LinearOpWithSpaces<Scalar>(
00037 ops[0].domain(), ops[0].range()
00038 )
00039 , ops_(ops)
00040 {
00041 TEUCHOS_TEST_FOR_EXCEPT(ops_.size() <= 1);
00042 for (int i=1; i<ops_.size(); i++)
00043 {
00044 TEUCHOS_TEST_FOR_EXCEPT(!(ops[i].range() == ops[0].range()));
00045 TEUCHOS_TEST_FOR_EXCEPT(!(ops[i].domain() == ops[0].domain()));
00046 }
00047 }
00048
00049
00050 template <class Scalar> inline
00051 void SimpleAddedOp<Scalar>::apply(Teuchos::ETransp transApplyType,
00052 const Vector<Scalar>& in,
00053 Vector<Scalar> out) const
00054 {
00055 Tabs tab(0);
00056 PLAYA_MSG2(this->verb(), tab << "SimpleAddedOp::apply()");
00057
00058 Vector<Scalar> tmp=out.copy();
00059 tmp.zero();
00060 for (int i=0; i<ops_.size(); i++)
00061 {
00062 Tabs tab1;
00063 PLAYA_MSG3(this->verb(), tab1 << "applying term i=" << i << " of "
00064 << ops_.size());
00065 if (transApplyType == Teuchos::NO_TRANS)
00066 tmp += ops_[i] * in;
00067 else if (transApplyType == Teuchos::TRANS)
00068 tmp += ops_[i].transpose() * in;
00069 else
00070 TEUCHOS_TEST_FOR_EXCEPT(transApplyType != Teuchos::TRANS && transApplyType != Teuchos::NO_TRANS);
00071 }
00072 out.acceptCopyOf(tmp);
00073
00074 PLAYA_MSG2(this->verb(), tab << "done SimpleAddedOp::apply()");
00075 }
00076
00077
00078 template <class Scalar> inline
00079 std::string SimpleAddedOp<Scalar>::description() const
00080 {
00081 std::string rtn="(";
00082 for (int i=0; i<ops_.size(); i++)
00083 {
00084 if (i > 0) rtn += "+";
00085 rtn += ops_[i].description();
00086 }
00087 rtn += ")";
00088 return rtn;
00089 }
00090
00091
00092
00093 template <class Scalar> inline
00094 LinearOperator<Scalar> addedOperator(
00095 const Array<LinearOperator<Scalar> >& ops)
00096 {
00097
00098 Array<LinearOperator<Scalar> > strippedOps;
00099
00100 for (int i=0; i<ops.size(); i++)
00101 {
00102 LinearOperator<Scalar> op_i = ops[i];
00103
00104
00105 const SimpleZeroOp<Scalar>* zPtr
00106 = dynamic_cast<const SimpleZeroOp<Scalar>*>(op_i.ptr().get());
00107
00108 if (zPtr != 0) continue;
00109
00110 strippedOps.append(op_i);
00111 }
00112
00113 TEUCHOS_TEST_FOR_EXCEPT(strippedOps.size() < 1);
00114 if (strippedOps.size()==1) return strippedOps[0];
00115
00116 RCP<LinearOperatorBase<Scalar> > op
00117 = rcp(new SimpleAddedOp<Scalar>(strippedOps));
00118
00119 return op;
00120 }
00121
00122 template <class Scalar> inline
00123 LinearOperator<Scalar> operator+(const LinearOperator<Scalar>& A,
00124 const LinearOperator<Scalar>& B)
00125 {
00126 return addedOperator(Array<LinearOperator<Scalar> >(tuple(A, B)));
00127 }
00128
00129 }
00130
00131 #endif