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     // FIXME (mfh 24 Mar 2014): I added the parentheses around the ||
00266     // below to silence a compiler warning.  I'm concerned that the
00267     // original author of that code didn't understand that && takes
00268     // precedence over ||, but I didn't want to change the meaning of
00269     // the original code.
00270     if (((! UseTranspose_ && ! applyInverse) || (UseTranspose_ && applyInverse))
00271         && range_vecs_.size () == 0) {
00272       //
00273       // Forward Mat-vec
00274       //
00275       // We need to create storage to hold:
00276       //
00277       //  T[k-1] = M[k]*T[k]
00278       //
00279       //    for k = num_Op-1...1
00280       //
00281       //      where: T[num_Op-1] = X (input vector)
00282       //
00283       range_vecs_.resize (num_Op - 1);
00284       for (int k = num_Op-1; k >= 1; --k) {
00285         range_vecs_[k-1] = Teuchos::rcp (new Epetra_Vector ((Op_trans_[k]==Teuchos::NO_TRANS) == (Op_inverse_[k]==APPLY_MODE_APPLY)
00286                                                             ? Op_[k]->OperatorRangeMap ()
00287                                                             : Op_[k]->OperatorDomainMap ()));
00288       }
00289     }
00290     // FIXME (mfh 24 Mar 2014): I added the parentheses around the ||
00291     // below to silence a compiler warning.  I'm concerned that the
00292     // original author of that code didn't understand that && takes
00293     // precedence over ||, but I didn't want to change the meaning of
00294     // the original code.
00295     else if (((UseTranspose_ && ! applyInverse) || (! UseTranspose_ && applyInverse))
00296              && domain_vecs_.size () == 0) {
00297       //
00298       // Adjoint Mat-vec
00299       //
00300       // We need to create storage to hold:
00301       //
00302       //   T[k] = M[k]'*T[k-1]
00303       //
00304       //     for k = 0...num_Op-2
00305       //
00306       //       where: T[-1]       = X (input vector)
00307       //
00308       domain_vecs_.resize (num_Op - 1);
00309       for (int k = 0; k <= num_Op - 2; ++k) {
00310         domain_vecs_[k] = Teuchos::rcp (new Epetra_Vector ((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 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines