Teko Version of the Day
Teko_EpetraBlockPreconditioner.cpp
00001 /*
00002 // @HEADER
00003 // 
00004 // ***********************************************************************
00005 // 
00006 //      Teko: A package for block and physics based preconditioning
00007 //                  Copyright 2010 Sandia Corporation 
00008 //  
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 //  
00012 // Redistribution and use in source and binary forms, with or without
00013 // modification, are permitted provided that the following conditions are
00014 // met:
00015 //  
00016 // 1. Redistributions of source code must retain the above copyright
00017 // notice, this list of conditions and the following disclaimer.
00018 //  
00019 // 2. Redistributions in binary form must reproduce the above copyright
00020 // notice, this list of conditions and the following disclaimer in the
00021 // documentation and/or other materials provided with the distribution.
00022 //  
00023 // 3. Neither the name of the Corporation nor the names of the
00024 // contributors may be used to endorse or promote products derived from
00025 // this software without specific prior written permission. 
00026 //  
00027 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00028 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00029 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00030 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00031 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00032 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00033 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00034 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00035 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
00036 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00037 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00038 //  
00039 // Questions? Contact Eric C. Cyr (eccyr@sandia.gov)
00040 // 
00041 // ***********************************************************************
00042 // 
00043 // @HEADER
00044 
00045 */
00046 
00047 #include "Teko_EpetraBlockPreconditioner.hpp"
00048 #include "Teko_Preconditioner.hpp"
00049 
00050 // Thyra includes
00051 #include "Thyra_DefaultLinearOpSource.hpp"
00052 #include "Thyra_EpetraLinearOp.hpp"
00053 
00054 // Teuchos includes
00055 #include "Teuchos_Time.hpp"
00056 
00057 // Teko includes
00058 #include "Teko_BasicMappingStrategy.hpp"
00059 
00060 namespace Teko {
00061 namespace Epetra {
00062 
00063 using Teuchos::RCP;
00064 using Teuchos::rcp;
00065 using Teuchos::rcpFromRef;
00066 using Teuchos::rcp_dynamic_cast;
00067 
00074 EpetraBlockPreconditioner::EpetraBlockPreconditioner(const Teuchos::RCP<const PreconditionerFactory> & bfp)
00075    : preconFactory_(bfp), firstBuildComplete_(false)
00076 { 
00077 }
00078 
00079 void EpetraBlockPreconditioner::initPreconditioner(bool clearOld)
00080 {
00081    if((not clearOld) && preconObj_!=Teuchos::null)
00082       return;
00083    preconObj_ = preconFactory_->createPrec();
00084 }
00085 
00098 void EpetraBlockPreconditioner::buildPreconditioner(const Teuchos::RCP<const Epetra_Operator> & A,bool clear)
00099 {
00100    Teko_DEBUG_SCOPE("EBP::buildPreconditioner",10);
00101 
00102    // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
00103    RCP<const Thyra::LinearOpBase<double> > thyraA = extractLinearOp(A);
00104 
00105    // set the mapping strategy
00106    // SetMapStrategy(rcp(new InverseMappingStrategy(eow->getMapStrategy())));
00107    SetMapStrategy(rcp(new InverseMappingStrategy(extractMappingStrategy(A))));
00108 
00109    // build preconObj_ 
00110    initPreconditioner(clear);
00111    
00112    // actually build the preconditioner
00113    RCP<const Thyra::LinearOpSourceBase<double> > lOpSrc = Thyra::defaultLinearOpSource(thyraA);
00114    preconFactory_->initializePrec(lOpSrc,&*preconObj_,Thyra::SUPPORT_SOLVE_UNSPECIFIED);
00115 
00116    // extract preconditioner operator
00117    RCP<const Thyra::LinearOpBase<double> > preconditioner = preconObj_->getUnspecifiedPrecOp();
00118 
00119    SetOperator(preconditioner,false);
00120 
00121    firstBuildComplete_ = true;
00122 
00123    TEUCHOS_ASSERT(preconObj_!=Teuchos::null);
00124    TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
00125    TEUCHOS_ASSERT(getMapStrategy()!=Teuchos::null);
00126    TEUCHOS_ASSERT(firstBuildComplete_==true);
00127 }
00128 
00142 void EpetraBlockPreconditioner::buildPreconditioner(const Teuchos::RCP<const Epetra_Operator> & A,const Epetra_MultiVector & epetra_mv,bool clear)
00143 {
00144    Teko_DEBUG_SCOPE("EBP::buildPreconditioner - with solution",10);
00145 
00146    // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
00147    RCP<const Thyra::LinearOpBase<double> > thyraA = extractLinearOp(A);
00148 
00149    // set the mapping strategy
00150    SetMapStrategy(rcp(new InverseMappingStrategy(extractMappingStrategy(A))));
00151 
00152    TEUCHOS_ASSERT(getMapStrategy()!=Teuchos::null);
00153    
00154    // build the thyra version of the source multivector
00155    RCP<Thyra::MultiVectorBase<double> > thyra_mv = Thyra::createMembers(thyraA->range(),epetra_mv.NumVectors());
00156    getMapStrategy()->copyEpetraIntoThyra(epetra_mv,thyra_mv.ptr());
00157 
00158    // build preconObj_ 
00159    initPreconditioner(clear);
00160 
00161    // actually build the preconditioner
00162    preconFactory_->initializePrec(Thyra::defaultLinearOpSource(thyraA),thyra_mv,&*preconObj_,Thyra::SUPPORT_SOLVE_UNSPECIFIED);
00163    RCP<const Thyra::LinearOpBase<double> > preconditioner = preconObj_->getUnspecifiedPrecOp();
00164 
00165    SetOperator(preconditioner,false);
00166 
00167    firstBuildComplete_ = true;
00168 
00169    TEUCHOS_ASSERT(preconObj_!=Teuchos::null);
00170    TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
00171    TEUCHOS_ASSERT(getMapStrategy()!=Teuchos::null);
00172    TEUCHOS_ASSERT(firstBuildComplete_==true);
00173 }
00174 
00188 void EpetraBlockPreconditioner::rebuildPreconditioner(const Teuchos::RCP<const Epetra_Operator> & A)
00189 {
00190    Teko_DEBUG_SCOPE("EBP::rebuildPreconditioner",10);
00191 
00192    // if the preconditioner hasn't been built yet, rebuild from scratch
00193    if(not firstBuildComplete_) {
00194       buildPreconditioner(A,false);
00195       return;
00196    }
00197    Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00198 
00199    // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
00200    Teko_DEBUG_EXPR(timer.start(true));
00201    RCP<const Thyra::LinearOpBase<double> > thyraA = extractLinearOp(A);
00202    Teko_DEBUG_EXPR(timer.stop());
00203    Teko_DEBUG_MSG("EBP::rebuild get thyraop time =  " << timer.totalElapsedTime(),2);
00204 
00205    // reinitialize the preconditioner
00206    Teko_DEBUG_EXPR(timer.start(true));
00207    preconFactory_->initializePrec(Thyra::defaultLinearOpSource(thyraA),&*preconObj_,Thyra::SUPPORT_SOLVE_UNSPECIFIED);
00208    RCP<const Thyra::LinearOpBase<double> > preconditioner = preconObj_->getUnspecifiedPrecOp();
00209    Teko_DEBUG_EXPR(timer.stop());
00210    Teko_DEBUG_MSG("EBP::rebuild initialize prec time =  " << timer.totalElapsedTime(),2);
00211 
00212    Teko_DEBUG_EXPR(timer.start(true));
00213    SetOperator(preconditioner,false);
00214    Teko_DEBUG_EXPR(timer.stop());
00215    Teko_DEBUG_MSG("EBP::rebuild set operator time =  " << timer.totalElapsedTime(),2);
00216 
00217    TEUCHOS_ASSERT(preconObj_!=Teuchos::null);
00218    TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
00219    TEUCHOS_ASSERT(firstBuildComplete_==true);
00220 }
00221 
00235 void EpetraBlockPreconditioner::rebuildPreconditioner(const Teuchos::RCP<const Epetra_Operator> & A,const Epetra_MultiVector & epetra_mv)
00236 {
00237    Teko_DEBUG_SCOPE("EBP::rebuildPreconditioner - with solution",10);
00238 
00239    // if the preconditioner hasn't been built yet, rebuild from scratch
00240    if(not firstBuildComplete_) {
00241       buildPreconditioner(A,epetra_mv,false);
00242       return;
00243    }
00244    Teko_DEBUG_EXPR(Teuchos::Time timer(""));
00245 
00246    // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
00247    Teko_DEBUG_EXPR(timer.start(true));
00248    RCP<const Thyra::LinearOpBase<double> > thyraA = extractLinearOp(A);
00249    Teko_DEBUG_EXPR(timer.stop());
00250    Teko_DEBUG_MSG("EBP::rebuild get thyraop time =  " << timer.totalElapsedTime(),2);
00251 
00252    // build the thyra version of the source multivector
00253    Teko_DEBUG_EXPR(timer.start(true));
00254    RCP<Thyra::MultiVectorBase<double> > thyra_mv = Thyra::createMembers(thyraA->range(),epetra_mv.NumVectors());
00255    getMapStrategy()->copyEpetraIntoThyra(epetra_mv,thyra_mv.ptr());
00256    Teko_DEBUG_EXPR(timer.stop());
00257    Teko_DEBUG_MSG("EBP::rebuild vector copy time =  " << timer.totalElapsedTime(),2);
00258 
00259    // reinitialize the preconditioner
00260    Teko_DEBUG_EXPR(timer.start(true));
00261    preconFactory_->initializePrec(Thyra::defaultLinearOpSource(thyraA),thyra_mv,&*preconObj_,Thyra::SUPPORT_SOLVE_UNSPECIFIED);
00262    RCP<const Thyra::LinearOpBase<double> > preconditioner = preconObj_->getUnspecifiedPrecOp();
00263    Teko_DEBUG_EXPR(timer.stop());
00264    Teko_DEBUG_MSG("EBP::rebuild initialize prec time =  " << timer.totalElapsedTime(),2);
00265 
00266    Teko_DEBUG_EXPR(timer.start(true));
00267    SetOperator(preconditioner,false);
00268    Teko_DEBUG_EXPR(timer.stop());
00269    Teko_DEBUG_MSG("EBP::rebuild set operator time =  " << timer.totalElapsedTime(),2);
00270 
00271    TEUCHOS_ASSERT(preconObj_!=Teuchos::null);
00272    TEUCHOS_ASSERT(getThyraOp()!=Teuchos::null);
00273    TEUCHOS_ASSERT(firstBuildComplete_==true);
00274 }
00275 
00285 Teuchos::RCP<PreconditionerState> EpetraBlockPreconditioner::getPreconditionerState()
00286 {
00287    Teuchos::RCP<Preconditioner> bp = rcp_dynamic_cast<Preconditioner>(preconObj_);
00288 
00289    if(bp!=Teuchos::null)
00290       return bp->getStateObject();
00291  
00292    return Teuchos::null;
00293 }
00294 
00304 Teuchos::RCP<const PreconditionerState> EpetraBlockPreconditioner::getPreconditionerState() const
00305 {
00306    Teuchos::RCP<const Preconditioner> bp = rcp_dynamic_cast<const Preconditioner>(preconObj_);
00307 
00308    if(bp!=Teuchos::null)
00309       return bp->getStateObject();
00310  
00311    return Teuchos::null;
00312 }
00313 
00314 Teuchos::RCP<const Thyra::LinearOpBase<double> > EpetraBlockPreconditioner::extractLinearOp(const Teuchos::RCP<const Epetra_Operator> & A) const
00315 {
00316    // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
00317    const RCP<const EpetraOperatorWrapper> & eow = rcp_dynamic_cast<const EpetraOperatorWrapper>(A);
00318   
00319    // if it is an EpetraOperatorWrapper, then get the Thyra operator
00320    if(eow!=Teuchos::null)
00321       return eow->getThyraOp(); 
00322 
00323    // otherwise wrap it up as a thyra operator 
00324    return Thyra::epetraLinearOp(A);
00325 }
00326 
00327 Teuchos::RCP<const MappingStrategy> EpetraBlockPreconditioner::extractMappingStrategy(const Teuchos::RCP<const Epetra_Operator> & A) const
00328 {
00329    // extract EpetraOperatorWrapper (throw on failure) and corresponding thyra operator
00330    const RCP<const EpetraOperatorWrapper> & eow = rcp_dynamic_cast<const EpetraOperatorWrapper>(A);
00331   
00332    // if it is an EpetraOperatorWrapper, then get the Thyra operator
00333    if(eow!=Teuchos::null)
00334       return eow->getMapStrategy(); 
00335 
00336    // otherwise wrap it up as a thyra operator 
00337    RCP<const Epetra_Map> range = rcpFromRef(A->OperatorRangeMap());
00338    RCP<const Epetra_Map> domain = rcpFromRef(A->OperatorDomainMap());
00339    return rcp(new BasicMappingStrategy(range,domain,A->Comm()));
00340 }
00341 
00342 } // end namespace Epetra
00343 } // end namespace Teko
 All Classes Files Functions Variables