EpetraExt_ProductOperator.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "EpetraExt_ProductOperator.h"
00030 #include "Epetra_Vector.h"
00031 #include "Epetra_Map.h"
00032 
00033 namespace EpetraExt {
00034 
00035 // Constructors / initializers / accessors
00036 
00037 ProductOperator::ProductOperator()
00038 {}
00039 
00040 ProductOperator::ProductOperator(
00041   const int                                            num_Op
00042   ,const Teuchos::RefCountPtr<const Epetra_Operator>   Op[]
00043   ,const Teuchos::ETransp                              Op_trans[]
00044   ,const EApplyMode                                    Op_inverse[]
00045   )
00046 {
00047   initialize(num_Op,Op,Op_trans,Op_inverse);
00048 }
00049 
00050 void ProductOperator::initialize(
00051   const int                                      num_Op
00052   ,const Teuchos::RefCountPtr<const Epetra_Operator>   Op[]
00053   ,const Teuchos::ETransp                        Op_trans[]
00054   ,const EApplyMode                              Op_inverse[]
00055   )
00056 {
00057 #ifdef _DEBUG
00058   TEST_FOR_EXCEPTION(
00059     num_Op < 1, std::invalid_argument
00060     ,"ProductOperator::initialize(...): Error!"
00061     );
00062   // ToDo: Validate maps for operators!
00063 #endif // _DEBUG
00064   Op_.resize(num_Op);
00065   Op_trans_.resize(num_Op);
00066   Op_inverse_.resize(num_Op);
00067   std::copy( Op, Op + num_Op, Op_.begin() );
00068   std::copy( Op_trans, Op_trans + num_Op, Op_trans_.begin() );
00069   std::copy( Op_inverse, Op_inverse + num_Op, Op_inverse_.begin() );
00070   UseTranspose_ = false;
00071   // Wipe cache vectors so that they will be reset just to be safe
00072   range_vecs_.resize(0);
00073   domain_vecs_.resize(0);
00074 }
00075 
00076 void ProductOperator::uninitialize(
00077   int                                      *num_Op
00078   ,Teuchos::RefCountPtr<const Epetra_Operator>   Op[]
00079   ,Teuchos::ETransp                        Op_trans[]
00080   ,EApplyMode                              Op_inverse[]
00081   )
00082 {
00083 #ifdef _DEBUG
00084   TEST_FOR_EXCEPTION(
00085     (Op!=NULL || Op_trans!=NULL || Op_inverse!=NULL) && num_Op==NULL
00086     ,std::invalid_argument
00087     ,"ProductOperator::uninitialize(...): Error!"
00088     );
00089 #endif // _DEBUG
00090   if(num_Op) {
00091     *num_Op = Op_.size();
00092     if(Op) std::copy( Op_.begin(), Op_.end(), Op );
00093     if(Op_trans) std::copy( Op_trans_.begin(), Op_trans_.end(), Op_trans );
00094     if(Op_inverse) std::copy( Op_inverse_.begin(), Op_inverse_.end(), Op_inverse );
00095   }
00096   UseTranspose_ = false;
00097   Op_.resize(0);
00098   Op_trans_.resize(0);
00099   Op_inverse_.resize(0);
00100   range_vecs_.resize(0);
00101   domain_vecs_.resize(0);
00102 }
00103 
00104 void ProductOperator::applyConstituent(
00105   const int                  k
00106   ,Teuchos::ETransp          Op_trans
00107   ,EApplyMode                Op_inverse
00108   ,const Epetra_MultiVector  &X_k
00109   ,Epetra_MultiVector        *Y_k
00110   ) const
00111 {
00112   Epetra_Operator &Op_k = const_cast<Epetra_Operator&>(*Op_[k]); // Okay since we put back UseTranspose!
00113   bool oldUseTranspose = Op_k.UseTranspose();
00114   Op_k.SetUseTranspose((Op_trans==Teuchos::NO_TRANS)!=(Op_trans_[k]==Teuchos::NO_TRANS));
00115   const bool applyInverse_k = (Op_inverse==APPLY_MODE_APPLY)!=(Op_inverse_[k]==APPLY_MODE_APPLY);
00116   const int err = !applyInverse_k ? Op_[k]->Apply(X_k,*Y_k) :  Op_[k]->ApplyInverse(X_k,*Y_k);
00117   Op_k.SetUseTranspose(oldUseTranspose);
00118   TEST_FOR_EXCEPTION(
00119     err!=0, std::runtime_error
00120     ,"ProductOperator::applyConstituent(...): Error, Op["<<k<<"]."
00121     <<(!applyInverse_k?"Apply":"ApplyInverse")<<"(...) returned "
00122     "err = "<<err<<" with Op["<<k<<"].UseTranspose() = "<<Op_[k]->UseTranspose()<<"!"
00123     );
00124 }
00125 
00126 // Overridden from Epetra_Operator
00127 
00128 int ProductOperator::SetUseTranspose(bool UseTranspose)
00129 {
00130   assertInitialized();
00131   UseTranspose_ = UseTranspose;
00132   return 0;
00133 }
00134 
00135 int ProductOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00136 {
00137   assertInitialized();
00138   const int num_Op = this->num_Op();
00139   // Setup the temporary vectors
00140   initializeTempVecs(false);
00141   // Apply the constituent operators one at a time!
00142   if( !UseTranspose_ ) {
00143     //
00144     // Forward Mat-vec: Y = M * X (See main documenation)
00145     //
00146     for( int k = num_Op-1; k >= 0; --k ) {
00147       const Epetra_MultiVector  &X_k = ( k==num_Op-1 ? X : *range_vecs_[k]   );
00148       Epetra_MultiVector        &Y_k = ( k==0        ? Y : *range_vecs_[k-1] );
00149       applyConstituent(k,Teuchos::NO_TRANS,APPLY_MODE_APPLY,X_k,&Y_k);
00150     }
00151   }
00152   else if( UseTranspose_ ) {
00153     //
00154     // Adjoint Mat-vec: Y = M' * X (See main documentation)
00155     //
00156     for( int k = 0; k <= num_Op-1; ++k ) {
00157       const Epetra_MultiVector  &X_k = ( k==0         ? X : *domain_vecs_[k-1] );
00158       Epetra_MultiVector        &Y_k = ( k==num_Op-1  ? Y : *domain_vecs_[k]   );
00159       applyConstituent(k,Teuchos::TRANS,APPLY_MODE_APPLY,X_k,&Y_k);
00160     }
00161   }
00162   return 0;
00163 }
00164 
00165 int ProductOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00166 {
00167   assertInitialized();
00168   const int num_Op = this->num_Op();
00169   // Setup the temporary vectors
00170   initializeTempVecs(true);
00171   // Apply the constituent operators one at a time!
00172   if( !UseTranspose_ ) {
00173     //
00174     // Forward Inverse Mat-vec: Y = inv(M) * X (See main documenation)
00175     //
00176     for( int k = 0; k <= num_Op-1; ++k ) {
00177       const Epetra_MultiVector  &X_k = ( k==0         ? X : *domain_vecs_[k-1] );
00178       Epetra_MultiVector        &Y_k = ( k==num_Op-1  ? Y : *domain_vecs_[k]   );
00179       applyConstituent(k,Teuchos::NO_TRANS,APPLY_MODE_APPLY_INVERSE,X_k,&Y_k);
00180     }
00181   }
00182   else if( UseTranspose_ ) {
00183     //
00184     // Adjoint Invese Mat-vec: Y = inv(M') * X (See main documentation)
00185     //
00186     for( int k = num_Op-1; k >= 0; --k ) {
00187       const Epetra_MultiVector  &X_k = ( k==num_Op-1 ? X : *range_vecs_[k]   );
00188       Epetra_MultiVector        &Y_k = ( k==0        ? Y : *range_vecs_[k-1] );
00189       applyConstituent(k,Teuchos::TRANS,APPLY_MODE_APPLY_INVERSE,X_k,&Y_k);
00190     }
00191   }
00192   return 0;
00193 }
00194 
00195 double ProductOperator::NormInf() const
00196 {
00197   assertInitialized();
00198   return -1.0;
00199 }
00200 
00201 const char* ProductOperator::Label() const
00202 {
00203   assertInitialized();
00204   return NULL;
00205 }
00206 
00207 bool ProductOperator::UseTranspose() const
00208 {
00209   assertInitialized();
00210   return UseTranspose_;
00211 }
00212 
00213 bool ProductOperator::HasNormInf() const
00214 {
00215   assertInitialized();
00216   return false;
00217 }
00218 
00219 const Epetra_Comm&
00220 ProductOperator::Comm() const
00221 {
00222   assertInitialized();
00223   return Op_.front()->OperatorRangeMap().Comm();
00224 }
00225 
00226 const Epetra_Map&
00227 ProductOperator::OperatorDomainMap() const
00228 {
00229   assertInitialized();
00230   return ( Op_trans_.back()==Teuchos::NO_TRANS
00231            ? Op_.back()->OperatorDomainMap()
00232            : Op_.back()->OperatorRangeMap()
00233            );
00234 }
00235 
00236 const Epetra_Map&
00237 ProductOperator::OperatorRangeMap() const
00238 {
00239   assertInitialized();
00240   return ( Op_trans_.front()==Teuchos::NO_TRANS
00241            ? Op_.front()->OperatorRangeMap()
00242            : Op_.front()->OperatorDomainMap()
00243            );
00244 }
00245 
00246 // private
00247 
00248 void ProductOperator::initializeTempVecs(bool applyInverse) const
00249 {
00250   const int num_Op = this->num_Op();
00251   if( num_Op > 0 ) {
00252     if( ( !UseTranspose_ && !applyInverse ) || ( UseTranspose_ && applyInverse )
00253         && range_vecs_.size()==0
00254       )
00255     {
00256       //
00257       // Forward Mat-vec
00258       //
00259       // We need to create storage to hold:
00260       //
00261       //  T[k-1] = M[k]*T[k]
00262       //
00263       //    for k = num_Op-1...1
00264       //
00265       //      where: T[num_Op-1] = X (input vector)
00266       //
00267       range_vecs_.resize(num_Op-1);
00268       for( int k = num_Op-1; k >= 1; --k ) {
00269         range_vecs_[k-1] = Teuchos::rcp(
00270           new Epetra_Vector(
00271             (Op_trans_[k]==Teuchos::NO_TRANS) == (Op_inverse_[k]==APPLY_MODE_APPLY)
00272             ? Op_[k]->OperatorRangeMap()
00273             : Op_[k]->OperatorDomainMap()
00274             )
00275           );
00276       }
00277     }
00278     else if( ( UseTranspose_ && !applyInverse ) || ( !UseTranspose_ && applyInverse )
00279              && domain_vecs_.size()==0
00280       )
00281     {
00282       //
00283       // Adjoint Mat-vec
00284       //
00285       // We need to create storage to hold:
00286       //
00287       //   T[k] = M[k]'*T[k-1]
00288       //
00289       //     for k = 0...num_Op-2
00290       //
00291       //       where: T[-1]       = X (input vector)
00292       //
00293       domain_vecs_.resize(num_Op-1);
00294       for( int k = 0; k <= num_Op-2; ++k ) {
00295         domain_vecs_[k] = Teuchos::rcp(
00296           new Epetra_Vector(
00297             (Op_trans_[k]==Teuchos::NO_TRANS) == (Op_inverse_[k]==APPLY_MODE_APPLY)
00298             ? Op_[k]->OperatorDomainMap()
00299             : Op_[k]->OperatorRangeMap()
00300             )
00301           );
00302       }
00303     }
00304   }
00305 }
00306 
00307 } // namespace EpetraExt

Generated on Tue Oct 20 12:45:30 2009 for EpetraExt by doxygen 1.4.7