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