00001 /* @HEADER@ */ 00002 // 00003 /* @HEADER@ */ 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 * Represent a sum of operators A_0 + A_1 + ... + A_n. 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 /* We will strip out any zero operators */ 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 /* Ignore any zero operators */ 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