Belos Version of the Day
BelosDGKSOrthoManager.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 
00047 #ifndef BELOS_DGKS_ORTHOMANAGER_HPP
00048 #define BELOS_DGKS_ORTHOMANAGER_HPP
00049 
00057 // #define ORTHO_DEBUG
00058 
00059 #include "BelosConfigDefs.hpp"
00060 #include "BelosMultiVecTraits.hpp"
00061 #include "BelosOperatorTraits.hpp"
00062 #include "BelosMatOrthoManager.hpp"
00063 
00064 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00065 #include "Teuchos_TimeMonitor.hpp"
00066 #endif // BELOS_TEUCHOS_TIME_MONITOR
00067 
00068 namespace Belos {
00069 
00073   template<class ScalarType>
00074   Teuchos::RCP<const Teuchos::ParameterList> 
00075   getDefaultDgksParameters()
00076   {
00077     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00078     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00079 
00080     // This part makes this class method non-reentrant.
00081     static Teuchos::RCP<Teuchos::ParameterList> params;
00082     if (! params.is_null())
00083       return params;
00084 
00085     // Default parameter values for DGKS orthogonalization.
00086     // Documentation will be embedded in the parameter list.
00087     const int defaultMaxNumOrthogPasses = 2;
00088     const magnitude_type one = STM::one();
00089     const magnitude_type eps = STM::eps();
00090     const magnitude_type squareRootTwo = STM::squareroot(one + one);
00091     const magnitude_type defaultBlkTol = magnitude_type(10) * STM::squareroot(eps);
00092     const magnitude_type defaultDepTol = one / squareRootTwo;
00093     const magnitude_type defaultSingTol = magnitude_type(10) * eps;
00094 
00095     params = Teuchos::parameterList();
00096     params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00097      "Maximum number of orthogonalization passes "
00098      "(includes the first).  Default is 2, since "
00099      "\"twice is enough\" for Krylov methods.");
00100     params->set ("blkTol", defaultBlkTol, 
00101      "Block reorthogonalization threshhold.");
00102     params->set ("depTol", defaultDepTol, 
00103      "(Non-block) reorthogonalization threshold.");
00104     params->set ("singTol", defaultSingTol, 
00105      "Singular block detection threshold.");
00106     return params;
00107   }
00108 
00112   template<class ScalarType>
00113   Teuchos::RCP<const Teuchos::ParameterList> 
00114   getFastDgksParameters()
00115   {
00116     using Teuchos::ParameterList;
00117     using Teuchos::RCP;
00118     using Teuchos::rcp;
00119     using Teuchos::ScalarTraits;
00120     typedef typename ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00121     typedef ScalarTraits<magnitude_type> STM;
00122 
00123     // This part makes this class method non-reentrant.
00124     static RCP<ParameterList> params;
00125     if (params.is_null())
00126       {
00127   RCP<const ParameterList> defaultParams = getDefaultDgksParameters<ScalarType>();
00128   // Start with a clone of the default parameters
00129   params = rcp (new ParameterList (*defaultParams));
00130 
00131   const int maxBlkOrtho = 1;
00132   params->set ("maxNumOrthogPasses", maxBlkOrtho);
00133 
00134   const magnitude_type blkTol = STM::zero();
00135   params->set ("blkTol", blkTol);
00136 
00137   const magnitude_type depTol = STM::zero();
00138   params->set ("depTol", depTol);
00139 
00140   const magnitude_type singTol = STM::zero();
00141   params->set ("singTol", singTol);
00142       }
00143     return params;
00144   }
00145 
00151   template<class ScalarType>
00152   void
00153   readDgksParameters (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00154           int& maxNumOrthogPasses,
00155           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& blkTol,
00156           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& depTol,
00157           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& singTol)
00158   {
00159     using Teuchos::ParameterList;
00160     using Teuchos::RCP;
00161     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00162     const magnitude_type zero = Teuchos::ScalarTraits<magnitude_type>::zero();
00163     RCP<const ParameterList> defaultParams = getDefaultDgksParameters<ScalarType>();
00164 
00165     // Using temporary variables and fetching all values before
00166     // setting the output arguments ensures the strong exception
00167     // guarantee for this function: if an exception is thrown, no
00168     // externally visible side effects (in this case, setting the
00169     // output arguments) have taken place.
00170     int _maxNumOrthogPasses;
00171     magnitude_type _blkTol, _depTol, _singTol;
00172     if (params.is_null())
00173       {
00174   _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00175   _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00176   _depTol = defaultParams->get<magnitude_type> ("depTol");  
00177   _singTol = defaultParams->get<magnitude_type> ("singTol");
00178       }
00179     else
00180       {
00181   try {
00182     _maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00183     if (_maxNumOrthogPasses < 1)
00184       _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00185   } catch (Teuchos::Exceptions::InvalidParameter&) {
00186     _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00187   }
00188 
00189   try {
00190     _blkTol = params->get<magnitude_type> ("blkTol");
00191     if (_blkTol < zero)
00192       _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00193   } catch (Teuchos::Exceptions::InvalidParameter&) {
00194     _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00195   }
00196 
00197   try {
00198     _depTol = params->get<magnitude_type> ("depTol");
00199     if (_depTol < zero)
00200       _depTol = defaultParams->get<magnitude_type> ("depTol");
00201   } catch (Teuchos::Exceptions::InvalidParameter&) {
00202     _depTol = defaultParams->get<magnitude_type> ("depTol");
00203   }
00204 
00205   try {
00206     _singTol = params->get<magnitude_type> ("singTol");
00207     if (_singTol < zero)
00208       _singTol = defaultParams->get<magnitude_type> ("singTol");
00209   } catch (Teuchos::Exceptions::InvalidParameter&) {
00210     _singTol = defaultParams->get<magnitude_type> ("singTol");
00211   }
00212       }
00213     maxNumOrthogPasses = _maxNumOrthogPasses;
00214     blkTol = _blkTol;
00215     depTol = _depTol;
00216     singTol = _singTol;
00217   }
00218 
00219   template<class ScalarType, class MV, class OP>
00220   class DGKSOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00221 
00222   private:
00223     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00224     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00225     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00226     typedef MultiVecTraits<ScalarType,MV>      MVT;
00227     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00228 
00229   public:
00230     
00232 
00233 
00234     DGKSOrthoManager( const std::string& label = "Belos",
00235                       Teuchos::RCP<const OP> Op = Teuchos::null,
00236           const int max_blk_ortho = 2,
00237           const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00238           const MagnitudeType dep_tol = MGT::one()/MGT::squareroot( 2*MGT::one() ),
00239           const MagnitudeType sing_tol = 10*MGT::eps() )
00240       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00241   max_blk_ortho_( max_blk_ortho ),
00242   blk_tol_( blk_tol ),
00243   dep_tol_( dep_tol ),
00244   sing_tol_( sing_tol ),
00245   label_( label )
00246     {
00247         std::string orthoLabel = label_ + ": Orthogonalization";
00248 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00249         timerOrtho_ = Teuchos::TimeMonitor::getNewTimer( orthoLabel );
00250 #endif
00251     }    
00252 
00254     ~DGKSOrthoManager() {}
00256 
00257 
00259 
00260 
00262     void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
00263 
00265     void setDepTol( const MagnitudeType dep_tol ) { dep_tol_ = dep_tol; }
00266 
00268     void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
00269 
00271     MagnitudeType getBlkTol() const { return blk_tol_; } 
00272 
00274     MagnitudeType getDepTol() const { return dep_tol_; } 
00275 
00277     MagnitudeType getSingTol() const { return sing_tol_; } 
00278 
00280 
00281 
00283 
00284 
00312     void project ( MV &X, Teuchos::RCP<MV> MX, 
00313                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00314                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00315 
00316 
00319     void project ( MV &X, 
00320                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00321                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00322       project(X,Teuchos::null,C,Q);
00323     }
00324 
00325 
00326  
00351     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00352                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00353 
00354 
00357     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00358       return normalize(X,Teuchos::null,B);
00359     }
00360 
00361 
00417     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00418                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00419                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00420                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00421 
00424     int projectAndNormalize ( MV &X, 
00425                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00426                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00427                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const {
00428       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00429     }
00430 
00432 
00434 
00435 
00441     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00442     orthonormError(const MV &X) const {
00443       return orthonormError(X,Teuchos::null);
00444     }
00445 
00452     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00453     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00454 
00458     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00459     orthogError(const MV &X1, const MV &X2) const {
00460       return orthogError(X1,Teuchos::null,X2);
00461     }
00462 
00467     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00468     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00469 
00471 
00473 
00474 
00477     void setLabel(const std::string& label);
00478 
00481     const std::string& getLabel() const { return label_; }
00482 
00484 
00485   private:
00486     
00488     int max_blk_ortho_;
00489     MagnitudeType blk_tol_;
00490     MagnitudeType dep_tol_;
00491     MagnitudeType sing_tol_;
00492 
00494     std::string label_;
00495 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00496     Teuchos::RCP<Teuchos::Time> timerOrtho_;
00497 #endif // BELOS_TEUCHOS_TIME_MONITOR
00498 
00500     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00501       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00502       bool completeBasis, int howMany = -1 ) const;
00503     
00505     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00506         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00507         Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00508 
00522     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00523            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00524            Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00525            Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;    
00526   };
00527 
00529   // Set the label for this orthogonalization manager and create new timers if it's changed
00530   template<class ScalarType, class MV, class OP>
00531   void DGKSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00532   {
00533     if (label != label_) {
00534       label_ = label;
00535       std::string orthoLabel = label_ + ": Orthogonalization";
00536 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00537       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00538 #endif
00539     }
00540   }
00541 
00543   // Compute the distance from orthonormality
00544   template<class ScalarType, class MV, class OP>
00545   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00546   DGKSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00547     const ScalarType ONE = SCT::one();
00548     int rank = MVT::GetNumberVecs(X);
00549     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00550     MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
00551     for (int i=0; i<rank; i++) {
00552       xTx(i,i) -= ONE;
00553     }
00554     return xTx.normFrobenius();
00555   }
00556 
00558   // Compute the distance from orthogonality
00559   template<class ScalarType, class MV, class OP>
00560   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00561   DGKSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00562     int r1 = MVT::GetNumberVecs(X1);
00563     int r2  = MVT::GetNumberVecs(X2);
00564     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00565     MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
00566     return xTx.normFrobenius();
00567   }
00568 
00570   // Find an Op-orthonormal basis for span(X) - span(W)
00571   template<class ScalarType, class MV, class OP>
00572   int DGKSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00573                                     MV &X, Teuchos::RCP<MV> MX, 
00574                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00575                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00576                                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const 
00577   {
00578     using Teuchos::Array;
00579     using Teuchos::null;
00580     using Teuchos::is_null;
00581     using Teuchos::RCP;
00582     using Teuchos::rcp;
00583     using Teuchos::SerialDenseMatrix;
00584     typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
00585     typedef typename Array< RCP< const MV > >::size_type size_type;
00586     
00587 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00588     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00589 #endif
00590 
00591     ScalarType    ONE  = SCT::one();
00592     ScalarType    ZERO  = SCT::zero();
00593 
00594     int nq = Q.size();
00595     int xc = MVT::GetNumberVecs( X );
00596     int xr = MVT::GetVecLength( X );
00597     int rank = xc;
00598 
00599     // If the user doesn't want to store the normalization
00600     // coefficients, allocate some local memory for them.  This will
00601     // go away at the end of this method.
00602     if (is_null (B)) {
00603       B = rcp (new serial_dense_matrix_type (xc, xc));
00604     }
00605     // Likewise, if the user doesn't want to store the projection
00606     // coefficients, allocate some local memory for them.  Also make
00607     // sure that all the entries of C are the right size.  We're going
00608     // to overwrite them anyway, so we don't have to worry about the
00609     // contents (other than to resize them if they are the wrong
00610     // size).
00611     if (C.size() < nq)
00612       C.resize (nq);
00613     for (size_type k = 0; k < nq; ++k)
00614       {
00615   const int numRows = MVT::GetNumberVecs (*Q[k]);
00616   const int numCols = xc; // Number of vectors in X
00617   
00618   if (is_null (C[k]))
00619     C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
00620   else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
00621     {
00622       int err = C[k]->reshape (numRows, numCols);
00623       TEST_FOR_EXCEPTION(err != 0, std::runtime_error, 
00624              "DGKS orthogonalization: failed to reshape "
00625              "C[" << k << "] (the array of block "
00626              "coefficients resulting from projecting X "
00627              "against Q[1:" << nq << "]).");
00628     }
00629       }
00630 
00631     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00632     if (this->_hasOp) {
00633       if (MX == Teuchos::null) {
00634         // we need to allocate space for MX
00635         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00636         OPT::Apply(*(this->_Op),X,*MX);
00637       }
00638     }
00639     else {
00640       // Op == I  -->  MX = X (ignore it if the user passed it in)
00641       MX = Teuchos::rcp( &X, false );
00642     }
00643 
00644     int mxc = MVT::GetNumberVecs( *MX );
00645     int mxr = MVT::GetVecLength( *MX );
00646 
00647     // short-circuit
00648     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
00649 
00650     int numbas = 0;
00651     for (int i=0; i<nq; i++) {
00652       numbas += MVT::GetNumberVecs( *Q[i] );
00653     }
00654 
00655     // check size of B
00656     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00657                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00658     // check size of X and MX
00659     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00660                         "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00661     // check size of X w.r.t. MX 
00662     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00663                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00664     // check feasibility
00665     //TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00666     //                    "Belos::DGKSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00667 
00668     // Some flags for checking dependency returns from the internal orthogonalization methods
00669     bool dep_flg = false;
00670 
00671     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00672     Teuchos::RCP<MV> tmpX, tmpMX;
00673     tmpX = MVT::CloneCopy(X);
00674     if (this->_hasOp) {
00675       tmpMX = MVT::CloneCopy(*MX);
00676     }
00677 
00678     // Use the cheaper block orthogonalization.
00679     dep_flg = blkOrtho( X, MX, C, Q );
00680 
00681     // If a dependency has been detected in this block, then perform
00682     // the more expensive single-vector orthogonalization.
00683     if (dep_flg) {
00684       rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00685 
00686       // Copy tmpX back into X.
00687       MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00688       if (this->_hasOp) {
00689   MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00690       }
00691     } 
00692     else {
00693       // There is no dependency, so orthonormalize new block X
00694       rank = findBasis( X, MX, B, false );
00695       if (rank < xc) {
00696   // A dependency was found during orthonormalization of X,
00697   // rerun orthogonalization using more expensive single-vector orthogonalization.
00698   rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00699 
00700   // Copy tmpX back into X.
00701   MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00702   if (this->_hasOp) {
00703     MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00704   }
00705       }    
00706     }
00707 
00708     // this should not raise an std::exception; but our post-conditions oblige us to check
00709     TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00710                         "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00711 
00712     // Return the rank of X.
00713     return rank;
00714   }
00715 
00716 
00717 
00719   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00720   template<class ScalarType, class MV, class OP>
00721   int DGKSOrthoManager<ScalarType, MV, OP>::normalize(
00722                                 MV &X, Teuchos::RCP<MV> MX, 
00723                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00724 
00725 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00726     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00727 #endif
00728 
00729     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00730     return findBasis(X, MX, B, true);
00731   }
00732 
00733 
00734 
00736   template<class ScalarType, class MV, class OP>
00737   void DGKSOrthoManager<ScalarType, MV, OP>::project(
00738                           MV &X, Teuchos::RCP<MV> MX, 
00739                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00740                           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00741     // For the inner product defined by the operator Op or the identity (Op == 0)
00742     //   -> Orthogonalize X against each Q[i]
00743     // Modify MX accordingly
00744     //
00745     // Note that when Op is 0, MX is not referenced
00746     //
00747     // Parameter variables
00748     //
00749     // X  : Vectors to be transformed
00750     //
00751     // MX : Image of the block vector X by the mass matrix
00752     //
00753     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00754     //
00755 
00756 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00757     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00758 #endif
00759 
00760     int xc = MVT::GetNumberVecs( X );
00761     int xr = MVT::GetVecLength( X );
00762     int nq = Q.size();
00763     std::vector<int> qcs(nq);
00764     // short-circuit
00765     if (nq == 0 || xc == 0 || xr == 0) {
00766       return;
00767     }
00768     int qr = MVT::GetVecLength ( *Q[0] );
00769     // if we don't have enough C, expand it with null references
00770     // if we have too many, resize to throw away the latter ones
00771     // if we have exactly as many as we have Q, this call has no effect
00772     C.resize(nq);
00773 
00774 
00775     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00776     if (this->_hasOp) {
00777       if (MX == Teuchos::null) {
00778         // we need to allocate space for MX
00779         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00780         OPT::Apply(*(this->_Op),X,*MX);
00781       }
00782     }
00783     else {
00784       // Op == I  -->  MX = X (ignore it if the user passed it in)
00785       MX = Teuchos::rcp( &X, false );
00786     }
00787     int mxc = MVT::GetNumberVecs( *MX );
00788     int mxr = MVT::GetVecLength( *MX );
00789 
00790     // check size of X and Q w.r.t. common sense
00791     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00792                         "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00793     // check size of X w.r.t. MX and Q
00794     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00795                         "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
00796 
00797     // tally up size of all Q and check/allocate C
00798     int baslen = 0;
00799     for (int i=0; i<nq; i++) {
00800       TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00801                           "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
00802       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00803       TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00804                           "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
00805       baslen += qcs[i];
00806 
00807       // check size of C[i]
00808       if ( C[i] == Teuchos::null ) {
00809         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00810       }
00811       else {
00812         TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00813                            "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
00814       }
00815     }
00816 
00817     // Use the cheaper block orthogonalization, don't check for rank deficiency.
00818     blkOrtho( X, MX, C, Q );
00819 
00820   }  
00821 
00823   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00824   // the rank is numvectors(X)
00825   template<class ScalarType, class MV, class OP>
00826   int DGKSOrthoManager<ScalarType, MV, OP>::findBasis(
00827                   MV &X, Teuchos::RCP<MV> MX, 
00828                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00829                   bool completeBasis, int howMany ) const {
00830     // For the inner product defined by the operator Op or the identity (Op == 0)
00831     //   -> Orthonormalize X 
00832     // Modify MX accordingly
00833     //
00834     // Note that when Op is 0, MX is not referenced
00835     //
00836     // Parameter variables
00837     //
00838     // X  : Vectors to be orthonormalized
00839     //
00840     // MX : Image of the multivector X under the operator Op
00841     //
00842     // Op  : Pointer to the operator for the inner product
00843     //
00844     //
00845 
00846     const ScalarType ONE  = SCT::one();
00847     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00848 
00849     int xc = MVT::GetNumberVecs( X );
00850     int xr = MVT::GetVecLength( X );
00851 
00852     if (howMany == -1) {
00853       howMany = xc;
00854     }
00855 
00856     /*******************************************************
00857      *  If _hasOp == false, we will not reference MX below *
00858      *******************************************************/
00859 
00860     // if Op==null, MX == X (via pointer)
00861     // Otherwise, either the user passed in MX or we will allocated and compute it
00862     if (this->_hasOp) {
00863       if (MX == Teuchos::null) {
00864         // we need to allocate space for MX
00865         MX = MVT::Clone(X,xc);
00866         OPT::Apply(*(this->_Op),X,*MX);
00867       }
00868     }
00869 
00870     /* if the user doesn't want to store the coefficienets, 
00871      * allocate some local memory for them 
00872      */
00873     if ( B == Teuchos::null ) {
00874       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00875     }
00876 
00877     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00878     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00879 
00880     // check size of C, B
00881     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00882                         "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
00883     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00884                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
00885     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00886                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
00887     TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00888                         "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
00889     TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
00890                         "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
00891 
00892     /* xstart is which column we are starting the process with, based on howMany
00893      * columns before xstart are assumed to be Op-orthonormal already
00894      */
00895     int xstart = xc - howMany;
00896 
00897     for (int j = xstart; j < xc; j++) {
00898 
00899       // numX is 
00900       // * number of currently orthonormal columns of X
00901       // * the index of the current column of X
00902       int numX = j;
00903       bool addVec = false;
00904 
00905       // Get a view of the vector currently being worked on.
00906       std::vector<int> index(1);
00907       index[0] = numX;
00908       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
00909       Teuchos::RCP<MV> MXj;
00910       if ((this->_hasOp)) {
00911         // MXj is a view of the current vector in MX
00912         MXj = MVT::CloneViewNonConst( *MX, index );
00913       }
00914       else {
00915         // MXj is a pointer to Xj, and MUST NOT be modified
00916         MXj = Xj;
00917       }
00918 
00919       // Get a view of the previous vectors.
00920       std::vector<int> prev_idx( numX );
00921       Teuchos::RCP<const MV> prevX, prevMX;
00922 
00923       if (numX > 0) {
00924         for (int i=0; i<numX; i++) {
00925           prev_idx[i] = i;
00926         }
00927         prevX = MVT::CloneView( X, prev_idx );
00928         if (this->_hasOp) {
00929           prevMX = MVT::CloneView( *MX, prev_idx );
00930         }
00931       } 
00932 
00933       // Make storage for these Gram-Schmidt iterations.
00934       Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
00935       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00936       //
00937       // Save old MXj vector and compute Op-norm
00938       //
00939       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00940       MVT::MvDot( *Xj, *MXj, oldDot );
00941       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
00942       TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
00943         "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
00944 
00945       if (numX > 0) {
00946   // Apply the first step of Gram-Schmidt
00947   
00948   // product <- prevX^T MXj
00949   MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,product);
00950   
00951   // Xj <- Xj - prevX prevX^T MXj   
00952   //     = Xj - prevX product
00953   MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
00954   
00955   // Update MXj
00956   if (this->_hasOp) {
00957     // MXj <- Op*Xj_new
00958     //      = Op*(Xj_old - prevX prevX^T MXj)
00959     //      = MXj - prevMX product
00960     MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
00961   }
00962   
00963   // Compute new Op-norm
00964   MVT::MvDot( *Xj, *MXj, newDot );
00965   
00966   // Check if a correction is needed.
00967   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(dep_tol_*oldDot[0]) ) {
00968     // Apply the second step of Gram-Schmidt
00969     // This is the same as above
00970     Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
00971     
00972     MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
00973     product += P2;
00974     MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00975     if ((this->_hasOp)) {
00976       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00977     }
00978   } // if (newDot[0] < dep_tol_*oldDot[0])
00979   
00980       } // if (numX > 0)
00981 
00982         // Compute Op-norm with old MXj
00983       MVT::MvDot( *Xj, *oldMXj, newDot );
00984 
00985       // Check to see if the new vector is dependent.
00986       if (completeBasis) {
00987   //
00988   // We need a complete basis, so add random vectors if necessary
00989   //
00990   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
00991     
00992     // Add a random vector and orthogonalize it against previous vectors in block.
00993     addVec = true;
00994 #ifdef ORTHO_DEBUG
00995     std::cout << "Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
00996 #endif
00997     //
00998     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
00999     Teuchos::RCP<MV> tempMXj;
01000     MVT::MvRandom( *tempXj );
01001     if (this->_hasOp) {
01002       tempMXj = MVT::Clone( X, 1 );
01003       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01004     } 
01005     else {
01006       tempMXj = tempXj;
01007     }
01008     MVT::MvDot( *tempXj, *tempMXj, oldDot );
01009     //
01010     for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
01011       MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,product);
01012       MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
01013       if (this->_hasOp) {
01014         MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
01015       }
01016     }
01017     // Compute new Op-norm
01018     MVT::MvDot( *tempXj, *tempMXj, newDot );
01019     //
01020     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01021       // Copy vector into current column of _basisvecs
01022       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
01023       if (this->_hasOp) {
01024         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
01025       }
01026     }
01027     else {
01028       return numX;
01029     } 
01030   }
01031       }
01032       else {
01033   //
01034   // We only need to detect dependencies.
01035   //
01036   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
01037     return numX;
01038   }
01039       }
01040       
01041       // If we haven't left this method yet, then we can normalize the new vector Xj.
01042       // Normalize Xj.
01043       // Xj <- Xj / std::sqrt(newDot)
01044       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01045 
01046       if (SCT::magnitude(diag) > ZERO) {      
01047         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01048         if (this->_hasOp) {
01049     // Update MXj.
01050     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01051         }
01052       }
01053 
01054       // If we've added a random vector, enter a zero in the j'th diagonal element.
01055       if (addVec) {
01056   (*B)(j,j) = ZERO;
01057       }
01058       else {
01059   (*B)(j,j) = diag;
01060       }
01061 
01062       // Save the coefficients, if we are working on the original vector and not a randomly generated one
01063       if (!addVec) {
01064   for (int i=0; i<numX; i++) {
01065     (*B)(i,j) = product(i,0);
01066   }
01067       }
01068 
01069     } // for (j = 0; j < xc; ++j)
01070 
01071     return xc;
01072   }
01073 
01075   // Routine to compute the block orthogonalization
01076   template<class ScalarType, class MV, class OP>
01077   bool 
01078   DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
01079                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01080                Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01081   {
01082     int nq = Q.size();
01083     int xc = MVT::GetNumberVecs( X );
01084     bool dep_flg = false;
01085     const ScalarType ONE  = SCT::one();
01086 
01087     std::vector<int> qcs( nq );
01088     for (int i=0; i<nq; i++) {
01089       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01090     }
01091 
01092     // Perform the Gram-Schmidt transformation for a block of vectors
01093     
01094     // Compute the initial Op-norms
01095     std::vector<ScalarType> oldDot( xc );
01096     MVT::MvDot( X, *MX, oldDot );
01097 
01098     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01099     // Define the product Q^T * (Op*X)
01100     for (int i=0; i<nq; i++) {
01101       // Multiply Q' with MX
01102       MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
01103       // Multiply by Q and subtract the result in X
01104       MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
01105 
01106       // Update MX, with the least number of applications of Op as possible
01107       if (this->_hasOp) {
01108         if (xc <= qcs[i]) {
01109           OPT::Apply( *(this->_Op), X, *MX);
01110         }
01111         else {
01112           // this will possibly be used again below; don't delete it
01113           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01114           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01115           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01116         }
01117       }
01118     }
01119 
01120     // Do as many steps of classical Gram-Schmidt as required by max_blk_ortho_
01121     for (int j = 1; j < max_blk_ortho_; ++j) {
01122       
01123       for (int i=0; i<nq; i++) {
01124   Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
01125         
01126   // Apply another step of classical Gram-Schmidt
01127   MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
01128   *C[i] += C2;
01129   MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
01130         
01131   // Update MX, with the least number of applications of Op as possible
01132   if (this->_hasOp) {
01133     if (MQ[i].get()) {
01134       // MQ was allocated and computed above; use it
01135       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01136     }
01137     else if (xc <= qcs[i]) {
01138       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01139       OPT::Apply( *(this->_Op), X, *MX);
01140     }
01141   }
01142       } // for (int i=0; i<nq; i++)
01143     } // for (int j = 0; j < max_blk_ortho; ++j)
01144   
01145     // Compute new Op-norms
01146     std::vector<ScalarType> newDot(xc);
01147     MVT::MvDot( X, *MX, newDot );
01148  
01149     // Check to make sure the new block of vectors are not dependent on previous vectors
01150     for (int i=0; i<xc; i++){
01151       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01152   dep_flg = true;
01153   break;
01154       }
01155     } // end for (i=0;...)
01156 
01157     return dep_flg;
01158   }
01159   
01160 
01161   template<class ScalarType, class MV, class OP>
01162   int
01163   DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01164                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01165                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01166                    Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
01167   {
01168     Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
01169 
01170     const ScalarType ONE  = SCT::one();
01171     const ScalarType ZERO  = SCT::zero();
01172     
01173     int nq = Q.size();
01174     int xc = MVT::GetNumberVecs( X );
01175     std::vector<int> indX( 1 );
01176     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01177 
01178     std::vector<int> qcs( nq );
01179     for (int i=0; i<nq; i++) {
01180       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01181     }
01182 
01183     // Create pointers for the previous vectors of X that have already been orthonormalized.
01184     Teuchos::RCP<const MV> lastQ;
01185     Teuchos::RCP<MV> Xj, MXj;
01186     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01187 
01188     // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
01189     for (int j=0; j<xc; j++) {
01190       
01191       bool dep_flg = false;
01192       
01193       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01194       if (j > 0) {
01195   std::vector<int> index( j );
01196   for (int ind=0; ind<j; ind++) {
01197     index[ind] = ind;
01198   }
01199   lastQ = MVT::CloneView( X, index );
01200 
01201   // Add these views to the Q and C arrays.
01202   Q.push_back( lastQ );
01203   C.push_back( B );
01204   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01205       }
01206       
01207       // Get a view of the current vector in X to orthogonalize.
01208       indX[0] = j;
01209       Xj = MVT::CloneViewNonConst( X, indX );
01210       if (this->_hasOp) {
01211   MXj = MVT::CloneViewNonConst( *MX, indX );
01212       }
01213       else {
01214   MXj = Xj;
01215       }
01216       
01217       // Compute the initial Op-norms
01218       MVT::MvDot( *Xj, *MXj, oldDot );
01219       
01220       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
01221       // Define the product Q^T * (Op*X)
01222       for (int i=0; i<Q.size(); i++) {
01223 
01224   // Get a view of the current serial dense matrix
01225   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01226 
01227   // Multiply Q' with MX
01228   MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,tempC);
01229   // Multiply by Q and subtract the result in Xj
01230   MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
01231   
01232   // Update MXj, with the least number of applications of Op as possible
01233   if (this->_hasOp) {
01234     if (xc <= qcs[i]) {
01235       OPT::Apply( *(this->_Op), *Xj, *MXj);
01236     }
01237     else {
01238       // this will possibly be used again below; don't delete it
01239       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01240       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01241       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01242     }
01243   }
01244       }
01245       
01246       // Compute the Op-norms
01247       MVT::MvDot( *Xj, *MXj, newDot );
01248       
01249       // Do one step of classical Gram-Schmidt orthogonalization 
01250       // with a second correction step if needed.
01251       
01252       if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
01253   
01254   for (int i=0; i<Q.size(); i++) {
01255     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01256     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01257     
01258     // Apply another step of classical Gram-Schmidt
01259     MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,C2);
01260     tempC += C2;
01261     MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
01262     
01263     // Update MXj, with the least number of applications of Op as possible
01264     if (this->_hasOp) {
01265       if (MQ[i].get()) {
01266         // MQ was allocated and computed above; use it
01267         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01268       }
01269       else if (xc <= qcs[i]) {
01270         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01271         OPT::Apply( *(this->_Op), *Xj, *MXj);
01272       }
01273     }
01274   } // for (int i=0; i<Q.size(); i++)
01275   
01276   // Compute the Op-norms after the correction step.
01277   MVT::MvDot( *Xj, *MXj, newDot );
01278   
01279       } // if ()
01280       
01281       // Check for linear dependence.
01282       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01283   dep_flg = true;
01284       }
01285       
01286       // Normalize the new vector if it's not dependent
01287       if (!dep_flg) {
01288   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01289   
01290   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01291   if (this->_hasOp) {
01292     // Update MXj.
01293     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01294   }
01295   
01296   // Enter value on diagonal of B.
01297   (*B)(j,j) = diag;
01298       }
01299       else {
01300   // Create a random vector and orthogonalize it against all previous columns of Q.
01301   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01302   Teuchos::RCP<MV> tempMXj;
01303   MVT::MvRandom( *tempXj );
01304   if (this->_hasOp) {
01305     tempMXj = MVT::Clone( X, 1 );
01306     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01307   } 
01308   else {
01309     tempMXj = tempXj;
01310   }
01311   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01312   //
01313   for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
01314     
01315     for (int i=0; i<Q.size(); i++) {
01316       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01317       
01318       // Apply another step of classical Gram-Schmidt
01319       MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*tempXj,tempMXj,product);
01320       MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
01321       
01322       // Update MXj, with the least number of applications of Op as possible
01323       if (this->_hasOp) {
01324         if (MQ[i].get()) {
01325     // MQ was allocated and computed above; use it
01326     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01327         }
01328         else if (xc <= qcs[i]) {
01329     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01330     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01331         }
01332       }
01333     } // for (int i=0; i<nq; i++)
01334     
01335   }
01336   
01337   // Compute the Op-norms after the correction step.
01338   MVT::MvDot( *tempXj, *tempMXj, newDot );
01339   
01340   // Copy vector into current column of Xj
01341   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01342     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01343     
01344     // Enter value on diagonal of B.
01345     (*B)(j,j) = ZERO;
01346 
01347     // Copy vector into current column of _basisvecs
01348     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01349     if (this->_hasOp) {
01350       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01351     }
01352   }
01353   else {
01354     return j;
01355   } 
01356       } // if (!dep_flg)
01357 
01358       // Remove the vectors from array
01359       if (j > 0) {
01360   Q.resize( nq );
01361   C.resize( nq );
01362   qcs.resize( nq );
01363       }
01364 
01365     } // for (int j=0; j<xc; j++)
01366 
01367     return xc;
01368   }
01369 
01370 } // namespace Belos
01371 
01372 #endif // BELOS_DGKS_ORTHOMANAGER_HPP
01373 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines