Belos Version of the Day
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 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_GCRODR_ITER_HPP
00043 #define BELOS_GCRODR_ITER_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 
00052 #include "BelosLinearProblem.hpp"
00053 #include "BelosMatOrthoManager.hpp"
00054 #include "BelosOutputManager.hpp"
00055 #include "BelosStatusTest.hpp"
00056 #include "BelosOperatorTraits.hpp"
00057 #include "BelosMultiVecTraits.hpp"
00058 
00059 #include "Teuchos_BLAS.hpp"
00060 #include "Teuchos_SerialDenseMatrix.hpp"
00061 #include "Teuchos_SerialDenseVector.hpp"
00062 #include "Teuchos_ScalarTraits.hpp"
00063 #include "Teuchos_ParameterList.hpp"
00064 #include "Teuchos_TimeMonitor.hpp"
00065 
00078 namespace Belos {
00079   
00081 
00082   
00087   template <class ScalarType, class MV>
00088   struct GCRODRIterState {
00093     int curDim;
00094     
00096     Teuchos::RCP<MV> V;
00097    
00099     Teuchos::RCP<MV> U, C;
00100 
00106     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H;
00107 
00110     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B;
00111 
00112     GCRODRIterState() : curDim(0), V(Teuchos::null), 
00113       U(Teuchos::null), C(Teuchos::null),
00114       H(Teuchos::null), B(Teuchos::null)
00115     {}
00116   };
00117   
00119   
00121 
00122   
00134   class GCRODRIterInitFailure : public BelosError {
00135     public:
00136       GCRODRIterInitFailure(const std::string& what_arg) : BelosError(what_arg) {}
00137   };
00138   
00145   class GCRODRIterOrthoFailure : public BelosError {
00146     public:
00147       GCRODRIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
00148   };
00149   
00151   
00152   
00153   template<class ScalarType, class MV, class OP>
00154   class GCRODRIter : virtual public Iteration<ScalarType,MV,OP> {
00155     
00156   public:
00157     
00158     //
00159     // Convenience typedefs
00160     //
00161     typedef MultiVecTraits<ScalarType,MV> MVT;
00162     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00163     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00164     typedef typename SCT::magnitudeType MagnitudeType;
00165     
00167 
00168     
00177     GCRODRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00178     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00179     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00180     const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00181     Teuchos::ParameterList &params );
00182     
00184     virtual ~GCRODRIter() {};
00186     
00187     
00189 
00190     
00212     void iterate();
00213     
00235     void initialize(GCRODRIterState<ScalarType,MV> newstate);
00236     
00240     void initialize() {
00241       GCRODRIterState<ScalarType,MV> empty;
00242       initialize(empty);
00243     }
00244     
00252     GCRODRIterState<ScalarType,MV> getState() const {
00253       GCRODRIterState<ScalarType,MV> state;
00254       state.curDim = curDim_;
00255       state.V = V_;
00256       state.U = U_;
00257       state.C = C_;
00258       state.H = H_;
00259       state.B = B_;
00260       return state;
00261     }
00262     
00264     
00265     
00267 
00268     
00270     int getNumIters() const { return iter_; }
00271     
00273     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00274     
00277     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
00278     
00280 
00285     Teuchos::RCP<MV> getCurrentUpdate() const;
00286     
00288 
00291     void updateLSQR( int dim = -1 );
00292     
00294     int getCurSubspaceDim() const { 
00295       if (!initialized_) return 0;
00296       return curDim_;
00297     };
00298     
00300     int getMaxSubspaceDim() const { return numBlocks_; }
00301     
00303     
00304     
00306 
00307     
00309     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00310     
00312     int getNumBlocks() const { return numBlocks_; }
00313     
00315     void setNumBlocks(int numBlocks) { setSize( recycledBlocks_, numBlocks ); };
00316     
00318     int getRecycledBlocks() const { return recycledBlocks_; }
00319 
00321     void setRecycledBlocks(int recycledBlocks) { setSize( recycledBlocks, numBlocks_ ); };
00322     
00324     int getBlockSize() const { return 1; }
00325     
00327     void setBlockSize(int blockSize) {
00328       TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,"Belos::GCRODRIter::setBlockSize(): Cannot use a block size that is not one.");
00329     }
00330 
00332     void setSize( int recycledBlocks, int numBlocks ) {
00333       // only call resize if size changed
00334       if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
00335         recycledBlocks_ = recycledBlocks;
00336         numBlocks_ = numBlocks;
00337         cs_.sizeUninitialized( numBlocks_+1 );
00338         sn_.sizeUninitialized( numBlocks_+1 );
00339         z_.sizeUninitialized( numBlocks_+1 );
00340         R_.shapeUninitialized( numBlocks_+1,numBlocks );
00341       }
00342     }
00343 
00345     bool isInitialized() { return initialized_; }
00346     
00348     
00349   private:
00350     
00351     //
00352     // Internal methods
00353     //
00354     
00355     //
00356     // Classes inputed through constructor that define the linear problem to be solved.
00357     //
00358     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00359     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00360     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00361     const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00362     
00363     //
00364     // Algorithmic parameters
00365     //  
00366     // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
00367     int numBlocks_; 
00368 
00369     // recycledBlocks_ is the size of the allocated space for the recycled subspace, in blocks.
00370     int recycledBlocks_; 
00371     
00372     // Storage for QR factorization of the least squares system.
00373     Teuchos::SerialDenseVector<int,ScalarType> sn_;
00374     Teuchos::SerialDenseVector<int,MagnitudeType> cs_;
00375 
00376     // 
00377     // Current solver state
00378     //
00379     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00380     // is capable of running; _initialize is controlled  by the initialize() member method
00381     // For the implications of the state of initialized_, please see documentation for initialize()
00382     bool initialized_;
00383     
00384     // Current subspace dimension, and number of iterations performed.
00385     int curDim_, iter_;
00386     
00387     // 
00388     // State Storage
00389     //
00390     // Krylov vectors.
00391     Teuchos::RCP<MV> V_;
00392     //
00393     // Recycled subspace vectors.
00394     Teuchos::RCP<MV> U_, C_;
00395     //
00396     // Projected matrices
00397     // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
00398     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
00399     //
00400     // B_ : Projected matrix from the recycled subspace B = C^H*A*V
00401     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
00402     //
00403     // QR decomposition of Projected matrices for solving the least squares system HY = B.
00404     // R_: Upper triangular reduction of H
00405     // z_: Q applied to right-hand side of the least squares system
00406     Teuchos::SerialDenseMatrix<int,ScalarType> R_;
00407     Teuchos::SerialDenseVector<int,ScalarType> z_;  
00408   };
00409   
00411   // Constructor.
00412   template<class ScalarType, class MV, class OP>
00413   GCRODRIter<ScalarType,MV,OP>::GCRODRIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00414              const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00415              const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00416              const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00417              Teuchos::ParameterList &params ):
00418     lp_(problem), om_(printer), stest_(tester), ortho_(ortho) {
00419     numBlocks_      = 0;
00420     recycledBlocks_ = 0;
00421     initialized_    = false;
00422     curDim_         = 0;
00423     iter_           = 0;
00424     V_              = Teuchos::null;
00425     U_              = Teuchos::null;
00426     C_              = Teuchos::null;
00427     H_              = Teuchos::null;
00428     B_              = Teuchos::null;
00429 
00430     // Get the maximum number of blocks allowed for this Krylov subspace
00431     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument, "Belos::GCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
00432     int nb = Teuchos::getParameter<int>(params, "Num Blocks");
00433 
00434     TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,"Belos::GCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
00435     int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
00436 
00437     TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument, "Belos::GCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
00438     TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument, "Belos::GCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
00439 
00440     numBlocks_ = nb;
00441     recycledBlocks_ = rb;
00442 
00443     cs_.sizeUninitialized( numBlocks_+1 );
00444     sn_.sizeUninitialized( numBlocks_+1 );
00445     z_.sizeUninitialized( numBlocks_+1 );
00446     R_.shapeUninitialized( numBlocks_+1,numBlocks_ );
00447 
00448   }
00449   
00451   // Get the current update from this subspace.
00452   template <class ScalarType, class MV, class OP>
00453   Teuchos::RCP<MV> GCRODRIter<ScalarType,MV,OP>::getCurrentUpdate() const {
00454     //
00455     // If this is the first iteration of the Arnoldi factorization, 
00456     // there is no update, so return Teuchos::null. 
00457     //
00458     Teuchos::RCP<MV> currentUpdate = Teuchos::null;
00459     if (curDim_==0) { 
00460       return currentUpdate; 
00461     } else {
00462       const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00463       const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00464       Teuchos::BLAS<int,ScalarType> blas;
00465       currentUpdate = MVT::Clone( *V_, 1 );
00466       //
00467       //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
00468       //
00469       Teuchos::SerialDenseMatrix<int,ScalarType> y( Teuchos::Copy, z_, curDim_, 1 );
00470       //
00471       //  Solve the least squares problem.
00472       //
00473       blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
00474      Teuchos::NON_UNIT_DIAG, curDim_, 1, one,  
00475      R_.values(), R_.stride(), y.values(), y.stride() );
00476       //
00477       //  Compute the current update from the Krylov basis; V(:,1:curDim_)*y.
00478       //
00479       std::vector<int> index(curDim_);
00480       for ( int i=0; i<curDim_; i++ ) index[i] = i;
00481       Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
00482       MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *currentUpdate );
00483       //
00484       //  Add in portion of update from recycled subspace U; U(:,1:recycledBlocks_)*B*y.
00485       //
00486       if (U_ != Teuchos::null) {
00487         Teuchos::SerialDenseMatrix<int,ScalarType> z(recycledBlocks_,1);
00488         Teuchos::SerialDenseMatrix<int,ScalarType> subB( Teuchos::View, *B_, recycledBlocks_, curDim_ );
00489         z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, y, zero );
00490         MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
00491       }
00492     }
00493     return currentUpdate;
00494   }
00495 
00496 
00498   // Get the native residuals stored in this iteration.  
00499   // Note: This method does not return a MultiVector of the residual vectors, only return Teuchos::null
00500   template <class ScalarType, class MV, class OP>
00501   Teuchos::RCP<const MV> GCRODRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const {
00502     //
00503     // NOTE: Make sure the incoming std::vector is the correct size!
00504     //
00505     if ( norms && (int)norms->size()==0 )                         
00506       norms->resize( 1 );                                          
00507     
00508     if (norms) {
00509       Teuchos::BLAS<int,ScalarType> blas;
00510       (*norms)[0] = blas.NRM2( 1, &z_(curDim_), 1);
00511     }
00512     return Teuchos::null;
00513   }
00514   
00515   
00516 
00518   // Initialize this iteration object
00519   template <class ScalarType, class MV, class OP>
00520   void GCRODRIter<ScalarType,MV,OP>::initialize(GCRODRIterState<ScalarType,MV> newstate) {
00521     
00522     if (newstate.V != Teuchos::null &&  newstate.H != Teuchos::null) {
00523       curDim_ = newstate.curDim;
00524       V_      = newstate.V;
00525       U_      = newstate.U;
00526       C_      = newstate.C;
00527       H_      = newstate.H;
00528       B_      = newstate.B;
00529     }
00530     else {
00531       TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): GCRODRIterState does not have V initialized.");
00532       TEUCHOS_TEST_FOR_EXCEPTION(newstate.H == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): GCRODRIterState does not have H initialized.");
00533     }
00534 
00535     // the solver is initialized
00536     initialized_ = true;
00537 
00538   }
00539 
00540 
00542   // Iterate until the status test informs us we should stop.
00543   template <class ScalarType, class MV, class OP>
00544   void GCRODRIter<ScalarType,MV,OP>::iterate() {
00545 
00546     TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, GCRODRIterInitFailure,"Belos::GCRODRIter::iterate(): GCRODRIter class not initialized." );
00547 
00548     // Force call to setsize to ensure internal storage is correct dimension
00549     setSize( recycledBlocks_, numBlocks_ );
00550     
00551     Teuchos::RCP<MV> Vnext;
00552     Teuchos::RCP<const MV> Vprev;
00553     std::vector<int> curind(1);
00554 
00555     // z_ must be zeroed out in order to compute Givens rotations correctly
00556     z_.putScalar(0.0);
00557 
00558     // Orthonormalize the new V_0
00559     curind[0] = 0;
00560     Vnext = MVT::CloneViewNonConst(*V_,curind);
00561     // Orthonormalize first column
00562     // First, get a monoelemental matrix to hold the orthonormalization coefficients
00563     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > z0 =
00564       Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(1,1) );
00565     int rank = ortho_->normalize( *Vnext, z0 );
00566     TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
00567     // Copy the scalar coefficient back into the z_ vector
00568     z_(0) = (*z0)(0,0);
00569 
00570     std::vector<int> prevind(numBlocks_+1);
00571 
00573     // iterate until the status test tells us to stop.
00574     //
00575     // also break if our basis is full
00576     //
00577     if (U_ == Teuchos::null) { // iterate without recycle space (because none is available yet)
00578       while (stest_->checkStatus(this) != Passed && curDim_+1 <= numBlocks_) {
00579         iter_++;
00580         int lclDim = curDim_ + 1;
00581 
00582         // Get next index in basis
00583         curind[0] = lclDim;
00584         Vnext = MVT::CloneViewNonConst(*V_,curind);
00585 
00586         // Get previous index in basis
00587         curind[0] = curDim_;
00588         Vprev = MVT::CloneView(*V_,curind);
00589 
00590         // Compute the next vector in the Krylov basis:  Vnext = Op*Vprev
00591         lp_->apply(*Vprev,*Vnext);
00592 
00593         // Remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.
00594 
00595         // Get a view of all the previous vectors
00596         prevind.resize(lclDim);
00597         for (int i=0; i<lclDim; i++) { prevind[i] = i; }
00598         Vprev = MVT::CloneView(*V_,prevind);
00599         Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
00600 
00601         // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00602         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00603           subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>  ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
00604         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
00605         AsubH.append( subH );
00606 
00607         // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00608         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00609           subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
00610 
00611         // Project out the previous Krylov vectors and normalize the next vector.
00612         int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
00613 
00614         // Copy over the coefficients to R just in case we run into an error.
00615         Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
00616         Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
00617         subR2.assign(subH2);
00618 
00619         TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
00620 
00621         // Update the QR factorization of the upper Hessenberg matrix
00622         updateLSQR();
00623 
00624         // Update basis dim
00625         curDim_++;
00626       } // end while (statusTest == false)
00627     }
00628     else { // iterate with recycle space
00629       while (stest_->checkStatus(this) != Passed && curDim_+1 <= numBlocks_) {
00630         iter_++;
00631         int lclDim = curDim_ + 1; 
00632 
00633         // Get the current part of the basis.
00634         curind[0] = lclDim;
00635         Vnext = MVT::CloneViewNonConst(*V_,curind);
00636 
00637         // Get a view of the previous vectors.
00638         // This is used for orthogonalization and for computing V^H K H
00639         curind[0] = curDim_;
00640         Vprev = MVT::CloneView(*V_,curind);
00641 
00642         // Compute the next std::vector in the Krylov basis:  Vnext = Op*Vprev
00643         lp_->apply(*Vprev,*Vnext);
00644         Vprev = Teuchos::null;
00645 
00646         // First, remove the recycled subspace (C) from Vnext and put coefficients in B.
00647         Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
00648         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >
00649       subB = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*B_,recycledBlocks_,1,0,curDim_ ) );
00650         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubB;
00651         AsubB.append( subB );
00652 
00653         // Project out the recycled subspace.
00654         ortho_->project( *Vnext, AsubB, C );
00655 
00656         // Now, remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.          
00657         // Get a view of all the previous vectors
00658         prevind.resize(lclDim);
00659         for (int i=0; i<lclDim; i++) { prevind[i] = i; }
00660         Vprev = MVT::CloneView(*V_,prevind);
00661         Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
00662   
00663         // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
00664         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > subH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>  ( Teuchos::View,*H_,lclDim,1,0,curDim_ ) );
00665         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > AsubH;
00666         AsubH.append( subH );
00667       
00668         // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
00669         Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> >  subR = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View,*H_,1,1,lclDim,curDim_ ) );
00670 
00671         // Project out the previous Krylov vectors and normalize the next vector.
00672         int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
00673 
00674         // Copy over the coefficients to R just in case we run into an error.
00675         Teuchos::SerialDenseMatrix<int,ScalarType> subR2( Teuchos::View,R_,lclDim+1,1,0,curDim_ );
00676         Teuchos::SerialDenseMatrix<int,ScalarType> subH2( Teuchos::View,*H_,lclDim+1,1,0,curDim_ );
00677         subR2.assign(subH2);
00678       
00679         TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
00680 
00681         // Update the QR factorization of the upper Hessenberg matrix
00682         updateLSQR();
00683 
00684         // Update basis dim
00685         curDim_++;
00686       
00687       } // end while (statusTest == false)
00688     } // end if (U_ == Teuchos::null)
00689    
00690   }
00691 
00692   
00693   template<class ScalarType, class MV, class OP>
00694   void GCRODRIter<ScalarType,MV,OP>::updateLSQR( int dim ) {
00695 
00696     int i;
00697     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00698     
00699     // Get correct dimension based on input "dim"
00700     // Remember that ortho failures result in an exit before updateLSQR() is called.
00701     // Therefore, it is possible that dim == curDim_.
00702     int curDim = curDim_;
00703     if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) )
00704       curDim = dim;
00705     
00706     Teuchos::BLAS<int, ScalarType> blas;
00707     //
00708     // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
00709     // system to upper-triangular form.
00710     //
00711     // QR factorization of Least-Squares system with Givens rotations
00712     //
00713     for (i=0; i<curDim; i++) {
00714       //
00715       // Apply previous Givens rotations to new column of Hessenberg matrix
00716       //
00717       blas.ROT( 1, &R_(i,curDim), 1, &R_(i+1, curDim), 1, &cs_[i], &sn_[i] );
00718     }
00719     //
00720     // Calculate new Givens rotation
00721     //
00722     blas.ROTG( &R_(curDim,curDim), &R_(curDim+1,curDim), &cs_[curDim], &sn_[curDim] );
00723     R_(curDim+1,curDim) = zero;
00724     //
00725     // Update RHS w/ new transformation
00726     //
00727     blas.ROT( 1, &z_(curDim), 1, &z_(curDim+1), 1, &cs_[curDim], &sn_[curDim] );
00728 
00729   } // end updateLSQR()
00730 
00731 } // end Belos namespace
00732 
00733 #endif /* BELOS_GCRODR_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines