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     resid_norm_(0.0),
00301     frob_mat_norm_(0.0),
00302     mat_cond_num_(0.0),
00303     mat_resid_norm_(0.0),
00304     sol_norm_(0.0),
00305     bnorm_(0.0)
00306   {
00307   }
00308 
00310   // Setup the state storage.
00311   template <class ScalarType, class MV, class OP>
00312   void LSQRIter<ScalarType,MV,OP>::setStateSize ()
00313   {
00314     if (!stateStorageInitialized_) {
00315       // Check if there is any multivector to clone from.
00316       Teuchos::RCP<const MV> rhsMV = lp_->getInitPrecResVec();
00317       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00318       if (lhsMV == Teuchos::null || rhsMV == Teuchos::null) {
00319   stateStorageInitialized_ = false;
00320   return;
00321       }
00322       else {
00323   // Initialize the state storage
00324   // If the subspace has not been initialized before, generate it
00325         // using the LHS and RHS from lp_.
00326   if (U_ == Teuchos::null) {
00327     // Get the multivectors.
00328     TEUCHOS_TEST_FOR_EXCEPTION(rhsMV == Teuchos::null, std::invalid_argument, "LSQRIter::setStateSize(): linear problem does not specify right hand multivector to clone from.");
00329     TEUCHOS_TEST_FOR_EXCEPTION(lhsMV == Teuchos::null, std::invalid_argument, "LSQRIter::setStateSize(): linear problem does not specify left hand multivector to clone from.");
00330 
00331     U_ = MVT::Clone( *rhsMV, 1 ); // LeftPrecond * rhs
00332     V_ = MVT::Clone( *lhsMV, 1 ); // zero, overwrittein in
00333     W_ = MVT::Clone( *lhsMV, 1 ); // zero, initializeLSQR
00334   }
00335 
00336   // State storage has now been initialized.
00337   stateStorageInitialized_ = true;
00338       }
00339     }
00340   }
00341 
00342 
00344   // Initialize this iteration object
00345   template <class ScalarType, class MV, class OP>
00346   void LSQRIter<ScalarType,MV,OP>::initializeLSQR(LSQRIterationState<ScalarType,MV> newstate)
00347   {
00348     using Teuchos::RCP;
00349 
00350     // Initialize the state storage if it isn't already.
00351     if (!stateStorageInitialized_)
00352       setStateSize();
00353 
00354     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00355            "LSQRIter::initialize(): Cannot initialize state storage!");
00356 
00357     std::string errstr("LSQRIter::initialize(): Specified multivectors must have a consistent length and width.");
00358 
00359 
00360     // Compute initial bidiagonalization vectors and search direction
00361     //
00362     RCP<const MV> lhsMV = lp_->getLHS(); // contains initial guess,
00363 
00364     const bool debugSerialLSQR = false;
00365 
00366     if( debugSerialLSQR )
00367       {
00368         std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> lhsNorm(1);
00369         MVT::MvNorm( *lhsMV, lhsNorm );
00370         std::cout << "initializeLSQR lhsNorm " << lhsNorm[0] << std::endl;
00371       }
00372 
00373     // LinearProlbem provides right-hand side vectors including RHS CurrRHSVec InitResVec.
00374     RCP<const MV> rhsMV = lp_->getInitPrecResVec();
00375     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00376     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00377     MVT::MvAddMv( one, *rhsMV, zero, *U_, *U_);
00378 
00379     RCP<const OP> M_left = lp_->getLeftPrec();
00380     RCP<const OP> A = lp_->getOperator();
00381     RCP<const OP> M_right = lp_->getRightPrec();
00382 
00383     if( debugSerialLSQR )
00384       {
00385         std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> rhsNorm(1);
00386         MVT::MvNorm( *U_, rhsNorm );
00387         std::cout << "initializeLSQR | U_ | : " << rhsNorm[0] << std::endl;
00388         //U_->Print(std::cout);
00389       }
00390 
00391     //MVT::MvScale( *V_, zero );
00392 
00393     // Apply the (conjugate) transpose of the preconditioned operator:
00394     //
00395     // V := (M_L A M_R)^* U, which means
00396     // V := M_R^* (A^* (M_L^* U)).
00397     //
00398     //OPT::Apply(*(lp_->getOperator()), *U_, *V_, CONJTRANS);
00399     if ( M_left.is_null())
00400       {
00401         OPT::Apply (*A, *U_, *V_, CONJTRANS); // V_ = A' U_
00402         //std::cout << "***************  V_ ****************" << std::endl;
00403         //V_->Print(std::cout);
00404       }
00405     else
00406       {
00407         RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
00408         OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS);
00409         OPT::Apply (*A, *tempInRangeOfA, *V_, CONJTRANS); // V_ = A' LeftPrec' U_
00410         //std::cout << "mLeft  V_ = " << *V_ << std::endl;
00411       }
00412     if (! M_right.is_null())
00413       {
00414         RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
00415 
00416         OPT::Apply (*M_right, *tempInDomainOfA, *V_, CONJTRANS); // V:= RtPrec' A' LeftPrec' U
00417         //std::cout << "mRight  V_ = " << *V_ << std::endl;
00418       }
00419 
00420     // W := V (copy the vector)
00421     MVT::MvAddMv( one, *V_, zero, *V_, *W_);
00422 
00423     frob_mat_norm_ = zero; // These are
00424     mat_cond_num_ = one;   // lower
00425     sol_norm_ = zero;      // bounds.
00426 
00427     // The solver is initialized.
00428     initialized_ = true;
00429   }
00430 
00431 
00433   // Iterate until the status test informs us we should stop.
00434   template <class ScalarType, class MV, class OP>
00435   void LSQRIter<ScalarType,MV,OP>::iterate()
00436   {
00437     //
00438     // Allocate/initialize data structures
00439     //
00440     if (initialized_ == false) {
00441       initialize();
00442     }
00443 
00444     // Create convenience variables for zero and one.
00445     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00446     const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00447     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00448 
00449     // Allocate memory for scalars.
00450     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> alpha(1);
00451     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> beta(1);
00452     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> xi(1);
00453     // xi is a dumb scalar used for storing inner products.
00454     // Eventually SDM will replace the "vectors".
00455     std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> wnorm2(1);
00456     ScalarType rhobar, phibar, cs1, phi, rho, cs, sn, theta, xxnorm = MTzero, common;
00457     ScalarType zetabar, sn1, psi, res = zero, ddnorm = zero, gamma, tau;
00458     ScalarType anorm2 = zero;
00459     ScalarType cs2 = -one, sn2 = zero, gammabar, zeta = zero, delta;
00460 
00461     // The pair of work vectors AV and AtU are
00462     Teuchos::RCP<MV> AV; // used in applying A to V_ and
00463     AV = MVT::Clone( *U_, 1);
00464     Teuchos::RCP<MV> AtU; // used in applying A^TRANS to U_ respectively.
00465     AtU = MVT::Clone( *V_, 1);
00466     const bool debugSerialLSQR = false;
00467 
00468     // Get the current solution vector.
00469     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00470 
00471 
00472     // Check that the current solution vector only has one column.
00473     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, LSQRIterateFailure,
00474                         "LSQRIter::iterate(): current linear system has more than one vector!" );
00475 
00476     // In initializeLSQR among other things V = A' U.
00477     // alpha and beta normalize these vectors.
00478     MVT::MvNorm( *U_, beta );
00479     if( SCT::real(beta[0]) >  zero  )
00480       {
00481         MVT::MvScale( *U_, one / beta[0] );
00482 
00483         //std::cout << "***************  U/beta ****************" << std::endl;
00484         //U_->Print(std::cout);
00485 
00486         MVT::MvScale( *V_, one / beta[0] );  // scale V = A'U to normalize U
00487 
00488         //std::cout << "***************  V/beta ****************" << std::endl;
00489         //V_->Print(std::cout);
00490       }
00491     MVT::MvNorm( *V_, alpha );
00492     if( debugSerialLSQR )
00493       {
00494          // used to compare with implementations
00495          // initializing mat_resid_norm to alpha/beta
00496          std::cout << iter_ << " First alpha " << alpha[0]   << " beta " << beta[0] << " lambda " << lambda_ << std::endl;
00497       }
00498     if( SCT::real(alpha[0]) >  zero  )
00499       {
00500         MVT::MvScale( *V_, one / alpha[0] ); // V alpha = A' U to normalize V
00501       }
00502     if( beta[0] * alpha[0] >  zero  )
00503       {
00504         MVT::MvScale( *W_, one / (beta[0] * alpha[0]) ); // W = V
00505       }
00506     else
00507       {
00508         MVT::MvScale( *W_, zero  );
00509       }
00510 
00511     using Teuchos::RCP;
00512     RCP<const OP> M_left = lp_->getLeftPrec();
00513     RCP<const OP> A = lp_->getOperator();
00514     RCP<const OP> M_right = lp_->getRightPrec();
00515 
00516     rhobar = alpha[0];
00517     phibar = beta[0];
00518     bnorm_ = beta[0];
00519     resid_norm_ = beta[0];
00520     mat_resid_norm_ = alpha[0] * beta[0];
00521 
00522 
00524     // Iterate until the status test tells us to stop.
00525     //
00526     while (stest_->checkStatus(this) != Belos::Passed) {
00527       // Increment the iteration
00528       iter_++;
00529 
00530       // Perform the next step of the bidiagonalization.
00531       // The next U_ and V_ vectors and scalars alpha and beta satisfy
00532       // U_ betaNew := AV - U_ alphaOld ...
00533 
00534       if ( M_right.is_null() )
00535         {
00536           OPT::Apply(*A, *V_, *AV); // AV := A * V_
00537         }
00538       else
00539         {
00540           RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
00541           OPT::Apply (*M_right, *V_, *tempInDomainOfA);
00542           OPT::Apply(*A, *tempInDomainOfA, *AV);
00543         }
00544 
00545       if (! M_left.is_null())
00546         {
00547           RCP<MV> tempInRangeOfA = MVT::CloneCopy (*AV);
00548           OPT::Apply (*M_left, *tempInRangeOfA, *AV); // AV may change
00549         }
00550 
00551 
00552       if ( !( M_left.is_null()  &&  M_right.is_null() )
00553            && debugSerialLSQR && iter_ == 1)
00554         {
00555           // In practice, LSQR may reveal bugs in transposed preconditioners.
00556           // This is the test that catches this type of bug.
00557           // 1. confirm that V alpha = A' U
00558 
00559           if (! M_left.is_null())
00560             {
00561               RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
00562               OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS);
00563               OPT::Apply (*A, *tempInRangeOfA, *AtU, CONJTRANS);   // AtU = B'L'U
00564             }
00565           else
00566             {
00567               OPT::Apply (*A, *U_, *AtU, CONJTRANS);   // AtU = B'U
00568             }
00569           if ( !( M_right.is_null() ) )
00570             {
00571               RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
00572               OPT::Apply (*M_right, *tempInDomainOfA, *AtU, CONJTRANS); // AtU := R' AtU
00573             }
00574 
00575           MVT::MvAddMv( one, *AtU, -alpha[0], *V_, *AtU );
00576           MVT::MvNorm( *AtU, xi );
00577           std::cout << "| V alpha - A' u |= "  << xi[0] << std::endl;
00578           // 2. confirm that U is a unit vector
00579           Teuchos::SerialDenseMatrix<int,ScalarType> uotuo(1,1);
00580           MVT::MvTransMv( one, *U_, *U_, uotuo );
00581           std::cout << "<U, U> = " << uotuo << std::endl;
00582           // 3. print alpha =  <V, A'U>
00583           std::cout << "alpha = "  << alpha[0] << std::endl;
00584           // 4. compute < AV, U> which ought to be alpha
00585           Teuchos::SerialDenseMatrix<int,ScalarType> utav(1,1);
00586           MVT::MvTransMv( one, *AV, *U_, utav );
00587           std::cout << "<AV, U> = alpha = " << utav << std::endl;
00588         }
00589 
00590       MVT::MvAddMv( one, *AV, -alpha[0], *U_, *U_ ); // uNew := Av - uOld alphaOld
00591       MVT::MvNorm( *U_, beta);
00592 
00593       anorm2 += alpha[0]*alpha[0] + beta[0]*beta[0] + lambda_*lambda_;
00594 
00595 
00596       if ( SCT::real(beta[0]) > zero )
00597         {
00598 
00599           MVT::MvScale( *U_, one / beta[0] );
00600 
00601           if (M_left.is_null())
00602             { // ... and V_ alphaNew := AtU - V_ betaNew
00603               OPT::Apply(*A, *U_, *AtU, CONJTRANS);
00604             }
00605           else
00606             {
00607               RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
00608               OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS);
00609               OPT::Apply(*A, *tempInRangeOfA, *AtU, CONJTRANS);
00610             }
00611           if (! M_right.is_null())
00612             {
00613               RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
00614               OPT::Apply (*M_right, *tempInDomainOfA, *AtU, CONJTRANS); // AtU may change
00615             }
00616 
00617           MVT::MvAddMv( one, *AtU, -beta[0], *V_, *V_ );
00618           MVT::MvNorm( *V_, alpha );
00619         }
00620       else  // beta = 0
00621         {
00622           alpha[0] = zero;
00623         }
00624 
00625       if ( SCT::real(alpha[0]) > zero )
00626         {
00627           MVT::MvScale( *V_, one / alpha[0] );
00628         }
00629 
00630       // Use a plane rotation to eliminate the damping parameter.
00631       // This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
00632       common = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_);
00633       cs1 = rhobar / common;
00634       sn1 = lambda_ / common;
00635       psi = sn1 * phibar;
00636       phibar = cs1 * phibar;
00637 
00638       // Use a plane rotation to eliminate the subdiagonal element (beta)
00639       // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
00640       rho = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_ + beta[0]*beta[0]);
00641       cs = common / rho;
00642       sn = beta[0] / rho;
00643       theta = sn * alpha[0];
00644       rhobar = -cs * alpha[0];
00645       phi = cs * phibar;
00646       phibar = sn * phibar; // If beta vanishes, so do sn, theta, phibar and eventually resid_norm
00647       tau = sn * phi;
00648 
00649       delta = sn2 * rho;
00650       gammabar = -cs2 * rho;
00651       zetabar = (phi - delta*zeta) / gammabar;
00652       sol_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(xxnorm + zetabar*zetabar);
00653       gamma = Teuchos::ScalarTraits< ScalarType >::squareroot(gammabar*gammabar + theta*theta);
00654       cs2 = gammabar / gamma;
00655       sn2 = theta / gamma;
00656       zeta = (phi - delta*zeta) / gamma;
00657       xxnorm += zeta*zeta;
00658 
00659       //The next task may be addressed by some form of lp_->updateSolution.
00660       if ( M_right.is_null())
00661         {
00662           MVT::MvAddMv( phi / rho, *W_, one, *cur_soln_vec, *cur_soln_vec);
00663         }
00664       else
00665         {
00666           RCP<MV> tempInDomainOfA = MVT::CloneCopy (*W_);
00667           OPT::Apply (*M_right, *W_, *tempInDomainOfA);
00668           MVT::MvAddMv( phi / rho, *tempInDomainOfA, one, *cur_soln_vec, *cur_soln_vec);
00669         }
00670 
00671       MVT::MvNorm( *W_, wnorm2 );
00672       ddnorm += (one / rho)*(one / rho) * wnorm2[0]*wnorm2[0];
00673       MVT::MvAddMv( one, *V_, -theta / rho, *W_, *W_ );
00674 
00675       frob_mat_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(anorm2);
00676       mat_cond_num_ = frob_mat_norm_ * Teuchos::ScalarTraits< ScalarType >::squareroot(ddnorm);
00677       res+= psi*psi;
00678       resid_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(phibar*phibar + res);
00679       mat_resid_norm_ = alpha[0] * Teuchos::ScalarTraits< ScalarType >::magnitude(tau);
00680 
00681     } // end while (sTest_->checkStatus(this) != Passed)
00682   } // iterate()
00683 
00684 } // end Belos namespace
00685 
00686 #endif /* BELOS_LSQR_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines