BelosGCRODRIter.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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef BELOS_GCRODR_ITER_HPP
00030 #define BELOS_GCRODR_ITER_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 
00039 #include "BelosLinearProblem.hpp"
00040 #include "BelosMatOrthoManager.hpp"
00041 #include "BelosOutputManager.hpp"
00042 #include "BelosStatusTest.hpp"
00043 #include "BelosOperatorTraits.hpp"
00044 #include "BelosMultiVecTraits.hpp"
00045 
00046 #include "Teuchos_BLAS.hpp"
00047 #include "Teuchos_SerialDenseMatrix.hpp"
00048 #include "Teuchos_SerialDenseVector.hpp"
00049 #include "Teuchos_ScalarTraits.hpp"
00050 #include "Teuchos_ParameterList.hpp"
00051 #include "Teuchos_TimeMonitor.hpp"
00052 
00065 namespace Belos {
00066   
00068 
00069   
00074   template <class ScalarType, class MV>
00075   struct GCRODRIterState {
00080     int curDim;
00081     
00083     Teuchos::RCP<MV> V;
00084    
00086     Teuchos::RCP<MV> U, C;
00087 
00093     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H;
00094 
00097     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B;
00098 
00099     GCRODRIterState() : curDim(0), V(Teuchos::null), 
00100       U(Teuchos::null), C(Teuchos::null),
00101       H(Teuchos::null), B(Teuchos::null)
00102     {}
00103   };
00104   
00106   
00108 
00109   
00121   class GCRODRIterInitFailure : public BelosError {
00122     public:
00123       GCRODRIterInitFailure(const std::string& what_arg) : BelosError(what_arg) {}
00124   };
00125   
00132   class GCRODRIterOrthoFailure : public BelosError {
00133     public:
00134       GCRODRIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
00135   };
00136   
00138   
00139   
00140   template<class ScalarType, class MV, class OP>
00141   class GCRODRIter : virtual public Iteration<ScalarType,MV,OP> {
00142     
00143   public:
00144     
00145     //
00146     // Convenience typedefs
00147     //
00148     typedef MultiVecTraits<ScalarType,MV> MVT;
00149     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00150     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00151     typedef typename SCT::magnitudeType MagnitudeType;
00152     
00154 
00155     
00164     GCRODRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00165     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00166     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00167     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00168     Teuchos::ParameterList &params );
00169     
00171     virtual ~GCRODRIter() {};
00173     
00174     
00176 
00177     
00199     void iterate();
00200     
00222     void initialize(GCRODRIterState<ScalarType,MV> newstate);
00223     
00227     void initialize() {
00228       GCRODRIterState<ScalarType,MV> empty;
00229       initialize(empty);
00230     }
00231     
00239     GCRODRIterState<ScalarType,MV> getState() const {
00240       GCRODRIterState<ScalarType,MV> state;
00241       state.curDim = curDim_;
00242       state.V = V_;
00243       state.U = U_;
00244       state.C = C_;
00245       state.H = H_;
00246       state.B = B_;
00247       return state;
00248     }
00249     
00251     
00252     
00254 
00255     
00257     int getNumIters() const { return iter_; }
00258     
00260     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00261     
00264     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00265     
00267 
00272     Teuchos::RCP<MV> getCurrentUpdate() const;
00273     
00275 
00278     void updateLSQR( int dim = -1 );
00279     
00281     int getCurSubspaceDim() const { 
00282       if (!initialized_) return 0;
00283       return curDim_;
00284     };
00285     
00287     int getMaxSubspaceDim() const { return numBlocks_; }
00288     
00290     
00291     
00293 
00294     
00296     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00297     
00299     int getNumBlocks() const { return numBlocks_; }
00300     
00302     void setNumBlocks(int numBlocks) { setSize( recycledBlocks_, numBlocks ); };
00303     
00305     int getRecycledBlocks() const { return recycledBlocks_; }
00306 
00308     void setRecycledBlocks(int recycledBlocks) { setSize( recycledBlocks, numBlocks_ ); };
00309     
00311     int getBlockSize() const { return 1; }
00312     
00314     void setBlockSize(int blockSize) {
00315       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,"Belos::GCRODRIter::setBlockSize(): Cannot use a block size that is not one.");
00316     }
00317 
00319     void setSize( int recycledBlocks, int numBlocks ) {
00320       // only call resize if size changed
00321       if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
00322         recycledBlocks_ = recycledBlocks;
00323         numBlocks_ = numBlocks;
00324         cs_.sizeUninitialized( numBlocks_+1 );
00325         sn_.sizeUninitialized( numBlocks_+1 );
00326         z_.sizeUninitialized( numBlocks_+1 );
00327         R_.shapeUninitialized( numBlocks_+1,numBlocks );
00328       }
00329     }
00330 
00332     bool isInitialized() { return initialized_; }
00333     
00335     
00336   private:
00337     
00338     //
00339     // Internal methods
00340     //
00341     
00342     //
00343     // Classes inputed through constructor that define the linear problem to be solved.
00344     //
00345     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00346     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00347     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00348     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00349     
00350     //
00351     // Algorithmic parameters
00352     //  
00353     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00354     int numBlocks_; 
00355 
00356     // recycledBlocks_ is the size of the allocated space for the recycled subspace, in blocks.
00357     int recycledBlocks_; 
00358     
00359     // Storage for QR factorization of the least squares system.
00360     Teuchos::SerialDenseVector<int,ScalarType> sn_;
00361     Teuchos::SerialDenseVector<int,MagnitudeType> cs_;
00362 
00363     // 
00364     // Current solver state
00365     //
00366     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00367     // is capable of running; _initialize is controlled  by the initialize() member method
00368     // For the implications of the state of initialized_, please see documentation for initialize()
00369     bool initialized_;
00370     
00371     // Current subspace dimension, and number of iterations performed.
00372     int curDim_, iter_;
00373     
00374     // 
00375     // State Storage
00376     //
00377     // Krylov vectors.
00378     Teuchos::RCP<MV> V_;
00379     //
00380     // Recycled subspace vectors.
00381     Teuchos::RCP<MV> U_, C_;
00382     //
00383     // Projected matrices
00384     // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00385     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00386     //
00387     // B_ : Projected matrix from the recycled subspace B = C^H*A*V
00388     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
00389     //
00390     // QR decomposition of Projected matrices for solving the least squares system HY = B.
00391     // R_: Upper triangular reduction of H
00392     // z_: Q applied to right-hand side of the least squares system
00393     Teuchos::SerialDenseMatrix<int,ScalarType> R_;
00394     Teuchos::SerialDenseVector<int,ScalarType> z_;  
00395   };
00396   
00398   // Constructor.
00399   template<class ScalarType, class MV, class OP>
00400   GCRODRIter<ScalarType,MV,OP>::GCRODRIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00401              const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00402              const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00403              const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00404              Teuchos::ParameterList &params ):
00405     lp_(problem), om_(printer), stest_(tester), ortho_(ortho) {
00406     numBlocks_      = 0;
00407     recycledBlocks_ = 0;
00408     initialized_    = false;
00409     curDim_         = 0;
00410     iter_           = 0;
00411     V_              = Teuchos::null;
00412     U_              = Teuchos::null;
00413     C_              = Teuchos::null;
00414     H_              = Teuchos::null;
00415     B_              = Teuchos::null;
00416 
00417     // Get the maximum number of blocks allowed for this Krylov subspace
00418     TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument, "Belos::GCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
00419     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00420 
00421     TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,"Belos::GCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
00422     int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
00423 
00424     TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument, "Belos::GCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
00425     TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument, "Belos::GCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
00426 
00427     numBlocks_ = nb;
00428     recycledBlocks_ = rb;
00429 
00430     cs_.sizeUninitialized( numBlocks_+1 );
00431     sn_.sizeUninitialized( numBlocks_+1 );
00432     z_.sizeUninitialized( numBlocks_+1 );
00433     R_.shapeUninitialized( numBlocks_+1,numBlocks_ );
00434 
00435   }
00436   
00438   // Get the current update from this subspace.
00439   template <class ScalarType, class MV, class OP>
00440   Teuchos::RCP<MV> GCRODRIter<ScalarType,MV,OP>::getCurrentUpdate() const {
00441     //
00442     // If this is the first iteration of the Arnoldi factorization, 
00443     // there is no update, so return Teuchos::null. 
00444     //
00445     RCP<MV> currentUpdate = Teuchos::null;
00446     if (curDim_==0) { 
00447       return currentUpdate; 
00448     } else {
00449       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00450       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00451       Teuchos::BLAS<int,ScalarType> blas;
00452       currentUpdate = MVT::Clone( *V_, 1 );
00453       //
00454       //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00455       //
00456       Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, z_, curDim_, 1 );
00457       //
00458       //  Solve the least squares problem.
00459       //
00460       blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00461      Teuchos::NON_UNIT_DIAG, curDim_, 1, one,  
00462      R_.values(), R_.stride(), y.values(), y.stride() );
00463       //
00464       //  Compute the current update from the Krylov basis; V(:,1:curDim_)*y.
00465       //
00466       std::vector<int> index(curDim_);
00467       for ( int i=0; i<curDim_; i++ ) index[i] = i;
00468       RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
00469       MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
00470       //
00471       //  Add in portion of update from recycled subspace U; U(:,1:recycledBlocks_)*B*y.
00472       //
00473       if (U_ != Teuchos::null) {
00474         Teuchos::SerialDenseMatrix<int,ScalarType> z(recycledBlocks_,1);
00475         Teuchos::SerialDenseMatrix<int,ScalarType> subB( Teuchos::View, *B_, recycledBlocks_, curDim_ );
00476         z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, y, zero );
00477         MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
00478       }
00479     }
00480     return currentUpdate;
00481   }
00482 
00483 
00485   // Get the native residuals stored in this iteration.  
00486   // Note: This method does not return a MultiVector of the residual vectors, only return Teuchos::null
00487   template <class ScalarType, class MV, class OP>
00488   Teuchos::RCP<const MV> GCRODRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const {
00489     //
00490     // NOTE: Make sure the incoming std::vector is the correct size!
00491     //
00492     if ( norms && (int)norms->size()==0 )                         
00493       norms->resize( 1 );                                          
00494     
00495     if (norms) {
00496       Teuchos::BLAS<int,ScalarType> blas;
00497       (*norms)[0] = blas.NRM2( 1, &z_(curDim_), 1);
00498     }
00499     return Teuchos::null;
00500   }
00501   
00502   
00503 
00505   // Initialize this iteration object
00506   template <class ScalarType, class MV, class OP>
00507   void GCRODRIter<ScalarType,MV,OP>::initialize(GCRODRIterState<ScalarType,MV> newstate) {
00508     
00509     if (newstate.V != Teuchos::null &&  newstate.H != Teuchos::null) {
00510       curDim_ = newstate.curDim;
00511       V_      = newstate.V;
00512       U_      = newstate.U;
00513       C_      = newstate.C;
00514       H_      = newstate.H;
00515       B_      = newstate.B;
00516     }
00517     else {
00518       TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): GCRODRIterState does not have V initialized.");
00519       TEST_FOR_EXCEPTION(newstate.H == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): GCRODRIterState does not have H initialized.");
00520     }
00521 
00522     // the solver is initialized
00523     initialized_ = true;
00524 
00525   }
00526 
00527 
00529   // Iterate until the status test informs us we should stop.
00530   template <class ScalarType, class MV, class OP>
00531   void GCRODRIter<ScalarType,MV,OP>::iterate() {
00532 
00533     TEST_FOR_EXCEPTION( initialized_ == false, GCRODRIterInitFailure,"Belos::GCRODRIter::iterate(): GCRODRIter class not initialized." );
00534 
00535     // Force call to setsize to ensure internal storage is correct dimension
00536     setSize( recycledBlocks_, numBlocks_ );
00537     
00538     Teuchos::RCP<MV> Vnext;
00539     Teuchos::RCP<const MV> Vprev;
00540     std::vector<int> curind(1);
00541 
00542     // z_ must be zeroed out in order to compute Givens rotations correctly
00543     z_.putScalar(0.0);
00544 
00545     // Orthonormalize the new V_0
00546     curind[0] = 0;
00547     Vnext = MVT::CloneViewNonConst(*V_,curind);
00548     // Orthonormalize first column
00549     // First, get a monoelemental matrix to hold the orthonormalization coefficients
00550     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z0 = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(1,1) );
00551     int rank = ortho_->normalize( *Vnext, z0 );
00552     TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
00553     // Copy the scalar coefficient back into the z_ vector
00554     z_(0) = (*z0)(0,0);
00555 
00556     std::vector<int> prevind(numBlocks_+1);
00557 
00559     // iterate until the status test tells us to stop.
00560     //
00561     // also break if our basis is full
00562     //
00563     if (U_ == Teuchos::null) { // iterate without recycle space (because none is available yet)
00564       while (stest_->checkStatus(this) != Passed && curDim_+1 <= numBlocks_) {
00565         iter_++;
00566         int lclDim = curDim_ + 1;
00567 
00568         // Get next index in basis
00569         curind[0] = lclDim;
00570         Vnext = MVT::CloneViewNonConst(*V_,curind);
00571 
00572         // Get previous index in basis
00573         curind[0] = curDim_;
00574         Vprev = MVT::CloneView(*V_,curind);
00575 
00576         // Compute the next vector in the Krylov basis:  Vnext = Op*Vprev
00577         lp_->apply(*Vprev,*Vnext);
00578 
00579         // Remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.
00580 
00581         // Get a view of all the previous vectors
00582         prevind.resize(lclDim);
00583         for (int i=0; i<lclDim; i++) { prevind[i] = i; }
00584         Vprev = MVT::CloneView(*V_,prevind);
00585         Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
00586 
00587         // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00588         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00589           subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>  ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
00590         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
00591         AsubH.append( subH );
00592 
00593         // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00594         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00595           subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
00596 
00597         // Project out the previous Krylov vectors and normalize the next vector.
00598         int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
00599 
00600         // Copy over the coefficients to R just in case we run into an error.
00601         Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
00602         Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
00603         subR2.assign(subH2);
00604 
00605         TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
00606 
00607         // Update the QR factorization of the upper Hessenberg matrix
00608         updateLSQR();
00609 
00610         // Update basis dim
00611         curDim_++;
00612       } // end while (statusTest == false)
00613     }
00614     else { // iterate with recycle space
00615       while (stest_->checkStatus(this) != Passed && curDim_+1 <= numBlocks_) {
00616         iter_++;
00617         int lclDim = curDim_ + 1; 
00618 
00619         // Get the current part of the basis.
00620         curind[0] = lclDim;
00621         Vnext = MVT::CloneViewNonConst(*V_,curind);
00622 
00623         // Get a view of the previous vectors.
00624         // This is used for orthogonalization and for computing V^H K H
00625         curind[0] = curDim_;
00626         Vprev = MVT::CloneView(*V_,curind);
00627 
00628         // Compute the next std::vector in the Krylov basis:  Vnext = Op*Vprev
00629         lp_->apply(*Vprev,*Vnext);
00630         Vprev = Teuchos::null;
00631 
00632         // First, remove the recycled subspace (C) from Vnext and put coefficients in B.
00633         Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
00634         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00635       subB = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*B_,recycledBlocks_,1,0,curDim_ ) );
00636         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubB;
00637         AsubB.append( subB );
00638 
00639         // Project out the recycled subspace.
00640         ortho_->project( *Vnext, AsubB, C );
00641 
00642         // Now, remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.          
00643         // Get a view of all the previous vectors
00644         prevind.resize(lclDim);
00645         for (int i=0; i<lclDim; i++) { prevind[i] = i; }
00646         Vprev = MVT::CloneView(*V_,prevind);
00647         Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
00648   
00649         // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00650         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>  ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
00651         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
00652         AsubH.append( subH );
00653       
00654         // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00655         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >  subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
00656 
00657         // Project out the previous Krylov vectors and normalize the next vector.
00658         int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
00659 
00660         // Copy over the coefficients to R just in case we run into an error.
00661         Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
00662         Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
00663         subR2.assign(subH2);
00664       
00665         TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
00666 
00667         // Update the QR factorization of the upper Hessenberg matrix
00668         updateLSQR();
00669 
00670         // Update basis dim
00671         curDim_++;
00672       
00673       } // end while (statusTest == false)
00674     } // end if (U_ == Teuchos::null)
00675    
00676   }
00677 
00678   
00679   template<class ScalarType, class MV, class OP>
00680   void GCRODRIter<ScalarType,MV,OP>::updateLSQR( int dim ) {
00681 
00682     int i;
00683     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00684     
00685     // Get correct dimension based on input "dim"
00686     // Remember that ortho failures result in an exit before updateLSQR() is called.
00687     // Therefore, it is possible that dim == curDim_.
00688     int curDim = curDim_;
00689     if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) )
00690       curDim = dim;
00691     
00692     Teuchos::BLAS<int, ScalarType> blas;
00693     //
00694     // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
00695     // system to upper-triangular form.
00696     //
00697     // QR factorization of Least-Squares system with Givens rotations
00698     //
00699     for (i=0; i<curDim; i++) {
00700       //
00701       // Apply previous Givens rotations to new column of Hessenberg matrix
00702       //
00703       blas.ROT( 1, &R_(i,curDim), 1, &R_(i+1, curDim), 1, &cs_[i], &sn_[i] );
00704     }
00705     //
00706     // Calculate new Givens rotation
00707     //
00708     blas.ROTG( &R_(curDim,curDim), &R_(curDim+1,curDim), &cs_[curDim], &sn_[curDim] );
00709     R_(curDim+1,curDim) = zero;
00710     //
00711     // Update RHS w/ new transformation
00712     //
00713     blas.ROT( 1, &z_(curDim), 1, &z_(curDim+1), 1, &cs_[curDim], &sn_[curDim] );
00714 
00715   } // end updateLSQR()
00716 
00717 } // end Belos namespace
00718 
00719 #endif /* BELOS_GCRODR_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:05 2011 for Belos by  doxygen 1.6.3