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 EPETRAEXT_PRODUCT_OPERATOR_H
00030 #define EPETRAEXT_PRODUCT_OPERATOR_H
00031
00032 #include "Epetra_Operator.h"
00033 #include "Teuchos_RefCountPtr.hpp"
00034 #include "Teuchos_BLAS_types.hpp"
00035
00036 class Epetra_Vector;
00037
00038 namespace EpetraExt {
00039
00119 class ProductOperator : public Epetra_Operator {
00120 public:
00121
00124
00126 enum EApplyMode { APPLY_MODE_APPLY, APPLY_MODE_APPLY_INVERSE };
00127
00129
00132
00134 ProductOperator();
00135
00137 ProductOperator(
00138 const int num_Op
00139 ,const Teuchos::RefCountPtr<const Epetra_Operator> Op[]
00140 ,const Teuchos::ETransp Op_trans[]
00141 ,const EApplyMode Op_inverse[]
00142 );
00143
00207 void initialize(
00208 const int num_Op
00209 ,const Teuchos::RefCountPtr<const Epetra_Operator> Op[]
00210 ,const Teuchos::ETransp Op_trans[]
00211 ,const EApplyMode Op_inverse[]
00212 );
00213
00220 void uninitialize(
00221 int *num_Op
00222 ,Teuchos::RefCountPtr<const Epetra_Operator> Op[]
00223 ,Teuchos::ETransp Op_trans[]
00224 ,EApplyMode p_inverse[]
00225 );
00226
00244 void applyConstituent(
00245 const int k
00246 ,Teuchos::ETransp Op_trans
00247 ,EApplyMode Op_inverse
00248 ,const Epetra_MultiVector &X_k
00249 ,Epetra_MultiVector *Y_k
00250 ) const;
00251
00257 int num_Op() const;
00258
00270 Teuchos::RefCountPtr<const Epetra_Operator> Op(int k) const;
00271
00278 Teuchos::ETransp Op_trans(int k) const;
00279
00286 EApplyMode Op_inverse(int k) const;
00287
00289
00292
00294 int SetUseTranspose(bool UseTranspose);
00296 int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00298 int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00300 double NormInf() const;
00302 const char * Label() const;
00304 bool UseTranspose() const;
00306 bool HasNormInf() const;
00308 const Epetra_Comm & Comm() const;
00310 const Epetra_Map & OperatorDomainMap() const;
00312 const Epetra_Map & OperatorRangeMap() const;
00313
00315
00316 private:
00317
00318
00319
00320
00321 typedef std::vector<Teuchos::RefCountPtr<const Epetra_Operator> > Op_t;
00322 typedef std::vector<Teuchos::ETransp> Op_trans_t;
00323 typedef std::vector<EApplyMode> Op_inverse_t;
00324 typedef std::vector<Teuchos::RefCountPtr<Epetra_Vector> > EV_t;
00325
00326
00327
00328
00329 bool UseTranspose_;
00330 Op_t Op_;
00331 Op_trans_t Op_trans_;
00332 Op_inverse_t Op_inverse_;
00333
00334 mutable EV_t range_vecs_;
00335 mutable EV_t domain_vecs_;
00336
00337
00338
00339
00340 void assertInitialized() const;
00341 void validateIndex(int k) const;
00342 void initializeTempVecs(bool applyInverse) const;
00343
00344 };
00345
00346
00347
00348
00349
00350
00351 inline
00352 int ProductOperator::num_Op() const
00353 {
00354 return Op_.size();
00355 }
00356
00357 inline
00358 Teuchos::RefCountPtr<const Epetra_Operator>
00359 ProductOperator::Op(int k) const
00360 {
00361 validateIndex(k);
00362 return Op_[k];
00363 }
00364
00365 inline
00366 Teuchos::ETransp
00367 ProductOperator::Op_trans(int k) const
00368 {
00369 validateIndex(k);
00370 return Op_trans_[k];
00371 }
00372
00373 inline
00374 ProductOperator::EApplyMode
00375 ProductOperator::Op_inverse(int k) const
00376 {
00377 validateIndex(k);
00378 return Op_inverse_[k];
00379 }
00380
00381
00382
00383
00384 inline
00385 void ProductOperator::assertInitialized() const
00386 {
00387 TEST_FOR_EXCEPTION(
00388 Op_.size()==0, std::logic_error
00389 ,"Epetra::ProductOperator: Error, Client has not called initialize(...) yet!"
00390 );
00391 }
00392
00393 inline
00394 void ProductOperator::validateIndex(int k) const
00395 {
00396 TEST_FOR_EXCEPTION(
00397 k < 0 || static_cast<int>(Op_.size())-1 < k, std::logic_error
00398 ,"Epetra::ProductOperator: Error, k = "<<k<< " is not in the range [0,"<<Op_.size()-1<<"]!"
00399 );
00400 }
00401
00402 }
00403
00404 #endif // EPETRAEXT_PRODUCT_OPERATOR_H