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 "BelosOutputManager.hpp"
00055 #include "BelosStatusTest.hpp"
00056 #include "BelosOperatorTraits.hpp"
00057 #include "BelosMultiVecTraits.hpp"
00058 //#include "BelosMatOrthoManager.hpp"  (needed for blocks)
00059 
00060 //#include "Teuchos_BLAS.hpp"  (needed for blocks)
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 
00074 namespace Belos {
00075 
00076 template<class ScalarType, class MV, class OP>
00077 class LSQRIter : virtual public Belos::Iteration<ScalarType,MV,OP> {
00078 
00079   public:
00080     
00081   //
00082   // Convenience typedefs
00083   //
00084   typedef Belos::MultiVecTraits<ScalarType,MV> MVT;
00085   typedef Belos::OperatorTraits<ScalarType,MV,OP> OPT;
00086   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00087   typedef typename SCT::magnitudeType MagnitudeType;
00088 
00090 
00091 
00098   LSQRIter( const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> > &problem, 
00099       const Teuchos::RCP<Belos::OutputManager<ScalarType> > &printer,
00100       const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> > &tester,
00101       Teuchos::ParameterList &params );
00102 
00103 // If either blocks or reorthogonalization exist, then 
00104 //                  const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00105 
00106 
00108   virtual ~LSQRIter() {};
00110 
00111 
00113 
00114   
00126   void iterate();
00127 
00137   void initializeLSQR(LSQRIterationState<ScalarType,MV> newstate);
00138 
00141   void initialize()
00142   {
00143     LSQRIterationState<ScalarType,MV> empty;
00144     initializeLSQR(empty);
00145   }
00146   
00154   LSQRIterationState<ScalarType,MV> getState() const {
00155     LSQRIterationState<ScalarType,MV> state;
00156     state.U = U_;  // right Lanczos vector
00157     state.V = V_;  // left  Lanczos vector
00158     state.W = W_;  // OP * V 
00159     state.lambda = lambda_;  
00160     state.resid_norm = resid_norm_;
00161     state.frob_mat_norm = frob_mat_norm_;
00162     state.mat_cond_num = mat_cond_num_;
00163     state.mat_resid_norm = mat_resid_norm_; 
00164     state.sol_norm = sol_norm_;
00165     state.bnorm = bnorm_;
00166     return state;
00167   }
00168 
00170 
00171   
00173 
00174 
00176   int getNumIters() const { return iter_; }
00177   
00179   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00180 
00183   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return Teuchos::null; }
00184 
00186 
00188   Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00189 
00191   
00193 
00194 
00196   const Belos::LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00197 
00199   int getBlockSize() const { return 1; }
00200 
00201 
00203   //This is unique to single vector methods.
00204   void setBlockSize(int blockSize) {
00205     TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00206            "LSQRIter::setBlockSize(): Cannot use a block size that is not one.");
00207   }
00208 
00210   bool isInitialized() { return initialized_; }
00211 
00213 
00214   private:
00215 
00216   //
00217   // Internal methods
00218   //
00220   void setStateSize();
00221   
00222   //
00223   // Classes inputed through constructor that define the linear problem to be solved.
00224   //
00225   const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> >    lp_;
00226   const Teuchos::RCP<Belos::OutputManager<ScalarType> >          om_;
00227   const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> >       stest_;
00228 //  const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00229 
00230   //  
00231   // Current solver state
00232   //
00233   // initialized_ specifies that the basis vectors have been initialized and the iterate()
00234   // routine is capable of running; _initialize is controlled  by the initialize() member
00235   // method.  For the implications of the state of initialized_, please see documentation
00236   // for initialize()
00237   bool initialized_;
00238 
00239   // stateStorageInitialized_ specifies that the state storage has been initialized.
00240   // This initialization may be postponed if the linear problem was generated without 
00241   // the right-hand side or solution vectors.
00242   bool stateStorageInitialized_;
00243 
00244   // Current number of iterations performed.
00245   int iter_;
00246   
00247   // 
00248   // State Storage
00249   // 
00250   //
00251   // Bidiagonalization vector
00252   Teuchos::RCP<MV> U_;
00253   //
00254   // Bidiagonalization vector
00255   Teuchos::RCP<MV> V_;
00256   //
00257   // Direction vector
00258   Teuchos::RCP<MV> W_;
00259   //
00260   // Damping value
00261   MagnitudeType lambda_;
00262   //
00263   // Residual norm estimate
00264   ScalarType resid_norm_;
00265   //
00266   // Frobenius norm estimate
00267   ScalarType frob_mat_norm_;
00268   //
00269   // Condition number estimate
00270   ScalarType mat_cond_num_;
00271   //
00272   // A^T*resid norm estimate
00273   ScalarType mat_resid_norm_;
00274   //
00275   // Solution norm estimate
00276   ScalarType sol_norm_;
00277   //
00278   // RHS norm
00279   ScalarType bnorm_;
00280 
00281 };
00282 
00284   // Constructor.
00285 
00286 //                                     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00287 
00288   template<class ScalarType, class MV, class OP>
00289   LSQRIter<ScalarType,MV,OP>::LSQRIter(const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> > &problem, 
00290                const Teuchos::RCP<Belos::OutputManager<ScalarType> > &printer,
00291                const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> > &tester,
00292                Teuchos::ParameterList &params ):
00293     lp_(problem),
00294     om_(printer),
00295     stest_(tester),
00296     initialized_(false),
00297     stateStorageInitialized_(false),
00298     iter_(0),
00299     lambda_(params.get<MagnitudeType> ("Lambda"))
00300   {
00301   }
00302 
00304   // Setup the state storage.
00305   template <class ScalarType, class MV, class OP>
00306   void LSQRIter<ScalarType,MV,OP>::setStateSize ()
00307   {
00308     if (!stateStorageInitialized_) {
00309       // Check if there is any multivector to clone from.
00310       Teuchos::RCP<const MV> rhsMV = lp_->getInitPrecResVec();
00311       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00312       if (lhsMV == Teuchos::null || rhsMV == Teuchos::null) {
00313   stateStorageInitialized_ = false;
00314   return;
00315       }
00316       else {
00317   // Initialize the state storage
00318   // If the subspace has not been initialized before, generate it
00319         // using the LHS and RHS from lp_.
00320   if (U_ == Teuchos::null) {
00321     // Get the multivectors.
00322     TEUCHOS_TEST_FOR_EXCEPTION(rhsMV == Teuchos::null, std::invalid_argument, "LSQRIter::setStateSize(): linear problem does not specify right hand multivector to clone from.");
00323     TEUCHOS_TEST_FOR_EXCEPTION(lhsMV == Teuchos::null, std::invalid_argument, "LSQRIter::setStateSize(): linear problem does not specify left hand multivector to clone from.");
00324 
00325     U_ = MVT::Clone( *rhsMV, 1 ); // LeftPrecond * rhs
00326     V_ = MVT::Clone( *lhsMV, 1 ); // zero, overwrittein in 
00327     W_ = MVT::Clone( *lhsMV, 1 ); // zero, initializeLSQR
00328   }
00329   
00330   // State storage has now been initialized.
00331   stateStorageInitialized_ = true;
00332       }
00333     }
00334   }
00335 
00336 
00338   // Initialize this iteration object
00339   template <class ScalarType, class MV, class OP>
00340   void LSQRIter<ScalarType,MV,OP>::initializeLSQR(LSQRIterationState<ScalarType,MV> newstate)
00341   {
00342     using Teuchos::RCP;
00343 
00344     // Initialize the state storage if it isn't already.
00345     if (!stateStorageInitialized_) 
00346       setStateSize();
00347 
00348     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00349            "LSQRIter::initialize(): Cannot initialize state storage!");
00350     
00351     std::string errstr("LSQRIter::initialize(): Specified multivectors must have a consistent length and width.");
00352 
00353 
00354     // Compute initial bidiagonalization vectors and search direction
00355     //
00356     RCP<const MV> lhsMV = lp_->getLHS(); // contains initial guess,
00357 
00358     bool debugSerialLSQR = false;
00359 
00360     if( debugSerialLSQR )
00361       {
00362         std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> lhsNorm(1);
00363         MVT::MvNorm( *lhsMV, lhsNorm );
00364         std::cout << "initializeLSQR lhsNorm " << lhsNorm[0] << std::endl;
00365       }
00366 
00367     // LinearProlbem provides right-hand side vectors including RHS CurrRHSVec InitResVec.
00368     RCP<const MV> rhsMV = lp_->getInitPrecResVec(); 
00369     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00370     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00371     MVT::MvAddMv( one, *rhsMV, zero, *U_, *U_);
00372 
00373     RCP<const OP> M_left = lp_->getLeftPrec();
00374     RCP<const OP> A = lp_->getOperator();
00375     RCP<const OP> M_right = lp_->getRightPrec();
00376 
00377     if( debugSerialLSQR )
00378       {
00379         std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> rhsNorm(1); 
00380         MVT::MvNorm( *U_, rhsNorm );
00381         std::cout << "initializeLSQR | U_ | : " << rhsNorm[0] << std::endl;
00382         //U_->Print(std::cout);
00383       }
00384 
00385     //MVT::MvScale( *V_, zero );
00386 
00387     // Apply the (conjugate) transpose of the preconditioned operator:
00388     //
00389     // V := (M_L A M_R)^* U, which means
00390     // V := M_R^* (A^* (M_L^* U)).
00391     //
00392     //OPT::Apply(*(lp_->getOperator()), *U_, *V_, CONJTRANS);
00393     if ( M_left.is_null())
00394       {
00395         OPT::Apply (*A, *U_, *V_, CONJTRANS); // V_ = A' U_
00396         //std::cout << "***************  V_ ****************" << std::endl;
00397         //V_->Print(std::cout);
00398       }
00399     else
00400       {
00401         RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
00402         OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS); 
00403         OPT::Apply (*A, *tempInRangeOfA, *V_, CONJTRANS); // V_ = A' LeftPrec' U_
00404         //std::cout << "mLeft  V_ = " << *V_ << std::endl;
00405       }
00406     if (! M_right.is_null())
00407       {
00408         RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
00409 
00410         OPT::Apply (*M_right, *tempInDomainOfA, *V_, CONJTRANS); // V:= RtPrec' A' LeftPrec' U
00411         //std::cout << "mRight  V_ = " << *V_ << std::endl;
00412       }
00413 
00414     // W := V (copy the vector) 
00415     MVT::MvAddMv( one, *V_, zero, *V_, *W_);
00416     
00417     frob_mat_norm_ = zero; // These are 
00418     mat_cond_num_ = one;   // lower
00419     sol_norm_ = zero;      // bounds. 
00420     
00421     // The solver is initialized.
00422     initialized_ = true;
00423   }
00424 
00425 
00427   // Iterate until the status test informs us we should stop.
00428   template <class ScalarType, class MV, class OP>
00429   void LSQRIter<ScalarType,MV,OP>::iterate()
00430   {
00431     //
00432     // Allocate/initialize data structures
00433     //
00434     if (initialized_ == false) {
00435       initialize();
00436     }
00437     
00438     // Create convenience variables for zero and one.
00439     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00440     const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00441     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00442 
00443     // Allocate memory for scalars.
00444     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> alpha(1);
00445     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> beta(1);
00446     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> xi(1);
00447     // xi is a dumb scalar used for storing inner products. 
00448     // Eventually SDM will replace the "vectors".
00449     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> wnorm2(1); 
00450     ScalarType rhobar, phibar, cs1, phi, rho, cs, sn, theta, xxnorm = MTzero, common;
00451     ScalarType zetabar, sn1, psi, res = zero, ddnorm = zero, gamma, tau;
00452     ScalarType anorm2 = zero;
00453     ScalarType cs2 = -one, sn2 = zero, gammabar, zeta = zero, delta;
00454     
00455     // The pair of work vectors AV and AtU are 
00456     Teuchos::RCP<MV> AV; // used in applying A to V_ and
00457     AV = MVT::Clone( *U_, 1);
00458     Teuchos::RCP<MV> AtU; // used in applying A^TRANS to U_ respectively.
00459     AtU = MVT::Clone( *V_, 1);
00460     bool debugSerialLSQR = false;
00461 
00462     // Get the current solution vector.
00463     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00464 
00465 
00466     // Check that the current solution vector only has one column. 
00467     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, LSQRIterateFailure,
00468                         "LSQRIter::iterate(): current linear system has more than one vector!" );
00469 
00470     // In initializeLSQR among other things V = A' U.
00471     // alpha and beta normalize these vectors.
00472     MVT::MvNorm( *U_, beta );
00473     if( SCT::real(beta[0]) >  zero  )
00474       {
00475         MVT::MvScale( *U_, one / beta[0] );
00476 
00477         //std::cout << "***************  U/beta ****************" << std::endl;
00478         //U_->Print(std::cout);
00479 
00480         MVT::MvScale( *V_, one / beta[0] );  // scale V = A'U to normalize U
00481 
00482         //std::cout << "***************  V/beta ****************" << std::endl;
00483         //V_->Print(std::cout);
00484       }
00485     MVT::MvNorm( *V_, alpha );
00486     if( debugSerialLSQR )
00487       {
00488          // used to compare with implementations 
00489          // initializing mat_resid_norm to alpha/beta
00490          std::cout << iter_ << " First alpha " << alpha[0]   << " beta " << beta[0] << " lambda " << lambda_ << std::endl;
00491       }
00492     if( SCT::real(alpha[0]) >  zero  )
00493       {
00494         MVT::MvScale( *V_, one / alpha[0] ); // V alpha = A' U to normalize V
00495       }
00496     if( beta[0] * alpha[0] >  zero  )
00497       { 
00498         MVT::MvScale( *W_, one / (beta[0] * alpha[0]) ); // W = V
00499       }
00500     else
00501       { 
00502         MVT::MvScale( *W_, zero  ); 
00503       }
00504 
00505     using Teuchos::RCP;
00506     RCP<const OP> M_left = lp_->getLeftPrec();
00507     RCP<const OP> A = lp_->getOperator();
00508     RCP<const OP> M_right = lp_->getRightPrec();
00509 
00510     rhobar = alpha[0];
00511     phibar = beta[0];
00512     bnorm_ = beta[0];
00513     resid_norm_ = beta[0];
00514     mat_resid_norm_ = alpha[0] * beta[0];
00515 
00516 
00518     // Iterate until the status test tells us to stop.
00519     //
00520     while (stest_->checkStatus(this) != Belos::Passed) {
00521       // Increment the iteration
00522       iter_++;
00523 
00524       // Perform the next step of the bidiagonalization.
00525       // The next U_ and V_ vectors and scalars alpha and beta satisfy 
00526       // U_ betaNew := AV - U_ alphaOld ...
00527 
00528       if ( M_right.is_null() )
00529         {
00530           OPT::Apply(*A, *V_, *AV); // AV := A * V_
00531         }
00532       else
00533         {
00534           RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
00535           OPT::Apply (*M_right, *V_, *tempInDomainOfA); 
00536           OPT::Apply(*A, *tempInDomainOfA, *AV); 
00537         }
00538 
00539       if (! M_left.is_null())
00540         {
00541           RCP<MV> tempInRangeOfA = MVT::CloneCopy (*AV);
00542           OPT::Apply (*M_left, *tempInRangeOfA, *AV); // AV may change
00543         }
00544 
00545 
00546       if ( !( M_left.is_null()  &&  M_right.is_null() )
00547            && debugSerialLSQR && iter_ == 1)
00548         {
00549           // In practice, LSQR may reveal bugs in transposed preconditioners.
00550           // This is the test that catches this type of bug.
00551           // 1. confirm that V alpha = A' U
00552 
00553           if (! M_left.is_null())
00554             {
00555               RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
00556               OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS);
00557               OPT::Apply (*A, *tempInRangeOfA, *AtU, CONJTRANS);   // AtU = B'L'U
00558             }
00559           else
00560             {
00561               OPT::Apply (*A, *U_, *AtU, CONJTRANS);   // AtU = B'U
00562             }
00563           if ( !( M_right.is_null() ) )
00564             {
00565               RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
00566               OPT::Apply (*M_right, *tempInDomainOfA, *AtU, CONJTRANS); // AtU := R' AtU
00567             }
00568 
00569           MVT::MvAddMv( one, *AtU, -alpha[0], *V_, *AtU );
00570           MVT::MvNorm( *AtU, xi );
00571           std::cout << "| V alpha - A' u |= "  << xi[0] << std::endl;
00572           // 2. confirm that U is a unit vector
00573           Teuchos::SerialDenseMatrix<int,ScalarType> uotuo(1,1);
00574           MVT::MvTransMv( one, *U_, *U_, uotuo );
00575           std::cout << "<U, U> = " << uotuo << std::endl;
00576           // 3. print alpha =  <V, A'U>
00577           std::cout << "alpha = "  << alpha[0] << std::endl;
00578           // 4. compute < AV, U> which ought to be alpha
00579           Teuchos::SerialDenseMatrix<int,ScalarType> utav(1,1);
00580           MVT::MvTransMv( one, *AV, *U_, utav );
00581           std::cout << "<AV, U> = alpha = " << utav << std::endl;
00582         }
00583 
00584       MVT::MvAddMv( one, *AV, -alpha[0], *U_, *U_ ); // uNew := Av - uOld alphaOld
00585       MVT::MvNorm( *U_, beta);
00586 
00587       anorm2 += alpha[0]*alpha[0] + beta[0]*beta[0] + lambda_*lambda_;
00588 
00589 
00590       if ( SCT::real(beta[0]) > zero )
00591         {
00592       
00593           MVT::MvScale( *U_, one / beta[0] );
00594 
00595           if (M_left.is_null())
00596             { // ... and V_ alphaNew := AtU - V_ betaNew
00597               OPT::Apply(*A, *U_, *AtU, CONJTRANS);
00598             }
00599           else
00600             {
00601               RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
00602               OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS);
00603               OPT::Apply(*A, *tempInRangeOfA, *AtU, CONJTRANS);
00604             }
00605           if (! M_right.is_null())
00606             {
00607               RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
00608               OPT::Apply (*M_right, *tempInDomainOfA, *AtU, CONJTRANS); // AtU may change
00609             }
00610     
00611           MVT::MvAddMv( one, *AtU, -beta[0], *V_, *V_ );
00612           MVT::MvNorm( *V_, alpha );
00613         }
00614       else  // beta = 0
00615         {
00616           alpha[0] = zero;
00617         }
00618 
00619       if ( SCT::real(alpha[0]) > zero )
00620         {
00621           MVT::MvScale( *V_, one / alpha[0] );
00622         }
00623 
00624       // Use a plane rotation to eliminate the damping parameter.  
00625       // This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
00626       common = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_);
00627       cs1 = rhobar / common;
00628       sn1 = lambda_ / common;
00629       psi = sn1 * phibar;
00630       phibar = cs1 * phibar;
00631 
00632       // Use a plane rotation to eliminate the subdiagonal element (beta)
00633       // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
00634       rho = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_ + beta[0]*beta[0]);
00635       cs = common / rho; 
00636       sn = beta[0] / rho;  
00637       theta = sn * alpha[0];
00638       rhobar = -cs * alpha[0];
00639       phi = cs * phibar;
00640       phibar = sn * phibar; // If beta vanishes, so do sn, theta, phibar and eventually resid_norm
00641       tau = sn * phi;
00642 
00643       delta = sn2 * rho;
00644       gammabar = -cs2 * rho;
00645       zetabar = (phi - delta*zeta) / gammabar;
00646       sol_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(xxnorm + zetabar*zetabar);
00647       gamma = Teuchos::ScalarTraits< ScalarType >::squareroot(gammabar*gammabar + theta*theta);
00648       cs2 = gammabar / gamma;
00649       sn2 = theta / gamma;
00650       zeta = (phi - delta*zeta) / gamma;
00651       xxnorm += zeta*zeta;
00652 
00653       //The next task may be addressed by some form of lp_->updateSolution.
00654       if ( M_right.is_null())
00655         {
00656           MVT::MvAddMv( phi / rho, *W_, one, *cur_soln_vec, *cur_soln_vec);
00657         }
00658       else
00659         {
00660           RCP<MV> tempInDomainOfA = MVT::CloneCopy (*W_);
00661           OPT::Apply (*M_right, *W_, *tempInDomainOfA); 
00662           MVT::MvAddMv( phi / rho, *tempInDomainOfA, one, *cur_soln_vec, *cur_soln_vec);
00663         }
00664 
00665       MVT::MvNorm( *W_, wnorm2 );
00666       ddnorm += (one / rho)*(one / rho) * wnorm2[0]*wnorm2[0];
00667       MVT::MvAddMv( one, *V_, -theta / rho, *W_, *W_ );
00668 
00669       frob_mat_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(anorm2);
00670       mat_cond_num_ = frob_mat_norm_ * Teuchos::ScalarTraits< ScalarType >::squareroot(ddnorm);
00671       res+= psi*psi;
00672       resid_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(phibar*phibar + res);
00673       mat_resid_norm_ = alpha[0] * Teuchos::ScalarTraits< ScalarType >::magnitude(tau);
00674 
00675     } // end while (sTest_->checkStatus(this) != Passed)
00676   } // iterate()
00677 
00678 } // end Belos namespace
00679 
00680 #endif /* BELOS_LSQR_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines