EpetraExt Package Browser (Single Doxygen Collection) Development
EpetraExt_ProductOperator.cpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ***********************************************************************
00003 //
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include "EpetraExt_ProductOperator.h"
00043 #include "Epetra_Vector.h"
00044 #include "Epetra_Map.h"
00045 
00046 namespace EpetraExt {
00047 
00048 // Constructors / initializers / accessors
00049 
00050 ProductOperator::ProductOperator()
00051 {}
00052 
00053 ProductOperator::ProductOperator(
00054   const int                                            num_Op
00055   ,const Teuchos::RefCountPtr<const Epetra_Operator>   Op[]
00056   ,const Teuchos::ETransp                              Op_trans[]
00057   ,const EApplyMode                                    Op_inverse[]
00058   )
00059 {
00060   initialize(num_Op,Op,Op_trans,Op_inverse);
00061 }
00062 
00063 void ProductOperator::initialize(
00064   const int                                      num_Op
00065   ,const Teuchos::RefCountPtr<const Epetra_Operator>   Op[]
00066   ,const Teuchos::ETransp                        Op_trans[]
00067   ,const EApplyMode                              Op_inverse[]
00068   )
00069 {
00070 #ifdef _DEBUG
00071   TEUCHOS_TEST_FOR_EXCEPTION(
00072     num_Op < 1, std::invalid_argument
00073     ,"ProductOperator::initialize(...): Error!"
00074     );
00075   // ToDo: Validate maps for operators!
00076 #endif // _DEBUG
00077   Op_.resize(num_Op);
00078   Op_trans_.resize(num_Op);
00079   Op_inverse_.resize(num_Op);
00080   std::copy( Op, Op + num_Op, Op_.begin() );
00081   std::copy( Op_trans, Op_trans + num_Op, Op_trans_.begin() );
00082   std::copy( Op_inverse, Op_inverse + num_Op, Op_inverse_.begin() );
00083   UseTranspose_ = false;
00084   // Wipe cache vectors so that they will be reset just to be safe
00085   range_vecs_.resize(0);
00086   domain_vecs_.resize(0);
00087 }
00088 
00089 void ProductOperator::uninitialize(
00090   int                                      *num_Op
00091   ,Teuchos::RefCountPtr<const Epetra_Operator>   Op[]
00092   ,Teuchos::ETransp                        Op_trans[]
00093   ,EApplyMode                              Op_inverse[]
00094   )
00095 {
00096 #ifdef _DEBUG
00097   TEUCHOS_TEST_FOR_EXCEPTION(
00098     (Op!=NULL || Op_trans!=NULL || Op_inverse!=NULL) && num_Op==NULL
00099     ,std::invalid_argument
00100     ,"ProductOperator::uninitialize(...): Error!"
00101     );
00102 #endif // _DEBUG
00103   if(num_Op) {
00104     *num_Op = Op_.size();
00105     if(Op) std::copy( Op_.begin(), Op_.end(), Op );
00106     if(Op_trans) std::copy( Op_trans_.begin(), Op_trans_.end(), Op_trans );
00107     if(Op_inverse) std::copy( Op_inverse_.begin(), Op_inverse_.end(), Op_inverse );
00108   }
00109   UseTranspose_ = false;
00110   Op_.resize(0);
00111   Op_trans_.resize(0);
00112   Op_inverse_.resize(0);
00113   range_vecs_.resize(0);
00114   domain_vecs_.resize(0);
00115 }
00116 
00117 void ProductOperator::applyConstituent(
00118   const int                  k
00119   ,Teuchos::ETransp          Op_trans
00120   ,EApplyMode                Op_inverse
00121   ,const Epetra_MultiVector  &X_k
00122   ,Epetra_MultiVector        *Y_k
00123   ) const
00124 {
00125   Epetra_Operator &Op_k = const_cast<Epetra_Operator&>(*Op_[k]); // Okay since we put back UseTranspose!
00126   bool oldUseTranspose = Op_k.UseTranspose();
00127   Op_k.SetUseTranspose((Op_trans==Teuchos::NO_TRANS)!=(Op_trans_[k]==Teuchos::NO_TRANS));
00128   const bool applyInverse_k = (Op_inverse==APPLY_MODE_APPLY)!=(Op_inverse_[k]==APPLY_MODE_APPLY);
00129   const int err = !applyInverse_k ? Op_[k]->Apply(X_k,*Y_k) :  Op_[k]->ApplyInverse(X_k,*Y_k);
00130   Op_k.SetUseTranspose(oldUseTranspose);
00131   TEUCHOS_TEST_FOR_EXCEPTION(
00132     err!=0, std::runtime_error
00133     ,"ProductOperator::applyConstituent(...): Error, Op["<<k<<"]."
00134     <<(!applyInverse_k?"Apply":"ApplyInverse")<<"(...) returned "
00135     "err = "<<err<<" with Op["<<k<<"].UseTranspose() = "<<Op_[k]->UseTranspose()<<"!"
00136     );
00137 }
00138 
00139 // Overridden from Epetra_Operator
00140 
00141 int ProductOperator::SetUseTranspose(bool UseTranspose)
00142 {
00143   assertInitialized();
00144   UseTranspose_ = UseTranspose;
00145   return 0;
00146 }
00147 
00148 int ProductOperator::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00149 {
00150   assertInitialized();
00151   const int num_Op = this->num_Op();
00152   // Setup the temporary vectors
00153   initializeTempVecs(false);
00154   // Apply the constituent operators one at a time!
00155   if( !UseTranspose_ ) {
00156     //
00157     // Forward Mat-vec: Y = M * X (See main documenation)
00158     //
00159     for( int k = num_Op-1; k >= 0; --k ) {
00160       const Epetra_MultiVector  &X_k = ( k==num_Op-1 ? X : *range_vecs_[k]   );
00161       Epetra_MultiVector        &Y_k = ( k==0        ? Y : *range_vecs_[k-1] );
00162       applyConstituent(k,Teuchos::NO_TRANS,APPLY_MODE_APPLY,X_k,&Y_k);
00163     }
00164   }
00165   else if( UseTranspose_ ) {
00166     //
00167     // Adjoint Mat-vec: Y = M' * X (See main documentation)
00168     //
00169     for( int k = 0; k <= num_Op-1; ++k ) {
00170       const Epetra_MultiVector  &X_k = ( k==0         ? X : *domain_vecs_[k-1] );
00171       Epetra_MultiVector        &Y_k = ( k==num_Op-1  ? Y : *domain_vecs_[k]   );
00172       applyConstituent(k,Teuchos::TRANS,APPLY_MODE_APPLY,X_k,&Y_k);
00173     }
00174   }
00175   return 0;
00176 }
00177 
00178 int ProductOperator::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
00179 {
00180   assertInitialized();
00181   const int num_Op = this->num_Op();
00182   // Setup the temporary vectors
00183   initializeTempVecs(true);
00184   // Apply the constituent operators one at a time!
00185   if( !UseTranspose_ ) {
00186     //
00187     // Forward Inverse Mat-vec: Y = inv(M) * X (See main documenation)
00188     //
00189     for( int k = 0; k <= num_Op-1; ++k ) {
00190       const Epetra_MultiVector  &X_k = ( k==0         ? X : *domain_vecs_[k-1] );
00191       Epetra_MultiVector        &Y_k = ( k==num_Op-1  ? Y : *domain_vecs_[k]   );
00192       applyConstituent(k,Teuchos::NO_TRANS,APPLY_MODE_APPLY_INVERSE,X_k,&Y_k);
00193     }
00194   }
00195   else if( UseTranspose_ ) {
00196     //
00197     // Adjoint Invese Mat-vec: Y = inv(M') * X (See main documentation)
00198     //
00199     for( int k = num_Op-1; k >= 0; --k ) {
00200       const Epetra_MultiVector  &X_k = ( k==num_Op-1 ? X : *range_vecs_[k]   );
00201       Epetra_MultiVector        &Y_k = ( k==0        ? Y : *range_vecs_[k-1] );
00202       applyConstituent(k,Teuchos::TRANS,APPLY_MODE_APPLY_INVERSE,X_k,&Y_k);
00203     }
00204   }
00205   return 0;
00206 }
00207 
00208 double ProductOperator::NormInf() const
00209 {
00210   assertInitialized();
00211   return -1.0;
00212 }
00213 
00214 const char* ProductOperator::Label() const
00215 {
00216   assertInitialized();
00217   return NULL;
00218 }
00219 
00220 bool ProductOperator::UseTranspose() const
00221 {
00222   assertInitialized();
00223   return UseTranspose_;
00224 }
00225 
00226 bool ProductOperator::HasNormInf() const
00227 {
00228   assertInitialized();
00229   return false;
00230 }
00231 
00232 const Epetra_Comm&
00233 ProductOperator::Comm() const
00234 {
00235   assertInitialized();
00236   return Op_.front()->OperatorRangeMap().Comm();
00237 }
00238 
00239 const Epetra_Map&
00240 ProductOperator::OperatorDomainMap() const
00241 {
00242   assertInitialized();
00243   return ( Op_trans_.back()==Teuchos::NO_TRANS
00244            ? Op_.back()->OperatorDomainMap()
00245            : Op_.back()->OperatorRangeMap()
00246            );
00247 }
00248 
00249 const Epetra_Map&
00250 ProductOperator::OperatorRangeMap() const
00251 {
00252   assertInitialized();
00253   return ( Op_trans_.front()==Teuchos::NO_TRANS
00254            ? Op_.front()->OperatorRangeMap()
00255            : Op_.front()->OperatorDomainMap()
00256            );
00257 }
00258 
00259 // private
00260 
00261 void ProductOperator::initializeTempVecs(bool applyInverse) const
00262 {
00263   const int num_Op = this->num_Op();
00264   if( num_Op > 0 ) {
00265     if( ( !UseTranspose_ && !applyInverse ) || ( UseTranspose_ && applyInverse )
00266         && range_vecs_.size()==0
00267       )
00268     {
00269       //
00270       // Forward Mat-vec
00271       //
00272       // We need to create storage to hold:
00273       //
00274       //  T[k-1] = M[k]*T[k]
00275       //
00276       //    for k = num_Op-1...1
00277       //
00278       //      where: T[num_Op-1] = X (input vector)
00279       //
00280       range_vecs_.resize(num_Op-1);
00281       for( int k = num_Op-1; k >= 1; --k ) {
00282         range_vecs_[k-1] = Teuchos::rcp(
00283           new Epetra_Vector(
00284             (Op_trans_[k]==Teuchos::NO_TRANS) == (Op_inverse_[k]==APPLY_MODE_APPLY)
00285             ? Op_[k]->OperatorRangeMap()
00286             : Op_[k]->OperatorDomainMap()
00287             )
00288           );
00289       }
00290     }
00291     else if( ( UseTranspose_ && !applyInverse ) || ( !UseTranspose_ && applyInverse )
00292              && domain_vecs_.size()==0
00293       )
00294     {
00295       //
00296       // Adjoint Mat-vec
00297       //
00298       // We need to create storage to hold:
00299       //
00300       //   T[k] = M[k]'*T[k-1]
00301       //
00302       //     for k = 0...num_Op-2
00303       //
00304       //       where: T[-1]       = X (input vector)
00305       //
00306       domain_vecs_.resize(num_Op-1);
00307       for( int k = 0; k <= num_Op-2; ++k ) {
00308         domain_vecs_[k] = Teuchos::rcp(
00309           new Epetra_Vector(
00310             (Op_trans_[k]==Teuchos::NO_TRANS) == (Op_inverse_[k]==APPLY_MODE_APPLY)
00311             ? Op_[k]->OperatorDomainMap()
00312             : Op_[k]->OperatorRangeMap()
00313             )
00314           );
00315       }
00316     }
00317   }
00318 }
00319 
00320 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines