BelosBlockCG.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 // This file contains an implementation of the Block CG algorithm
00030 // for solving real symmetric positive definite linear systems of 
00031 // equations AX = B, where B is a matrix containing one or more right-hand
00032 // sides, and X is the matrix of corresponding solutions. This implementation 
00033 // allows the user to solve systems involving any number of right-hand sides.
00034 // The block size used in the solver is user specified, and is independent
00035 // of the number of right-hand sides. Thus, a system involving many 
00036 // right-hand sides can be processed by solving for only some number  
00037 // (the block size) of the right-hand sides simultaneously. Several passes
00038 // through the block solver are used to solve for all of them. A single
00039 // right-hand side system can be solved in the traditional way by choosing
00040 // the block size equal to one, or it can be solved using a block 
00041 // implementation (choosing a block size greater than one).
00042 //
00043 //
00044 //
00045 #ifndef BELOS_BLOCK_CG_HPP
00046 #define BELOS_BLOCK_CG_HPP
00047 
00055 #include "BelosConfigDefs.hpp"
00056 #include "BelosIterativeSolver.hpp"
00057 #include "BelosLinearProblem.hpp"
00058 #include "BelosOutputManager.hpp"
00059 #include "BelosOperator.hpp"
00060 #include "BelosStatusTest.hpp"
00061 #include "BelosMultiVecTraits.hpp"
00062 #include "BelosCG.hpp"
00063 
00064 #include "Teuchos_LAPACK.hpp"
00065 #include "Teuchos_SerialDenseMatrix.hpp"
00066 #include "Teuchos_SerialDenseVector.hpp"
00067 #include "Teuchos_ScalarTraits.hpp"
00068 #include "Teuchos_TimeMonitor.hpp"
00069 
00079 namespace Belos {
00080 
00081 template <class ScalarType, class MV, class OP>
00082 class BlockCG : public IterativeSolver<ScalarType,MV,OP> { 
00083 public:
00084   //
00085   // Convenience typedefs
00086   //
00087   typedef MultiVecTraits<ScalarType,MV> MVT;
00088   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00089   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00090   typedef typename SCT::magnitudeType MagnitudeType;
00091 
00093 
00094 
00096   BlockCG(const RefCountPtr<LinearProblem<ScalarType,MV,OP> > &lp,
00097     const RefCountPtr<StatusTest<ScalarType,MV,OP> > &stest,
00098     const RefCountPtr<OutputManager<ScalarType> > &om,    
00099     const RefCountPtr<ParameterList> &pl = Teuchos::null
00100     );
00101 
00103   virtual ~BlockCG() {};
00105   
00107 
00108   
00110   int GetNumIters() const { return( _iter ); }
00111   
00113   int GetNumRestarts() const { return(0); }
00114    
00116 
00119   RefCountPtr<const MV> GetNativeResiduals( std::vector<MagnitudeType> *normvec ) const;
00120   
00122 
00127   RefCountPtr<MV> GetCurrentSoln() { return MVT::CloneCopy( *_cur_block_sol ); };
00128   
00130 
00132   RefCountPtr<LinearProblem<ScalarType,MV,OP> > GetLinearProblem() const { return( _lp ); }
00133 
00134   RefCountPtr<StatusTest<ScalarType,MV,OP> > GetStatusTest() const { return( _stest ); }
00135 
00137   
00139 
00140   
00145   void Solve();
00147  
00149 
00150 
00152   std::string description() const;
00153 
00155 
00156 private:
00157 
00158   void SetCGBlkTols();
00159   bool QRFactorDef(MV&, Teuchos::SerialDenseMatrix<int,ScalarType>&, std::vector<int>* );
00160   void CheckOrthogonality(MV&, MV&);
00161   void CheckAOrthogonality(MV&, MV&);
00162   void PrintCGIterInfo(const std::vector<int> &ind);
00163   void BlockIteration();
00164 
00166   RefCountPtr<LinearProblem<ScalarType,MV,OP> > _lp; 
00167 
00169   RefCountPtr<StatusTest<ScalarType,MV,OP> > _stest; 
00170 
00172   RefCountPtr<OutputManager<ScalarType> > _om;
00173 
00175   RefCountPtr<ParameterList> _pl;     
00176   
00178   RefCountPtr<MV> _cur_block_sol;
00179 
00181   RefCountPtr<MV> _cur_block_rhs; 
00182 
00184   RefCountPtr<MV> _residvecs;
00185 
00187   RefCountPtr<ostream> _os;
00188 
00190   int _blocksize, _iter, _new_blk;
00191 
00193   MagnitudeType _prec, _dep_tol;
00194   
00196   bool _restartTimers;
00197 
00199   Teuchos::RefCountPtr<Teuchos::Time> _timerOp, _timerPrec, 
00200                                       _timerOrtho, _timerTotal;
00201 };
00202   
00203   //
00204   // Implementation
00205   //
00206   
00207   template <class ScalarType, class MV, class OP>
00208   BlockCG<ScalarType,MV,OP>::BlockCG(const RefCountPtr<LinearProblem<ScalarType,MV,OP> > &lp,
00209              const RefCountPtr<StatusTest<ScalarType,MV,OP> > &stest,
00210              const RefCountPtr<OutputManager<ScalarType> > &om,
00211              const RefCountPtr<ParameterList> &pl
00212              ) : 
00213     _lp(lp), 
00214     _stest(stest),
00215     _om(om),
00216     _pl(pl),
00217     _os(om->GetOStream()),
00218     _blocksize(0), 
00219     _iter(0),
00220     _new_blk(1),
00221     _prec(1.0), 
00222     _dep_tol(1.0),
00223     _restartTimers(true),
00224     _timerOp(Teuchos::TimeMonitor::getNewTimer("Operation Op*x")),
00225     _timerPrec(Teuchos::TimeMonitor::getNewTimer("Operation Prec*x")),
00226     _timerOrtho(Teuchos::TimeMonitor::getNewTimer("Orthogonalization")),
00227     _timerTotal(Teuchos::TimeMonitor::getNewTimer("Total time"))
00228   { 
00229     //
00230     // Set the block orthogonality tolerances
00231     //
00232     SetCGBlkTols();
00233   }
00234   
00235   template <class ScalarType, class MV, class OP>
00236   void 
00237   BlockCG<ScalarType,MV,OP>::SetCGBlkTols() 
00238   {
00239     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00240     const MagnitudeType two = 2.0;
00241     const MagnitudeType eps = SCT::eps();
00242     _prec = eps;
00243     _dep_tol = MGT::one()/MGT::squareroot(two);
00244   }
00245   
00246   template <class ScalarType, class MV, class OP>
00247   RefCountPtr<const MV> 
00248   BlockCG<ScalarType,MV,OP>::GetNativeResiduals( std::vector<MagnitudeType> *normvec ) const 
00249   {
00250     int i;
00251     std::vector<int> index( _blocksize );
00252     if (_new_blk)
00253       for (i=0; i<_blocksize; i++) { index[i] = i; }
00254     else
00255       for (i=0; i<_blocksize; i++) { index[i] = _blocksize + i; }
00256     RefCountPtr<MV> ResidMV = MVT::CloneView( *_residvecs, index );
00257     return ResidMV;
00258   }
00259   
00260 
00261   template <class ScalarType, class MV, class OP>
00262   void 
00263   BlockCG<ScalarType,MV,OP>::Solve () 
00264   {
00265     Teuchos::TimeMonitor LocalTimer(*_timerTotal,_restartTimers);
00266 
00267     if ( _restartTimers ) {
00268       _timerOp->reset();
00269       _timerPrec->reset();
00270       _timerOrtho->reset();
00271     }
00272     //
00273     // Get the current ostream from the OutputManager
00274     //
00275     _os = _om->GetOStream();
00276     //
00277     // Retrieve the first linear system to be solved.
00278     //
00279     _cur_block_sol = _lp->GetCurrLHSVec();
00280     _cur_block_rhs = _lp->GetCurrRHSVec();
00281     //
00282     //  Start executable statements. 
00283     //
00284     while (_cur_block_sol.get() && _cur_block_rhs.get() ) {
00285       //
00286       // Get the blocksize for this set of linear systems.
00287       //
00288       _blocksize = _lp->GetCurrBlockSize();
00289       if (_blocksize == 1 ) {
00290   //
00291   // Create single vector CG solver for this linear system.
00292   //
00293   Belos::CG<ScalarType,MV,OP> CGSolver(_lp, _stest, _om);
00294   CGSolver.Solve();
00295   //
00296       } else {
00297   BlockIteration();
00298       }
00299       //
00300       // Get the next block of linear systems, if it returns the null pointer we are done.
00301       //
00302       _cur_block_sol = _lp->GetCurrLHSVec();
00303       _cur_block_rhs = _lp->GetCurrRHSVec();
00304       //
00305     } // end while ( _cur_block_sol && _cur_block_rhs )
00306     // **********************************************************************************
00307     //
00308     // Print timing details 
00309 
00310     // Stop timer.
00311     _timerTotal->stop();
00312 
00313     // Reset format that will be used to print the summary
00314     Teuchos::TimeMonitor::format().setPageWidth(54);
00315 
00316     if (_om->isVerbosity( Belos::TimingDetails )) {
00317       if (_om->doPrint())
00318         *_os <<"********************TIMING DETAILS********************"<<endl;
00319       Teuchos::TimeMonitor::summarize( *_os );
00320       if (_om->doPrint())
00321         *_os <<"******************************************************"<<endl;
00322     }
00323     //
00324   } // end Solve()
00325   //  
00326   
00327   template <class ScalarType, class MV, class OP>
00328   void 
00329   BlockCG<ScalarType,MV,OP>::BlockIteration ( ) 
00330   {
00331     //
00332     int i, j, k, info, num_ind;
00333     int ind_blksz, prev_ind_blksz;
00334     bool exit_flg = false;
00335     bool isPrec = (_lp->GetLeftPrec().get()!=NULL);
00336     char UPLO = 'U';
00337     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00338     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00339     Teuchos::LAPACK<int,ScalarType> lapack;
00340     RefCountPtr<MV> P_prev, R_prev, AP_prev, P_new, R_new;
00341     //
00342     // Make additional space needed during iteration
00343     //
00344     std::vector<int> index2( _blocksize );
00345     std::vector<int> index( _blocksize );
00346     std::vector<int> ind_idx( _blocksize );
00347     std::vector<int> cols( _blocksize );
00348     //
00349     // Make room for the direction, residual, and operator applied to direction vectors.
00350     // We save 3 blocks of these vectors.  Also create 2 vectors of _blocksize to store
00351     // The residual norms used to determine valid direction/residual vectors.
00352     //
00353     RefCountPtr<MV> _basisvecs = MVT::Clone( *_cur_block_sol, 2*_blocksize ); 
00354     RefCountPtr<MV> temp_blk = MVT::Clone( *_cur_block_sol, _blocksize );
00355     std::vector<MagnitudeType> _cur_resid_norms( _blocksize );
00356     std::vector<MagnitudeType> _init_resid_norms( _blocksize );
00357     _residvecs = MVT::Clone( *_cur_block_sol, 2*_blocksize ); 
00358     //
00359     ind_blksz = _blocksize;
00360     prev_ind_blksz = _blocksize;
00361     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( _blocksize, _blocksize );
00362     Teuchos::SerialDenseMatrix<int,ScalarType> beta( _blocksize, _blocksize );
00363     Teuchos::SerialDenseMatrix<int,ScalarType> T2( _blocksize, _blocksize );
00364     //
00365     for (i=0;i<_blocksize;i++){
00366       index[i] = i;
00367       ind_idx[i] = i; 
00368       cols[i] = i;
00369     }
00370     //
00371     if (_om->isVerbosityAndPrint( IterationDetails )) {
00372       *_os << endl;
00373       *_os << "===================================================" << endl;
00374       *_os << "Solving linear system(s):  " << _lp->GetRHSIndex() 
00375      << " through " << _lp->GetRHSIndex()+_lp->GetNumToSolve() << endl;
00376       *_os << endl;
00377     } 
00378     //
00379     //
00380     // ************ Compute the initial residuals ********************************
00381     //
00382     // Associate the first block of _basisvecs with P_prev and the
00383     // first block of _residvecs with R_prev
00384     //
00385     P_prev = MVT::CloneView( *_basisvecs, ind_idx );
00386     R_prev = MVT::CloneView( *_residvecs, ind_idx );
00387     AP_prev = MVT::CloneView( *temp_blk, ind_idx );
00388     //
00389     // Store initial guesses to AX = B in 1st block of _basisvecs
00390     //         P_prev = one*cur_block_sol + zero*P_prev
00391     //
00392     MVT::MvAddMv(one, *_cur_block_sol, zero, *_cur_block_sol, *P_prev);
00393     //
00394     // Multiply by A and store in AP_prev
00395     //       AP_prev = A*P_prev
00396     //
00397     {
00398       Teuchos::TimeMonitor OpTimer(*_timerOp);
00399       _lp->ApplyOp( *P_prev, *AP_prev );
00400     }
00401     //
00402     // Compute initial residual block and store in 1st block of _residvecs
00403     //     R_prev = cur_block_rhs - A*P_prev
00404     //
00405     MVT::MvAddMv(one, *_cur_block_rhs, -one, *AP_prev, *R_prev );
00406     //
00407     //-------Compute and save the initial residual norms----------
00408     //
00409     MVT::MvNorm(*R_prev, &_init_resid_norms);
00410     //
00411     // Update indices of current (independent) blocks.
00412     // If a residual is too small, it will be dropped from
00413     // the current block, thus, from future computations
00414     //
00415     j = 0;
00416     ind_idx.resize( 0 );
00417     for (i=0; i<_blocksize; i++){
00418       _cur_resid_norms[i] = _init_resid_norms[i];
00419       if (_init_resid_norms[i] > _prec)
00420   ind_idx.push_back( i );
00421     }
00422     ind_blksz = ind_idx.size(); 
00423     //
00424     if (ind_blksz > 0) { 
00425       //
00426       // All initial residuals have not converged -- continue Block CG  
00427       // Compute the initial block of direciton vectors
00428       //
00429       // Associate current blocks of residuals, directions, and solution block
00430       // with R_prev, P_prev, and cur_sol
00431       //
00432       R_prev = MVT::CloneView( *_residvecs, ind_idx );
00433       P_prev = MVT::CloneView( *_basisvecs, ind_idx );
00434       //
00435       //----------------Compute initial direction vectors--------------------------
00436       // Initially, they are set to the preconditioned residuals
00437       //
00438       if ( isPrec ) {
00439   Teuchos::TimeMonitor PrecTimer(*_timerPrec);
00440   _lp->ApplyLeftPrec( *R_prev, *P_prev );
00441       } else {
00442   MVT::MvAddMv( one , *R_prev, zero, *R_prev, *P_prev); 
00443       }
00444       //
00445       // Compute an orthonormal block of initial direction vectors,
00446       // and check for dependencies, adjusting indices of independent
00447       // vectors if needed
00448       //
00449       Teuchos::SerialDenseMatrix<int,ScalarType> G(ind_blksz, ind_blksz);
00450       exit_flg = false;
00451       {
00452   Teuchos::TimeMonitor OrthoTimer(*_timerOrtho);
00453   exit_flg = QRFactorDef(*P_prev, G, &cols );
00454       }
00455       num_ind = cols.size();
00456       //
00457       if ( exit_flg ) {
00458   if (_om->isVerbosityAndPrint( Errors )) {
00459     *_os << " Exiting Block CG iteration " << endl; 
00460     *_os << " ERROR: No independent initial direction vectors" << endl;
00461   }
00462       }   
00463       if (num_ind < ind_blksz) {
00464   // The initial block of direction vectors are linearly dependent
00465   if (_om->isVerbosityAndPrint( Warnings )) {
00466     *_os << " Initial direction vectors are dependent" << endl;
00467     *_os << " Adjusting blocks and indices for iteration" << endl;
00468   }
00469   //
00470   // Adjust block and indices for iteration.
00471   //
00472   ind_blksz = num_ind;
00473   std::vector<int> temp_idx( ind_blksz );
00474   for (i=0; i< ind_blksz; i++)
00475     temp_idx[i] = ind_idx[cols[i]];
00476   ind_idx.resize( ind_blksz );
00477   for (i=0; i< ind_blksz; i++)      
00478     ind_idx[i] = temp_idx[i];
00479   
00480       }  // end if (num < ind_blksz)
00481     }  // end if (ind_blksz > 0)
00482     //    
00483     else {  // all initial residuals have converged
00484       if (_om->isVerbosityAndPrint( Warnings )) {
00485   *_os << " Exiting Block CG iteration " 
00486        << " -- Iteration# " << _iter << endl;
00487   *_os << "Reason: All initial residuals have converged" << endl;
00488       }
00489       exit_flg = true;
00490     }
00491     
00492     // ***************************************************************************
00493     // ************************Main CG Loop***************************************
00494     // ***************************************************************************
00495     // 
00496     _new_blk = 1;
00497     for (_iter=0; _stest->CheckStatus(this) == Unconverged && !exit_flg; _iter++) {
00498       //
00499       //----------------Compute the new blocks of iterates and residuals------------------
00500       //
00501       // Get views of the previous blocks of residuals, direction vectors, etc.
00502       //
00503       if (_new_blk){
00504   P_prev = MVT::CloneView( *_basisvecs, ind_idx );
00505       }
00506       else {
00507   index2.resize( ind_blksz );
00508   for (i=0; i< ind_blksz; i++) {
00509     index2[i] = _blocksize + ind_idx[i];
00510   } 
00511   P_prev = MVT::CloneView( *_basisvecs, index2 );
00512       }
00513       //
00514       index2.resize( _blocksize );
00515       for (i=0; i < _blocksize; i++){
00516   index2[i] = _blocksize + i;
00517       }
00518       if (_new_blk){
00519   R_prev = MVT::CloneView( *_residvecs, index );
00520   R_new = MVT::CloneView( *_residvecs, index2 );
00521       }
00522       else {
00523   R_prev = MVT::CloneView( *_residvecs, index2 );
00524   R_new = MVT::CloneView( *_residvecs, index );
00525       }
00526       //
00527       // Compute the coefficient matrix alpha
00528       //
00529       // P_prev^T * A * P_prev * alpha = P_prev^T * R_prev
00530       // 1) Compute P_prev^T * A * P_prev = T2 and P_prev^T * R_prev = T1
00531       // 2) Compute the Cholesky Factorization of T2
00532       // 3) Back and Forward Solves for alpha
00533       //
00534       // 1)
00535       if ( ind_blksz < prev_ind_blksz ) {  
00536   //
00537   // The number of independent direction vectors has changed,
00538   // so the dimension of the application multivectors needs to be resized.
00539   //
00540   AP_prev = MVT::CloneView( *temp_blk, ind_idx ); 
00541   alpha.reshape( ind_blksz, _blocksize );
00542   T2.reshape( ind_blksz, ind_blksz );
00543       }
00544       {
00545   Teuchos::TimeMonitor OpTimer(*_timerOp);
00546   _lp->ApplyOp( *P_prev, *AP_prev );
00547       }
00548       MVT::MvTransMv( one, *P_prev, *R_prev, alpha);   
00549       MVT::MvTransMv( one, *P_prev, *AP_prev, T2);
00550       //
00551       // 2)
00552       lapack.POTRF(UPLO, ind_blksz, T2.values(), ind_blksz, &info);
00553       if (info != 0) {
00554   if(_om->isVerbosityAndPrint( Errors )){
00555     *_os << " Exiting Block CG iteration "
00556          << " -- Iteration# " << _iter << endl;
00557     *_os << " ERROR: Cannot compute coefficient matrix alpha" << endl;
00558     *_os << " P_prev'* A*P_prev is singular" << endl;
00559     *_os << " Solution will be updated upon exiting loop" << endl;
00560   }
00561   break;
00562       }
00563       // 3)
00564       lapack.POTRS(UPLO, ind_blksz, _blocksize, T2.values(), ind_blksz, 
00565        alpha.values(), ind_blksz, &info);
00566       // Note: solution returned in alpha
00567       if (info != 0) {
00568   if(_om->isVerbosityAndPrint( Errors  )){
00569     *_os << " Exiting Block CG iteration "
00570          << " -- Iteration# " << _iter << endl;
00571     *_os << " ERROR: Cannot compute coefficient matrix alpha" << endl;
00572     *_os << " Solution will be updated upon exiting loop" << endl;
00573   }
00574   break;
00575       }
00576       //
00577       // Update the solution: cur_block_sol = cur_block_sol + P_prev*alpha
00578       // 
00579       MVT::MvTimesMatAddMv( one, *P_prev, alpha, one, *_cur_block_sol );
00580       _lp->SolutionUpdated();
00581       //
00582       // Update the residual vectors: R_new = R_prev - A*P_prev*alpha
00583       //
00584       MVT::MvAddMv(one, *R_prev, zero, *R_prev, *R_new);
00585       MVT::MvTimesMatAddMv(-one, *AP_prev, alpha, one, *R_new);
00586       //
00587       // ****Compute the Current Relative Residual Norms and the Block Error****
00588       //
00589       MVT::MvNorm( *R_new, &_cur_resid_norms );
00590       //
00591       prev_ind_blksz = ind_blksz; // Save old ind_blksz of P_prev
00592       //
00593       // Update the number of current residuals that correspond
00594       // to linearly independent direction vectors. Note that
00595       // ind_idx are a subset of cur_idx.
00596       //
00597       k = 0;
00598       std::vector<int> temp_idx;
00599       for (i=0; i< ind_blksz; i++){
00600   if (_cur_resid_norms[ ind_idx[i] ] / _init_resid_norms[ ind_idx[i] ] > _prec)
00601     temp_idx.push_back( ind_idx[i] );
00602       }
00603       ind_blksz = temp_idx.size();
00604       if (ind_blksz < prev_ind_blksz ) {
00605   ind_idx.resize( ind_blksz );
00606   for (i=0; i< ind_blksz; i++)
00607     ind_idx[i] = temp_idx[i];
00608       }
00609       //
00610       // ****************Print iteration information*****************************
00611       //
00612       if (_om->isVerbosity( IterationDetails )) {
00613   PrintCGIterInfo( ind_idx );
00614       }
00615       //
00616       // ****************Test for breakdown*************************************
00617       //
00618       if (ind_blksz <= 0){
00619   if (_om->isVerbosityAndPrint( Errors )) {
00620     *_os << " Exiting Block CG iteration " 
00621          << " -- Iteration# " << _iter << endl;
00622     *_os << " ERROR: No more independent direction vectors" << endl;
00623     *_os << " Solution will be updated upon exiting loop" << endl;
00624   }
00625   break;
00626       }
00627       //
00628       // **************Compute the new block of direction vectors****************
00629       //
00630       // Get views of the new blocks of independent direction vectors and
00631       // the corresponding residuals. 
00632       //
00633       index2.resize( ind_blksz );
00634       for (i=0; i<ind_blksz; i++){
00635   index2[i] = _blocksize + ind_idx[i];
00636       }
00637       
00638       if (_new_blk) {
00639   R_new = MVT::CloneView( *_residvecs, index2 );
00640   P_new = MVT::CloneView( *_basisvecs, index2 );
00641       }
00642       else {
00643   R_new = MVT::CloneView( *_residvecs, ind_idx );
00644   P_new = MVT::CloneView( *_basisvecs, ind_idx );
00645       }
00646       //
00647       // Put the current preconditioned initial residual into P_new 
00648       // since P_new = precond_resid + P_prev * beta
00649       //
00650       if ( isPrec ) {
00651   Teuchos::TimeMonitor PrecTimer(*_timerPrec);
00652   _lp->ApplyLeftPrec( *R_new, *P_new );
00653       } else { 
00654   MVT::MvAddMv( one, *R_new, zero, *R_new, *P_new ); 
00655       }
00656       // 
00657       // Compute coefficient matrix beta
00658       //
00659       // P_prev^T A * P_prev * beta = P_prev^T A * precond_resid
00660       // 1) Compute P_prev^T A * P_prev = T2 and P_prev^T * A * precond_resid = T3
00661       //                                     or (A*P_prev)^T * precond_resid (A SPD)
00662       // 2) Compute the Cholesky Factorization of T2
00663       // 3) Back and Forward Solves for beta
00664       //
00665       beta.reshape(prev_ind_blksz,ind_blksz);
00666       //
00667       // 1 & 2)  Note: we already have computed T2 and its Cholesky
00668       //         factorization during computation of alpha
00669       MVT::MvTransMv(-one, *AP_prev, *P_new, beta);
00670       // 3)
00671       lapack.POTRS(UPLO, prev_ind_blksz, ind_blksz, T2.values(), prev_ind_blksz, 
00672        beta.values(), prev_ind_blksz, &info);
00673       // Note: Solution returned in beta
00674       if (info != 0) {
00675   if (_om->isVerbosityAndPrint( Errors )) {
00676     *_os << " Exiting Block CG iteration " 
00677          << " -- Iteration# " << _iter << endl;
00678     *_os << "ERROR: Cannot compute coefficient matrix beta" << endl;
00679     *_os << "Solution will be updated upon exiting loop" << endl;
00680   }
00681   break;
00682       }
00683       //
00684       // Compute: P_new = precond_resid + P_prev * beta
00685       // ( remember P_new already = precond_resid )
00686       MVT::MvTimesMatAddMv(one, *P_prev, beta, one, *P_new);
00687       //
00688       // Check A-orthogonality of new and previous blocks of direction vectors
00689       //
00690       if (_om->isVerbosity( OrthoDetails )) {
00691   if (_om->doPrint()) *_os << "Orthogonality check" << endl;
00692   CheckAOrthogonality(*P_prev, *P_new);   
00693       }
00694       //
00695       // Compute orthonormal block of direction vectors,
00696       // and check for dependencies, adjusting indices of
00697       // independent vectors if needed
00698       //
00699       Teuchos::SerialDenseMatrix<int,ScalarType> G(ind_blksz,ind_blksz);
00700       {
00701   Teuchos::TimeMonitor OrthoTimer(*_timerOrtho);
00702   exit_flg = QRFactorDef(*P_new, G, &cols);
00703       }
00704       //
00705       // Check if the orthogonalization has failed.
00706       //
00707       num_ind = cols.size();
00708       if ( exit_flg ) {
00709   ind_blksz = num_ind;
00710   if (_om->isVerbosityAndPrint( Errors )) {
00711     *_os  << " Exiting Block CG iteration "  
00712     << " -- Iteration# " << _iter << endl;
00713     *_os << " ERROR: No more linearly independent direction vectors" << endl;
00714     *_os << " Solution will be updated upon exiting loop" << endl;
00715   }
00716   break;
00717       }
00718       //
00719       // Check if the new block of direction vectors are linearly dependent
00720       //
00721       if (num_ind < ind_blksz) {
00722   if (_om->isVerbosityAndPrint( OrthoDetails )) {
00723     *_os << "The new block of direction vectors are dependent " << endl;
00724     *_os << "# independent direction vectors: " << num_ind << endl;
00725     *_os << " Independent indices: " << endl;
00726     for (i=0; i<num_ind ; i++)
00727       *_os << cols[i] << " ";
00728     *_os << endl << endl;
00729   }
00730   ind_blksz = num_ind;
00731   std::vector<int> temp_idx( ind_blksz );
00732   for (i=0; i<ind_blksz; i++)
00733     temp_idx[i] = ind_idx[cols[i]];
00734   ind_idx.resize( ind_blksz );
00735   for (i=0; i<ind_blksz; i++)
00736     ind_idx[i] = temp_idx[i];
00737       }
00738       //
00739       // Check A-orthogonality after orthonormalization
00740       //
00741       if (_om->isVerbosity( OrthoDetails )) {
00742   if (_om->doPrint())  *_os << "Orthogonality check after orthonormalization " << endl; 
00743   CheckAOrthogonality(*P_prev, *P_new);
00744       }
00745       // 
00746       // *****Update index of new blocks*****************************************
00747       //
00748       if (_new_blk)
00749   _new_blk--;
00750       else
00751   _new_blk++;
00752       //
00753     } // end of the main CG loop -- for(_iter = 0;...)
00754     //
00755     // Inform the linear problem manager that we are done with the current block of linear systems.
00756     //
00757     _lp->SetCurrLSVec();
00758     //
00759     // Print out solver status.
00760     //
00761     if (_om->isVerbosityAndPrint( FinalSummary )) {
00762       _stest->Print(*_os);
00763   }  
00764     //
00765   } // end BlockIteration()
00766   //
00767   
00768   template<class ScalarType, class MV, class OP>
00769   bool 
00770   BlockCG<ScalarType,MV,OP>::QRFactorDef (MV& VecIn, 
00771             Teuchos::SerialDenseMatrix<int,ScalarType>& R, 
00772             std::vector<int>* cols) 
00773   {
00774     int i, j, k;
00775     int num_orth, num_dep = 0;
00776     int nb = MVT::GetNumberVecs( VecIn );
00777     std::vector<int> index, dep_idx;
00778     std::vector<int> index2( 1 );
00779     RefCountPtr<MV> qj, Qj;
00780     Teuchos::SerialDenseVector<int,ScalarType> rj;
00781     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00782     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00783     std::vector<MagnitudeType> NormVecIn( nb );
00784     std::vector<MagnitudeType> qNorm( 1 );
00785     //
00786     // Zero out the array that will contain the Fourier coefficients.
00787     // 
00788     R.putScalar();
00789     //
00790     // Reset the size of the independent columns back to zero.
00791     //
00792     cols->resize( 0 );
00793     //
00794     // Compute the Frobenius norm of VecIn -- this will be used to
00795     // determine rank deficiency of VecIn
00796     //
00797     MVT::MvNorm( VecIn, &NormVecIn );
00798     ScalarType FroNorm = zero;
00799     for (j=0; j<nb; j++) {
00800       FroNorm += NormVecIn[j] * NormVecIn[j];
00801     }
00802     FroNorm = Teuchos::ScalarTraits<ScalarType>::squareroot(FroNorm);
00803     //
00804     // Start the loop to orthogonalize the nb columns of VecIn.
00805     //
00806     //
00807     for ( j=0; j<nb; j++ ) {
00808       //
00809       // Grab the j-th column of VecIn (the first column is indexed to 
00810       // be the zero-th one).
00811       //
00812       index2[0] = j;
00813       qj = MVT::CloneView( VecIn, index2 );
00814       if ( j ) {
00815   //
00816   // Grab the first j columns of VecIn (that are now an orthogonal
00817   // basis for first j columns of the entering VecIn).
00818   //
00819   index.resize( j );
00820   for (i=0; i<j; i++)
00821     index[i] = i;
00822   Qj = MVT::CloneView( VecIn, index );
00823   rj.size(j);
00824   //
00825   // Enter a for loop that does two (num_orth) steps of classical 
00826   // Gram-Schmidt orthogonalization.
00827   //
00828   for (num_orth=0; num_orth<2; num_orth++) {
00829     //
00830     // Determine the Fourier coefficients for orthogonalizing column
00831     // j of VecIn against columns 0:j-1 of VecIn. In other words, 
00832     // result = trans(Qj)*qj.
00833     //
00834     MVT::MvTransMv( one, *Qj, *qj, rj );
00835     //
00836     // Sum result[0:j-1] into column j of R.
00837     //
00838     for (k=0; k<num_dep; k++) {
00839       rj[dep_idx[k]] = zero;
00840     }
00841     //
00842     for ( k=0; k<j; k++ ) {
00843       R(k,j) += rj[k];
00844     }
00845     //
00846     //   Compute qj <- qj - Qj * rj.
00847     //
00848     MVT::MvTimesMatAddMv(-one, *Qj, rj, one, *qj);
00849   }
00850   
00851       } // if (j)
00852       //
00853       // Compute the norm of column j of VecIn (=qj).
00854       //
00855       MVT::MvNorm( *qj, &qNorm );
00856       R(j,j) = qNorm[0];
00857       //
00858       if ( NormVecIn[j] > _prec && SCT::magnitude(R(j,j)) > (_prec * NormVecIn[j]) ) {
00859   //
00860   // Normalize qj to make it into a unit vector.
00861   //
00862   ScalarType rjj = one / R(j,j);
00863   MVT::MvAddMv( rjj, *qj, zero, *qj, *qj );
00864   cols->push_back( j );  
00865       }
00866       else {
00867   // 
00868   if (_om->isVerbosityAndPrint( OrthoDetails )){
00869     *_os << "Rank deficiency at column index: " << j << endl;
00870   }
00871   //
00872   // Don't normalize qj, enter one on diagonal of R,
00873   // and zeros in the row to the right of the diagonal -- this
00874   // requires updating the indices of the dependent columns
00875   //
00876   R(j,j) = one;
00877   dep_idx.push_back( j );
00878   num_dep++;
00879       } 
00880       
00881     } // for (j=0; j<nb; j++)
00882     //
00883     // Return true if we could not create any independent direction vectors (failure).
00884     //
00885     if (cols->size() == 0) return true;
00886     return false;
00887     //
00888   } // end QRFactorDef
00889   
00890   
00891   template<class ScalarType, class MV, class OP>
00892   void 
00893   BlockCG<ScalarType,MV,OP>::CheckOrthogonality(MV& P1, MV& P2) 
00894   {
00895     //
00896     // This routine computes P2^T * P1
00897     //
00898     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00899     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00900     int i, k;
00901     //
00902     int numvecs1 = MVT::GetNumberVecs( P1 );
00903     int numvecs2 = MVT::GetNumberVecs( P2 );
00904     //
00905     Teuchos::SerialDenseMatrix<int,ScalarType> PtP(numvecs2, numvecs1);
00906     MVT::MvTransMv( one, P2, P1, PtP);
00907     //
00908     ScalarType* ptr = PtP.values();
00909     ScalarType column_sum;
00910     //
00911     for (k=0; k<numvecs1; k++) {
00912       column_sum = zero;
00913       for (i=0; i<numvecs2; i++) {
00914   column_sum += ptr[i];
00915       }
00916       *_os << " P2^T*P1 for column "
00917      << k << " is  " << Teuchos::ScalarTraits<ScalarType>::magnitude(column_sum) << endl;
00918       ptr += numvecs2;
00919     }
00920     *_os << endl;
00921     //
00922   } // end check_orthog
00923   //
00924   
00925   template<class ScalarType, class MV, class OP>
00926   void 
00927   BlockCG<ScalarType,MV,OP>::CheckAOrthogonality(MV& P1, MV& P2) 
00928   {
00929     //
00930     // This routine computes P2^T * A * P1
00931     // Checks the orthogonality wrt A between any two blocks of multivectors with the same length
00932     //
00933     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00934     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00935     int i, k;
00936     //
00937     int numvecs1 = MVT::GetNumberVecs( P1 );
00938     int numvecs2 = MVT::GetNumberVecs( P2 );
00939     //
00940     RefCountPtr<MV> AP = MVT::Clone( P1, numvecs1 );
00941     {
00942       Teuchos::TimeMonitor OpTimer(*_timerOp);
00943       _lp->ApplyOp(P1, *AP);
00944     }
00945     //
00946     Teuchos::SerialDenseMatrix<int,ScalarType> PAP(numvecs2, numvecs1);
00947     MVT::MvTransMv( one, P2, *AP, PAP);
00948     //
00949     ScalarType* ptr = PAP.values();
00950     ScalarType column_sum;
00951     //
00952     for (k=0; k<numvecs1; k++) {
00953       column_sum = zero;
00954       for (i=0; i<numvecs2; i++) {
00955   column_sum += ptr[i];
00956     }
00957       *_os << " P2^T*A*P1 for column "
00958      << k << " is  " << Teuchos::ScalarTraits<ScalarType>::magnitude(column_sum) << endl;
00959       ptr += numvecs2;
00960     }
00961     *_os << " " << endl;
00962     //
00963   } // end check_orthog
00964   
00965   
00966   template<class ScalarType, class MV, class OP>
00967   void 
00968   BlockCG<ScalarType,MV,OP>::PrintCGIterInfo( const std::vector<int> &ind )
00969   {
00970     //
00971     int i;
00972     int indsz = ind.size();
00973     *_os << "# of independent direction vectors: " << indsz << endl;    
00974     *_os << " Independent indices: " << endl;
00975     for (i=0; i<indsz; i++){
00976       *_os << ind[i] << " ";
00977     }
00978     *_os << endl << endl;
00979     //
00980   } // end Print_CGiter_info
00981 
00982 
00983   // Overridden from Teuchos::Describable
00984   template<class ScalarType, class MV, class OP>
00985   std::string 
00986   BlockCG<ScalarType,MV,OP>::description() const
00987   {
00988     std::ostringstream oss;
00989     oss << "Belos::BlockCG<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00990     oss << "{";
00991     oss << "}";
00992     return oss.str();
00993   }
00994 //                                 
00995 } // end namespace Belos
00996 
00997 #endif
00998 // End of file BelosBlockCG.hpp
00999 
01000 

Generated on Thu Sep 18 12:30:12 2008 for Belos by doxygen 1.3.9.1