Belos Version of the Day
BelosLSQRIter.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_LSQR_ITER_HPP
00043 #define BELOS_LSQR_ITER_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 #include "BelosLSQRIteration.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_SerialDenseMatrix.hpp"
00062 #include "Teuchos_SerialDenseVector.hpp"
00063 #include "Teuchos_ScalarTraits.hpp"
00064 #include "Teuchos_ParameterList.hpp"
00065 #include "Teuchos_TimeMonitor.hpp"
00066 
00078 namespace Belos {
00079 
00080 template<class ScalarType, class MV, class OP>
00081 class LSQRIter : virtual public Belos::Iteration<ScalarType,MV,OP> {
00082 
00083   public:
00084     
00085   //
00086   // Convenience typedefs
00087   //
00088   typedef Belos::MultiVecTraits<ScalarType,MV> MVT;
00089   typedef Belos::OperatorTraits<ScalarType,MV,OP> OPT;
00090   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00091   typedef typename SCT::magnitudeType MagnitudeType;
00092 
00094 
00095 
00102   LSQRIter( const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> > &problem, 
00103       const Teuchos::RCP<Belos::OutputManager<ScalarType> > &printer,
00104       const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> > &tester,
00105       Teuchos::ParameterList &params );
00106 
00107 //                  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00108 
00109 
00111   virtual ~LSQRIter() {};
00113 
00114 
00116 
00117   
00130   void iterate();
00131 
00141   void initializeLSQR(LSQRIterationState<ScalarType,MV> newstate);
00142 
00145   void initialize()
00146   {
00147     LSQRIterationState<ScalarType,MV> empty;
00148     initializeLSQR(empty);
00149   }
00150   
00158   LSQRIterationState<ScalarType,MV> getState() const {
00159     LSQRIterationState<ScalarType,MV> state;
00160     state.U = U_;
00161     state.V = V_;
00162     state.W = W_;
00163     state.lambda = lambda_;
00164     state.resid_norm = resid_norm_;
00165     state.frob_mat_norm = frob_mat_norm_;
00166     state.mat_cond_num = mat_cond_num_;
00167     state.mat_resid_norm = mat_resid_norm_; 
00168     state.sol_norm = sol_norm_;
00169     state.bnorm = bnorm_;
00170     return state;
00171   }
00172 
00174 
00175   
00177 
00178 
00180   int getNumIters() const { return iter_; }
00181   
00183   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00184 
00187   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return Teuchos::null; }
00188 
00190 
00192   Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00193 
00195   
00197 
00198 
00200   const Belos::LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00201 
00203   int getBlockSize() const { return 1; }
00204 
00206   //This is unique to single vector methods.
00207   void setBlockSize(int blockSize) {
00208     TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00209            "LSQRIter::setBlockSize(): Cannot use a block size that is not one.");
00210   }
00211 
00213   bool isInitialized() { return initialized_; }
00214 
00216 
00217   private:
00218 
00219   //
00220   // Internal methods
00221   //
00223   void setStateSize();
00224   
00225   //
00226   // Classes inputed through constructor that define the linear problem to be solved.
00227   //
00228   const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> >    lp_;
00229   const Teuchos::RCP<Belos::OutputManager<ScalarType> >          om_;
00230   const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> >       stest_;
00231 //  const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00232 
00233   //  
00234   // Current solver state
00235   //
00236   // initialized_ specifies that the basis vectors have been initialized and the iterate()
00237   // routine is capable of running; _initialize is controlled  by the initialize() member
00238   // method.  For the implications of the state of initialized_, please see documentation
00239   // for initialize()
00240   bool initialized_;
00241 
00242   // stateStorageInitialized_ specifies that the state storage has been 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 number of iterations performed.
00248   int iter_;
00249   
00250   // 
00251   // State Storage
00252   // 
00253   //
00254   // Bidiagonalization vector
00255   Teuchos::RCP<MV> U_;
00256   //
00257   // Bidiagonalization vector
00258   Teuchos::RCP<MV> V_;
00259   //
00260   // Direction vector
00261   Teuchos::RCP<MV> W_;
00262   //
00263   // Damping value
00264   ScalarType lambda_;
00265   //
00266   // Residual norm estimate
00267   ScalarType resid_norm_;
00268   //
00269   // Frobenius norm estimate
00270   ScalarType frob_mat_norm_;
00271   //
00272   // Condition number estimate
00273   ScalarType mat_cond_num_;
00274   //
00275   // A^T*resid norm estimate
00276   ScalarType mat_resid_norm_;
00277   //
00278   // Solution norm estimate
00279   ScalarType sol_norm_;
00280   //
00281   // RHS norm
00282   ScalarType bnorm_;
00283 
00284 };
00285 
00287   // Constructor.
00288 
00289 //                                     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00290 
00291   template<class ScalarType, class MV, class OP>
00292   LSQRIter<ScalarType,MV,OP>::LSQRIter(const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> > &problem, 
00293                const Teuchos::RCP<Belos::OutputManager<ScalarType> > &printer,
00294                const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> > &tester,
00295                Teuchos::ParameterList &params ):
00296     lp_(problem),
00297     om_(printer),
00298     stest_(tester),
00299     initialized_(false),
00300     stateStorageInitialized_(false),
00301     iter_(0),
00302     lambda_(params.get("Lambda", 0.0))
00303   {
00304   }
00305 
00307   // Setup the state storage.
00308   template <class ScalarType, class MV, class OP>
00309   void LSQRIter<ScalarType,MV,OP>::setStateSize ()
00310   {
00311     if (!stateStorageInitialized_) {
00312       // Check if there is any multivector to clone from.
00313       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00314       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00315       if (lhsMV == Teuchos::null || rhsMV == Teuchos::null) {
00316   stateStorageInitialized_ = false;
00317   return;
00318       }
00319       else {
00320   // Initialize the state storage
00321   // If the subspace has not been initialized before, generate it
00322         // using the LHS and RHS from lp_.
00323   if (U_ == Teuchos::null) {
00324     // Get the multivectors.
00325     TEST_FOR_EXCEPTION(rhsMV == Teuchos::null, std::invalid_argument, "LSQRIter::setStateSize(): linear problem does not specify right hand multivector to clone from.");
00326     TEST_FOR_EXCEPTION(lhsMV == Teuchos::null, std::invalid_argument, "LSQRIter::setStateSize(): linear problem does not specify left hand multivector to clone from.");
00327 
00328     U_ = MVT::Clone( *rhsMV, 1 );
00329     V_ = MVT::Clone( *lhsMV, 1 );
00330     W_ = MVT::Clone( *lhsMV, 1 );
00331   }
00332   
00333   // State storage has now been initialized.
00334   stateStorageInitialized_ = true;
00335       }
00336     }
00337   }
00338 
00339 
00341   // Initialize this iteration object
00342   template <class ScalarType, class MV, class OP>
00343   void LSQRIter<ScalarType,MV,OP>::initializeLSQR(LSQRIterationState<ScalarType,MV> newstate)
00344   {
00345     // Initialize the state storage if it isn't already.
00346     if (!stateStorageInitialized_) 
00347       setStateSize();
00348 
00349     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00350            "LSQRIter::initialize(): Cannot initialize state storage!");
00351     
00352     std::string errstr("LSQRIter::initialize(): Specified multivectors must have a consistent length and width.");
00353 
00354     // Create convenience variables for zero and one.
00355     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00356     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00357 
00358     // Compute initial bidiagonalization vectors and search direction
00359     //
00360     Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00361     Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00362 
00363     OPT::Apply(*(lp_->getOperator()), *lhsMV, *U_);
00364     MVT::MvAddMv( one, *rhsMV, -one, *U_, *U_);
00365     OPT::Apply(*(lp_->getOperator()), *U_, *V_, CONJTRANS);
00366     MVT::MvAddMv( one, *V_, zero, *V_, *W_);
00367     
00368     frob_mat_norm_ = zero;
00369     mat_cond_num_ = zero;
00370     sol_norm_ = zero;
00371     
00372     // The solver is initialized
00373     initialized_ = true;
00374   }
00375 
00376 
00378   // Iterate until the status test informs us we should stop.
00379   template <class ScalarType, class MV, class OP>
00380   void LSQRIter<ScalarType,MV,OP>::iterate()
00381   {
00382     //
00383     // Allocate/initialize data structures
00384     //
00385     if (initialized_ == false) {
00386       initialize();
00387     }
00388     
00389     // Create convenience variables for zero and one.
00390     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00391     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00392 
00393     // Allocate memory for scalars.
00394     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> alpha(1);
00395     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> beta(1);
00396     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> wnorm2(1);
00397     ScalarType rhobar, phibar, cs1, phi, rho, cs, sn, theta, xxnorm = zero, common;
00398     ScalarType zetabar, sn1, psi, res = zero, bbnorm = zero, ddnorm = zero, gamma, tau;
00399     ScalarType cs2 = -one, sn2 = zero, gammabar, zeta = zero, delta;
00400     
00401     // Allocate memory for working vectors.
00402     // Operator applied to bidiagonalization vector
00403     Teuchos::RCP<MV> AV;
00404     // Transpose of operator applied to bidiagonalization vector
00405     Teuchos::RCP<MV> AtU;
00406     AV = MVT::Clone( *U_, 1);
00407     AtU = MVT::Clone( *V_, 1);
00408 
00409     // Get the current solution vector.
00410     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00411 
00412     // Check that the current solution vector only has one column. 
00413     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, LSQRIterateFailure,
00414                         "LSQRIter::iterate(): current linear system has more than one vector!" );
00415 
00416     // Compute alpha and beta and scale bidiagonalization vectors
00417     MVT::MvNorm( *U_, beta );
00418     MVT::MvScale( *U_, one / beta[0] );
00419     MVT::MvScale( *V_, one / beta[0] );
00420     MVT::MvNorm( *V_, alpha );
00421     MVT::MvScale( *V_, one / alpha[0] );
00422     MVT::MvScale( *W_, one / (beta[0] * alpha[0]) );
00423 
00424     rhobar = alpha[0];
00425     phibar = beta[0];
00426     
00427     resid_norm_ = beta[0];
00428     mat_resid_norm_ = alpha[0] * beta[0];
00429     bnorm_ = beta[0];
00430 
00432     // Iterate until the status test tells us to stop.
00433     //
00434     while (stest_->checkStatus(this) != Belos::Passed) {
00435       // Increment the iteration
00436       iter_++;
00437 
00438       // Perform the next step of the bidiagonalization to obtain the next U_ and V_ vectors
00439       // and the scalars alpha and beta.  These satisfy the relations
00440       // beta*U_ = AV - alpha*U_
00441       // alpha*V_ = AtU - beta*V_
00442 
00443       // AV := A * V_
00444       OPT::Apply(*(lp_->getOperator()), *V_, *AV);
00445       MVT::MvAddMv( one, *AV, -alpha[0], *U_, *U_ );
00446       MVT::MvNorm( *U_, beta);
00447       // Check that beta is a positive number!
00448       TEST_FOR_EXCEPTION( SCT::real(beta[0]) <= zero, LSQRIterateFailure, "LSQRIter::iterate(): non-positive value for beta encountered!");
00449       bbnorm += alpha[0]*alpha[0] + beta[0]*beta[0] + lambda_*lambda_;
00450       MVT::MvScale( *U_, one / beta[0] );
00451 
00452       // AtU := A^T * U_
00453       OPT::Apply(*(lp_->getOperator()), *U_, *AtU, CONJTRANS);
00454       MVT::MvAddMv( one, *AtU, -beta[0], *V_, *V_ );
00455       MVT::MvNorm( *V_, alpha );
00456       // Check that alpha is a positive number!
00457       TEST_FOR_EXCEPTION( SCT::real(alpha[0]) <= zero, LSQRIterateFailure, "LSQRIter::iterate(): non-positive value for alpha encountered!");
00458       MVT::MvScale( *V_, one / alpha[0] );
00459 
00460       // Use a plane rotation to eliminate the damping parameter.  
00461       // This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
00462       common = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_);
00463       cs1 = rhobar / common;
00464       sn1 = lambda_ / common;
00465       psi = sn1 * phibar;
00466       phibar = cs1 * phibar;
00467 
00468       // Use a plane rotation to eliminate the subdiagonal element (beta)
00469       // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
00470       rho = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_ + beta[0]*beta[0]);
00471       cs = common / rho; 
00472       sn = beta[0] / rho;
00473       theta = sn * alpha[0];
00474       rhobar = -cs * alpha[0];
00475       phi = cs * phibar;
00476       phibar = sn * phibar;
00477       tau = sn * phi;
00478 
00479       delta = sn2 * rho;
00480       gammabar = -cs2 * rho;
00481       zetabar = (phi - delta*zeta) / gammabar;
00482       sol_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(xxnorm + zetabar*zetabar);
00483       gamma = Teuchos::ScalarTraits< ScalarType >::squareroot(gammabar*gammabar + theta*theta);
00484       cs2 = gammabar / gamma;
00485       sn2 = theta / gamma;
00486       zeta = (phi - delta*zeta) / gamma;
00487       xxnorm += zeta*zeta;
00488 
00489       // Update the solution vector and search direction vector
00490       MVT::MvAddMv( phi / rho, *W_, one, *cur_soln_vec, *cur_soln_vec);
00491       lp_->updateSolution();
00492       MVT::MvNorm( *W_, wnorm2 );
00493       ddnorm += (one / rho)*(one / rho) * wnorm2[0]*wnorm2[0];
00494       MVT::MvAddMv( one, *V_, -theta / rho, *W_, *W_ );
00495 
00496       frob_mat_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(bbnorm);
00497       mat_cond_num_ = frob_mat_norm_ * Teuchos::ScalarTraits< ScalarType >::squareroot(ddnorm);
00498       res+= psi*psi;
00499       resid_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(phibar*phibar + res);
00500       mat_resid_norm_ = alpha[0] * Teuchos::ScalarTraits< ScalarType >::magnitude(tau);
00501 
00502     } // end while (sTest_->checkStatus(this) != Passed)
00503   } // iterate()
00504 
00505 } // end Belos namespace
00506 
00507 #endif /* BELOS_LSQR_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines