00001
00002
00003
00004
00005 #ifndef PLAYA_LINEAROPERATORIMPL_HPP
00006 #define PLAYA_LINEAROPERATORIMPL_HPP
00007
00008 #include "PlayaDefs.hpp"
00009 #include "PlayaLinearOperatorDecl.hpp"
00010 #include "Teuchos_RefCountPtr.hpp"
00011 #include "PlayaVectorDecl.hpp"
00012 #include "PlayaVectorSpaceDecl.hpp"
00013 #include "PlayaInverseOperatorDecl.hpp"
00014 #include "PlayaSimpleTransposedOpDecl.hpp"
00015 #include "PlayaBlockOperatorBaseDecl.hpp"
00016 #include "PlayaVectorType.hpp"
00017 #include "PlayaOut.hpp"
00018
00019 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00020 #include "PlayaVectorImpl.hpp"
00021 #endif
00022
00023
00024
00025 using namespace Playa;
00026 using namespace Teuchos;
00027
00028
00029 template <class Scalar>
00030 class InverseOperator;
00031
00032
00033
00034 template <class Scalar>
00035 LinearOperator<Scalar>::LinearOperator()
00036 : Handle<LinearOperatorBase<Scalar> >() {;}
00037
00038
00039
00040 template <class Scalar>
00041 LinearOperator<Scalar>::LinearOperator(const RCP<LinearOperatorBase<Scalar> >& smartPtr)
00042 : Handle<LinearOperatorBase<Scalar> >(smartPtr) {;}
00043
00044
00045
00046
00047
00048 template <class Scalar> inline
00049 void LinearOperator<Scalar>::apply(const Vector<Scalar>& in,
00050 Vector<Scalar>& out) const
00051 {
00052 Tabs tab(0);
00053 PLAYA_MSG1(this->verb(), tab << "Operator=" << this->description()
00054 << ", calling apply() function");
00055 Tabs tab1;
00056
00057 if (this->verb() > 2)
00058 {
00059 Tabs tab2;
00060 Out::os() << tab2 << "input vector = " << in << std::endl;
00061 }
00062 else if (this->verb() > 1)
00063 {
00064 Tabs tab2;
00065 Out::os() << tab2 << "input vector = " << in.description() << std::endl;
00066 }
00067
00068
00069
00070 if (out.ptr().get()==0)
00071 {
00072 Tabs tab2;
00073 PLAYA_MSG3(this->verb(), tab2 << "allocating output vector");
00074 out = this->range().createMember();
00075 }
00076 else
00077 {
00078 Tabs tab2;
00079 PLAYA_MSG3(this->verb(), tab2 << "using preallocated output vector");
00080 }
00081
00082 this->ptr()->apply(Teuchos::NO_TRANS, in, out);
00083
00084 if (this->verb() > 2)
00085 {
00086 Tabs tab2;
00087 Out::os() << tab2 << "output vector = " << out << std::endl;
00088 }
00089 else if (this->verb() > 1)
00090 {
00091 Tabs tab2;
00092 Out::os() << tab2 << "output vector = " << out.description() << std::endl;
00093 }
00094
00095 PLAYA_MSG1(this->verb(), tab << "Operator=" << this->description()
00096 << ", done with apply() function");
00097
00098 }
00099
00100
00101
00102
00103
00104 template <class Scalar> inline
00105 void LinearOperator<Scalar>::applyTranspose(const Vector<Scalar>& in,
00106 Vector<Scalar>& out) const
00107 {
00108 Tabs tab(0);
00109 PLAYA_MSG1(this->verb(), tab << "Operator=" << this->description()
00110 << ", calling applyTranspose() function");
00111 Tabs tab1;
00112
00113 if (this->verb() > 2)
00114 {
00115 Tabs tab2;
00116 Out::os() << tab2 << "input vector = " << in << std::endl;
00117 }
00118 else if (this->verb() > 1)
00119 {
00120 Tabs tab2;
00121 Out::os() << tab2 << "input vector = " << in.description() << std::endl;
00122 }
00123
00124
00125
00126
00127 if (out.ptr().get()==0)
00128 {
00129 Tabs tab2;
00130 PLAYA_MSG3(this->verb(), tab2 << "allocating output vector");
00131 out = this->domain().createMember();
00132 }
00133 else
00134 {
00135 Tabs tab2;
00136 PLAYA_MSG3(this->verb(), tab2 << "using preallocated output vector");
00137 }
00138
00139 this->ptr()->apply(Teuchos::TRANS, in, out);
00140
00141
00142
00143 if (this->verb() > 2)
00144 {
00145 Tabs tab2;
00146 Out::os() << tab2 << "output vector = " << out << std::endl;
00147 }
00148 else if (this->verb() > 1)
00149 {
00150 Tabs tab2;
00151 Out::os() << tab2 << "output vector = " << out.description() << std::endl;
00152 }
00153
00154 PLAYA_MSG1(this->verb(), tab << "Operator=" << this->description()
00155 << ", done with applyTranpose() function");
00156
00157 }
00158
00159
00160
00161 template <class Scalar>
00162 RCP<Time>& LinearOperator<Scalar>::opTimer()
00163 {
00164 static RCP<Time> rtn
00165 = TimeMonitor::getNewTimer("Low-level vector operations");
00166 return rtn;
00167 }
00168
00169
00170 template <class Scalar>
00171 LinearOperator<Scalar> LinearOperator<Scalar>::transpose() const
00172 {
00173 LinearOperator<Scalar> op = transposedOperator(*this);
00174 return op;
00175 }
00176
00177
00178
00179
00180
00181
00182 template <class Scalar>
00183 RCP<LoadableMatrix<Scalar> > LinearOperator<Scalar>::matrix()
00184 {
00185 RCP<LoadableMatrix<Scalar> > rtn
00186 = rcp_dynamic_cast<LoadableMatrix<Scalar> >(this->ptr());
00187 return rtn;
00188 }
00189
00190
00191 template <class Scalar>
00192 void LinearOperator<Scalar>::getRow(const int& row,
00193 Teuchos::Array<int>& indices,
00194 Teuchos::Array<Scalar>& values) const
00195 {
00196 const RowAccessibleOp<Scalar>* val =
00197 dynamic_cast<const RowAccessibleOp<Scalar>* >(this->ptr().get());
00198 TEUCHOS_TEST_FOR_EXCEPTION(val == 0, std::runtime_error,
00199 "Operator not row accessible; getRow() not defined.");
00200 val->getRow(row, indices, values);
00201 }
00202
00203
00204 template <class Scalar>
00205 int LinearOperator<Scalar>::numBlockRows() const
00206 {
00207 const BlockOperatorBase<Scalar>* b = dynamic_cast<const BlockOperatorBase<Scalar>* >(this->ptr().get());
00208 if (b==0) return 1;
00209 return b->numBlockRows();
00210 }
00211
00212
00213 template <class Scalar>
00214 int LinearOperator<Scalar>::numBlockCols() const
00215 {
00216 const BlockOperatorBase<Scalar>* b = dynamic_cast<const BlockOperatorBase<Scalar>* >(this->ptr().get());
00217 if (b==0) return 1;
00218 return b->numBlockCols();
00219 }
00220
00221
00222
00223 template <class Scalar>
00224 const VectorSpace<Scalar>
00225 LinearOperator<Scalar>::range() const
00226 {return this->ptr()->range();}
00227
00228
00229
00230 template <class Scalar>
00231 void LinearOperator<Scalar>::setBlock(int i, int j,
00232 const LinearOperator<Scalar>& sub)
00233 {
00234 SetableBlockOperatorBase<Scalar>* b =
00235 dynamic_cast<SetableBlockOperatorBase<Scalar>* >(this->ptr().get());
00236
00237 TEUCHOS_TEST_FOR_EXCEPTION(b == 0, std::runtime_error,
00238 "Can't call setBlock since operator not SetableBlockOperatorBase");
00239
00240 b->setBlock(i, j, sub);
00241 }
00242
00243
00244
00245
00246 template <class Scalar>
00247 const VectorSpace<Scalar>
00248 LinearOperator<Scalar>::domain() const
00249 {return this->ptr()->domain();}
00250
00251
00252
00253
00254 template <class Scalar>
00255 LinearOperator<Scalar> LinearOperator<Scalar>::getBlock(const int &i,
00256 const int &j) const
00257 {
00258 const BlockOperatorBase<Scalar>* b =
00259 dynamic_cast<const BlockOperatorBase<Scalar>* >(this->ptr().get());
00260
00261 if (b==0)
00262 {
00263 TEUCHOS_TEST_FOR_EXCEPTION(i != 0 || j != 0, std::runtime_error,
00264 "nonzero block index (" << i << "," << j << ") into "
00265 "non-block operator");
00266 return *this;
00267 }
00268 return b->getBlock(i, j);
00269 }
00270
00271
00272
00273 template <class Scalar>
00274 LinearOperator<Scalar> LinearOperator<Scalar>::getNonconstBlock(const int &i,
00275 const int &j)
00276 {
00277 BlockOperatorBase<Scalar>* b =
00278 dynamic_cast<BlockOperatorBase<Scalar>* >(this->ptr().get());
00279
00280 if (b==0)
00281 {
00282 TEUCHOS_TEST_FOR_EXCEPTION(i != 0 || j != 0, std::runtime_error,
00283 "nonzero block index (" << i << "," << j << ") into "
00284 "non-block operator");
00285 return *this;
00286 }
00287 return b->getNonconstBlock(i, j);
00288 }
00289
00290
00291
00292
00293 template <class Scalar>
00294 void LinearOperator<Scalar>::endBlockFill()
00295 {
00296 SetableBlockOperatorBase<Scalar>* b =
00297 dynamic_cast<SetableBlockOperatorBase<Scalar>* >(this->ptr().get());
00298
00299 TEUCHOS_TEST_FOR_EXCEPTION(b == 0, std::runtime_error,
00300 "Can't call endBlockFill because operator is not a SetableBlockOperator");
00301
00302
00303 b->endBlockFill();
00304 }
00305
00306
00307
00308
00309
00310
00311 #endif