00001
00002
00003
00004
00005 #ifndef PLAYA_SIMPLE_COMPOSED_OP_IMPL_HPP
00006 #define PLAYA_SIMPLE_COMPOSED_OP_IMPL_HPP
00007
00008
00009
00010 #include "PlayaSimpleComposedOpDecl.hpp"
00011 #include "PlayaSimpleIdentityOpDecl.hpp"
00012 #include "PlayaSimpleZeroOpDecl.hpp"
00013 #include "PlayaOut.hpp"
00014 #include "PlayaTabs.hpp"
00015
00016 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00017 #include "PlayaLinearOperatorImpl.hpp"
00018 #include "PlayaSimpleIdentityOpImpl.hpp"
00019 #include "PlayaSimpleZeroOpImpl.hpp"
00020 #endif
00021
00022
00023 namespace Playa
00024 {
00025 using namespace Teuchos;
00026
00027
00028
00029
00030
00031
00032
00033 template <class Scalar> inline
00034 SimpleComposedOp<Scalar>::SimpleComposedOp(const Array<LinearOperator<Scalar> >& ops)
00035 : LinearOpWithSpaces<Scalar>(
00036 ops[ops.size()-1].domain(), ops[0].range()
00037 )
00038 , ops_(ops)
00039 {
00040 TEUCHOS_TEST_FOR_EXCEPT(ops_.size() <= 1);
00041 for (int i=1; i<ops_.size(); i++)
00042 {
00043 TEUCHOS_TEST_FOR_EXCEPT(!(ops[i].range() == ops[i-1].domain()));
00044 }
00045 }
00046
00047
00048
00049 template <class Scalar> inline
00050 void SimpleComposedOp<Scalar>::apply(Teuchos::ETransp transApplyType,
00051 const Vector<Scalar>& in,
00052 Vector<Scalar> out) const
00053 {
00054 Tabs tab(0);
00055 PLAYA_MSG2(this->verb(), tab << "SimpleComposedOp::apply()");
00056 if (transApplyType == Teuchos::NO_TRANS)
00057 {
00058 Vector<Scalar> tmpIn = in.copy();
00059 for (int i=0; i<ops_.size(); i++)
00060 {
00061 Tabs tab1;
00062 Vector<Scalar> tmpOut;
00063 int j = ops_.size()-1-i;
00064 PLAYA_MSG2(this->verb(), tab1 << "applying op #" << j
00065 << " of " << ops_.size());
00066 ops_[j].apply(tmpIn, tmpOut);
00067 tmpIn = tmpOut;
00068 }
00069 out.acceptCopyOf(tmpIn);
00070 }
00071
00072 else if (transApplyType == Teuchos::TRANS)
00073 {
00074 Vector<Scalar> tmpIn = in.copy();
00075 for (int i=0; i<ops_.size(); i++)
00076 {
00077 Tabs tab1;
00078 Vector<Scalar> tmpOut;
00079 PLAYA_MSG2(this->verb(), tab1 << "applying transpose op #" << i
00080 << " of " << ops_.size());
00081 ops_[i].applyTranspose(tmpIn, tmpOut);
00082 tmpIn = tmpOut;
00083 }
00084 out.acceptCopyOf(tmpIn);
00085 }
00086 else
00087 {
00088 TEUCHOS_TEST_FOR_EXCEPT(transApplyType != Teuchos::TRANS && transApplyType != Teuchos::NO_TRANS);
00089 }
00090 PLAYA_MSG2(this->verb(), tab << "done SimpleComposedOp::apply()");
00091 }
00092
00093
00094
00095
00096 template <class Scalar> inline
00097 std::string SimpleComposedOp<Scalar>::description() const
00098 {
00099 std::string rtn="(";
00100 for (int i=0; i<ops_.size(); i++)
00101 {
00102 if (i > 0) rtn += "*";
00103 rtn += ops_[i].description();
00104 }
00105 rtn += ")";
00106 return rtn;
00107 }
00108
00109
00110 template <class Scalar> inline
00111 void SimpleComposedOp<Scalar>::print(std::ostream& os) const
00112 {
00113 Tabs tab(0);
00114 os << tab << "ComposedOperator[" << std::endl;
00115 for (int i=0; i<ops_.size(); i++)
00116 {
00117 Tabs tab1;
00118 os << tab1 << "factor #" << i << std::endl;
00119 Tabs tab2;
00120 os << tab2 << ops_[i].description() << std::endl;
00121 os << std::endl;
00122 }
00123 os << tab << "]" << std::endl;
00124 }
00125
00126
00127 template <class Scalar> inline
00128 LinearOperator<Scalar> composedOperator(
00129 const Array<LinearOperator<Scalar> >& ops)
00130 {
00131
00132
00133 Array<LinearOperator<Scalar> > strippedOps;
00134
00135 for (int i=0; i<ops.size(); i++)
00136 {
00137 LinearOperator<Scalar> op_i = ops[i];
00138
00139
00140
00141 const SimpleZeroOp<Scalar>* zPtr
00142 = dynamic_cast<const SimpleZeroOp<Scalar>*>(op_i.ptr().get());
00143
00144 if (zPtr != 0)
00145 {
00146 VectorSpace<Scalar> r = ops[0].range();
00147 VectorSpace<Scalar> d = ops[ops.size()-1].domain();
00148 return zeroOperator(d, r);
00149 }
00150
00151
00152 const SimpleIdentityOp<Scalar>* IPtr
00153 = dynamic_cast<const SimpleIdentityOp<Scalar>*>(op_i.ptr().get());
00154 if (IPtr != 0)
00155 {
00156 continue;
00157 }
00158
00159 strippedOps.append(op_i);
00160 }
00161
00162 TEUCHOS_TEST_FOR_EXCEPT(strippedOps.size() < 1);
00163 if (strippedOps.size()==1) return strippedOps[0];
00164
00165 RCP<LinearOperatorBase<Scalar> > op
00166 = rcp(new SimpleComposedOp<Scalar>(strippedOps));
00167 return op;
00168 }
00169
00170 template <class Scalar> inline
00171 LinearOperator<Scalar> operator*(const LinearOperator<Scalar>& A,
00172 const LinearOperator<Scalar>& B)
00173 {
00174 return composedOperator(Array<LinearOperator<Scalar> >(tuple(A,B)));
00175 }
00176
00177
00178 }
00179
00180 #endif