Belos Version of the Day
BelosBlockCGIter.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 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 #ifndef BELOS_BLOCK_CG_ITER_HPP
00043 #define BELOS_BLOCK_CG_ITER_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 #include "BelosCGIteration.hpp"
00052 
00053 #include "BelosLinearProblem.hpp"
00054 #include "BelosMatOrthoManager.hpp"
00055 #include "BelosOutputManager.hpp"
00056 #include "BelosStatusTest.hpp"
00057 #include "BelosOperatorTraits.hpp"
00058 #include "BelosMultiVecTraits.hpp"
00059 
00060 #include "Teuchos_BLAS.hpp"
00061 #include "Teuchos_LAPACK.hpp"
00062 #include "Teuchos_SerialDenseMatrix.hpp"
00063 #include "Teuchos_SerialDenseVector.hpp"
00064 #include "Teuchos_SerialSymDenseMatrix.hpp"
00065 #include "Teuchos_SerialSpdDenseSolver.hpp"
00066 #include "Teuchos_ScalarTraits.hpp"
00067 #include "Teuchos_ParameterList.hpp"
00068 #include "Teuchos_TimeMonitor.hpp"
00069 
00080 namespace Belos {
00081 
00082 template<class ScalarType, class MV, class OP>
00083 class BlockCGIter : virtual public CGIteration<ScalarType,MV,OP> {
00084 
00085   public:
00086 
00087   //
00088   // Convenience typedefs
00089   //
00090   typedef MultiVecTraits<ScalarType,MV> MVT;
00091   typedef MultiVecTraitsExt<ScalarType,MV> MVText;
00092   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00093   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00094   typedef typename SCT::magnitudeType MagnitudeType;
00095 
00097 
00098 
00104   BlockCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00105                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00106                const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00107                const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00108                Teuchos::ParameterList &params );
00109 
00111   virtual ~BlockCGIter() {};
00113 
00114 
00116 
00117 
00130   void iterate();
00131 
00146   void initializeCG(CGIterationState<ScalarType,MV> newstate);
00147 
00151   void initialize()
00152   {
00153     CGIterationState<ScalarType,MV> empty;
00154     initializeCG(empty);
00155   }
00156 
00163   CGIterationState<ScalarType,MV> getState() const {
00164     CGIterationState<ScalarType,MV> state;
00165     state.R = R_;
00166     state.P = P_;
00167     state.AP = AP_;
00168     state.Z = Z_;
00169     return state;
00170   }
00171 
00173 
00174 
00176 
00177 
00179   int getNumIters() const { return iter_; }
00180 
00182   void resetNumIters( int iter=0 ) { iter_ = iter; }
00183 
00186   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00187 
00189 
00191   Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00192 
00194 
00196 
00197 
00199   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00200 
00202   int getBlockSize() const { return blockSize_; }
00203 
00205   void setBlockSize(int blockSize);
00206 
00208   bool isInitialized() { return initialized_; }
00209 
00211 
00212   private:
00213 
00214   //
00215   // Internal methods
00216   //
00218   void setStateSize();
00219 
00220   //
00221   // Classes inputed through constructor that define the linear problem to be solved.
00222   //
00223   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00224   const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00225   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00226   const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00227 
00228   //
00229   // Algorithmic parameters
00230   //
00231   // blockSize_ is the solver block size.
00232   int blockSize_;
00233 
00234   //
00235   // Current solver state
00236   //
00237   // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00238   // is capable of running; _initialize is controlled  by the initialize() member method
00239   // For the implications of the state of initialized_, please see documentation for initialize()
00240   bool initialized_;
00241 
00242   // stateStorageInitialized_ specified that the state storage has be initialized.
00243   // This initialization may be postponed if the linear problem was generated without
00244   // the right-hand side or solution vectors.
00245   bool stateStorageInitialized_;
00246 
00247   // Current subspace dimension, and number of iterations performed.
00248   int iter_;
00249 
00250   //
00251   // State Storage
00252   //
00253   // Residual
00254   Teuchos::RCP<MV> R_;
00255   //
00256   // Preconditioned residual
00257   Teuchos::RCP<MV> Z_;
00258   //
00259   // Direction std::vector
00260   Teuchos::RCP<MV> P_;
00261   //
00262   // Operator applied to direction std::vector
00263   Teuchos::RCP<MV> AP_;
00264 
00265 };
00266 
00268   // Constructor.
00269   template<class ScalarType, class MV, class OP>
00270   BlockCGIter<ScalarType,MV,OP>::BlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
00271                                              const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00272                                              const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00273                                              const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00274                                              Teuchos::ParameterList &params ):
00275     lp_(problem),
00276     om_(printer),
00277     stest_(tester),
00278     ortho_(ortho),
00279     blockSize_(0),
00280     initialized_(false),
00281     stateStorageInitialized_(false),
00282     iter_(0)
00283   {
00284     // Set the block size and allocate data
00285     int bs = params.get("Block Size", 1);
00286     setBlockSize( bs );
00287   }
00288 
00290   // Setup the state storage.
00291   template <class ScalarType, class MV, class OP>
00292   void BlockCGIter<ScalarType,MV,OP>::setStateSize ()
00293   {
00294     if (!stateStorageInitialized_) {
00295 
00296       // Check if there is any multivector to clone from.
00297       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00298       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00299       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00300         stateStorageInitialized_ = false;
00301         return;
00302       }
00303       else {
00304 
00305         // Initialize the state storage
00306         // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00307         if (R_ == Teuchos::null || MVT::GetNumberVecs(*R_)!=blockSize_) {
00308           // Get the multivector that is not null.
00309           Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00310           TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00311                              "Belos::BlockCGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00312           R_ = MVT::Clone( *tmp, blockSize_ );
00313           Z_ = MVT::Clone( *tmp, blockSize_ );
00314           P_ = MVT::Clone( *tmp, blockSize_ );
00315           AP_ = MVT::Clone( *tmp, blockSize_ );
00316         }
00317 
00318         // State storage has now been initialized.
00319         stateStorageInitialized_ = true;
00320       }
00321     }
00322   }
00323 
00325   // Set the block size and make necessary adjustments.
00326   template <class ScalarType, class MV, class OP>
00327   void BlockCGIter<ScalarType,MV,OP>::setBlockSize (int blockSize)
00328   {
00329     // This routine only allocates space; it doesn't not perform any computation
00330     // any change in size will invalidate the state of the solver.
00331 
00332     TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument, "Belos::BlockGmresIter::setBlockSize was passed a non-positive argument.");
00333     if (blockSize == blockSize_) {
00334       // do nothing
00335       return;
00336     }
00337 
00338     if (blockSize!=blockSize_)
00339       stateStorageInitialized_ = false;
00340 
00341     blockSize_ = blockSize;
00342 
00343     initialized_ = false;
00344 
00345     // Use the current blockSize_ to initialize the state storage.
00346     setStateSize();
00347 
00348   }
00349 
00351   // Initialize this iteration object
00352   template <class ScalarType, class MV, class OP>
00353   void BlockCGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00354   {
00355     // Initialize the state storage if it isn't already.
00356     if (!stateStorageInitialized_)
00357       setStateSize();
00358 
00359     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00360                        "Belos::BlockCGIter::initialize(): Cannot initialize state storage!");
00361 
00362     // NOTE:  In BlockCGIter R_, the initial residual, is required!!!
00363     //
00364     std::string errstr("Belos::BlockCGIter::initialize(): Specified multivectors must have a consistent length and width.");
00365 
00366     // Create convenience variables for zero and one.
00367     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00368     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00369 
00370     if (newstate.R != Teuchos::null) {
00371 
00372       TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*newstate.R) != MVText::GetGlobalLength(*R_),
00373                           std::invalid_argument, errstr );
00374       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != blockSize_,
00375                           std::invalid_argument, errstr );
00376 
00377       // Copy basis vectors from newstate into V
00378       if (newstate.R != R_) {
00379         // copy over the initial residual (unpreconditioned).
00380         MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00381       }
00382       // Compute initial direction vectors
00383       // Initially, they are set to the preconditioned residuals
00384       //
00385       if ( lp_->getLeftPrec() != Teuchos::null ) {
00386         lp_->applyLeftPrec( *R_, *Z_ );
00387         if ( lp_->getRightPrec() != Teuchos::null ) {
00388           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, blockSize_ );
00389           lp_->applyRightPrec( *Z_, *tmp );
00390           Z_ = tmp;
00391         }
00392       }
00393       else if ( lp_->getRightPrec() != Teuchos::null ) {
00394         lp_->applyRightPrec( *R_, *Z_ );
00395       }
00396       else {
00397         Z_ = R_;
00398       }
00399       MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00400     }
00401     else {
00402 
00403       TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00404                          "Belos::BlockCGIter::initialize(): BlockCGStateIterState does not have initial residual.");
00405     }
00406 
00407     // The solver is initialized
00408     initialized_ = true;
00409   }
00410 
00411 
00413   // Iterate until the status test informs us we should stop.
00414   template <class ScalarType, class MV, class OP>
00415   void BlockCGIter<ScalarType,MV,OP>::iterate()
00416   {
00417     //
00418     // Allocate/initialize data structures
00419     //
00420     if (initialized_ == false) {
00421       initialize();
00422     }
00423     // Allocate data needed for LAPACK work.
00424     int info = 0;
00425     //char UPLO = 'U';
00426     //(void) UPLO; // silence "unused variable" compiler warnings
00427     bool uplo = true;
00428     Teuchos::LAPACK<int,ScalarType> lapack;
00429 
00430     // Allocate memory for scalars.
00431     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( blockSize_, blockSize_ );
00432     Teuchos::SerialDenseMatrix<int,ScalarType> beta( blockSize_, blockSize_ );
00433     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( blockSize_, blockSize_ ),
00434       rHz_old( blockSize_, blockSize_ ), pAp( blockSize_, blockSize_ );
00435     Teuchos::SerialSymDenseMatrix<int,ScalarType> pApHerm(Teuchos::View, uplo, pAp.values(), blockSize_, blockSize_);
00436 
00437     // Create dense spd solver.
00438     Teuchos::SerialSpdDenseSolver<int,ScalarType> lltSolver;
00439 
00440     // Create convenience variables for zero and one.
00441     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00442 
00443     // Get the current solution std::vector.
00444     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00445 
00446     // Check that the current solution std::vector has blockSize_ columns.
00447     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != blockSize_, CGIterateFailure,
00448                         "Belos::BlockCGIter::iterate(): current linear system does not have the right number of vectors!" );
00449     int rank = ortho_->normalize( *P_, Teuchos::null );
00450     TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,CGIterationOrthoFailure,
00451                          "Belos::BlockCGIter::iterate(): Failed to compute initial block of orthonormal direction vectors.");
00452 
00453 
00455     // Iterate until the status test tells us to stop.
00456     //
00457     while (stest_->checkStatus(this) != Passed) {
00458 
00459       // Increment the iteration
00460       iter_++;
00461 
00462       // Multiply the current direction std::vector by A and store in Ap_
00463       lp_->applyOp( *P_, *AP_ );
00464 
00465       // Compute alpha := <P_,R_> / <P_,AP_>
00466       // 1) Compute P^T * A * P = pAp and P^T * R
00467       // 2) Compute the Cholesky Factorization of pAp
00468       // 3) Back and forward solves to compute alpha
00469       //
00470       MVT::MvTransMv( one, *P_, *R_, alpha );
00471       MVT::MvTransMv( one, *P_, *AP_, pAp );
00472 
00473       // Compute Cholesky factorization of pAp
00474       lltSolver.setMatrix( Teuchos::rcp(&pApHerm, false) );
00475       lltSolver.factorWithEquilibration( true );
00476       info = lltSolver.factor();
00477       TEUCHOS_TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
00478                          "Belos::BlockCGIter::iterate(): Failed to compute Cholesky factorization using LAPACK routine POTRF.");
00479 
00480       // Compute alpha by performing a back and forward solve with the Cholesky factorization in pAp.
00481       lltSolver.setVectors( Teuchos::rcp( &alpha, false ), Teuchos::rcp( &alpha, false ) );
00482       info = lltSolver.solve();
00483       TEUCHOS_TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
00484                          "Belos::BlockCGIter::iterate(): Failed to compute alpha using Cholesky factorization (POTRS).");
00485 
00486       //
00487       // Update the solution std::vector X := X + alpha * P_
00488       //
00489       MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
00490       lp_->updateSolution();
00491       //
00492       // Compute the new residual R_ := R_ - alpha * AP_
00493       //
00494       MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
00495       //
00496       // Compute the new preconditioned residual, Z_.
00497       if ( lp_->getLeftPrec() != Teuchos::null ) {
00498         lp_->applyLeftPrec( *R_, *Z_ );
00499         if ( lp_->getRightPrec() != Teuchos::null ) {
00500           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, blockSize_ );
00501           lp_->applyRightPrec( *Z_, *tmp );
00502           Z_ = tmp;
00503         }
00504       }
00505       else if ( lp_->getRightPrec() != Teuchos::null ) {
00506         lp_->applyRightPrec( *R_, *Z_ );
00507       }
00508       else {
00509         Z_ = R_;
00510       }
00511       //
00512       // Compute beta := <AP_,Z_> / <P_,AP_>
00513       // 1) Compute AP_^T * Z_
00514       // 2) Compute the Cholesky Factorization of pAp (already have)
00515       // 3) Back and forward solves to compute beta
00516 
00517       // Compute <AP_,Z>
00518       MVT::MvTransMv( -one, *AP_, *Z_, beta );
00519       //
00520       lltSolver.setVectors( Teuchos::rcp( &beta, false ), Teuchos::rcp( &beta, false ) );
00521       info = lltSolver.solve();
00522       TEUCHOS_TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
00523                          "Belos::BlockCGIter::iterate(): Failed to compute beta using Cholesky factorization (POTRS).");
00524       //
00525       // Compute the new direction vectors P_ = Z_ + P_ * beta
00526       //
00527       Teuchos::RCP<MV> Pnew = MVT::CloneCopy( *Z_ );
00528       MVT::MvTimesMatAddMv(one, *P_, beta, one, *Pnew);
00529       P_ = Pnew;
00530 
00531       // Compute orthonormal block of new direction vectors.
00532       rank = ortho_->normalize( *P_, Teuchos::null );
00533       TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,CGIterationOrthoFailure,
00534                          "Belos::BlockCGIter::iterate(): Failed to compute block of orthonormal direction vectors.");
00535 
00536     } // end while (sTest_->checkStatus(this) != Passed)
00537   }
00538 
00539 } // end Belos namespace
00540 
00541 #endif /* BELOS_BLOCK_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines