Belos Package Browser (Single Doxygen Collection) Development
BelosICGSOrthoManager.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_ICGS_ORTHOMANAGER_HPP
00048 #define BELOS_ICGS_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 namespace Belos {
00065 
00066   template<class ScalarType, class MV, class OP>
00067   class ICGSOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00068 
00069   private:
00070     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00071     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00072     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00073     typedef MultiVecTraits<ScalarType,MV>      MVT;
00074     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00075 
00076   public:
00077     
00079 
00080 
00081     ICGSOrthoManager( const std::string& label = "Belos",
00082                       Teuchos::RCP<const OP> Op = Teuchos::null,
00083           const int max_ortho_steps = 2,
00084           const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00085           const MagnitudeType sing_tol = 10*MGT::eps() )
00086       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00087   max_ortho_steps_( max_ortho_steps ),
00088   blk_tol_( blk_tol ),
00089   sing_tol_( sing_tol ),
00090         label_( label )
00091     {
00092         std::string orthoLabel = label_ + ": Orthogonalization";
00093         timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00094 
00095         std::string updateLabel = label_ + ": Ortho (Update)";
00096         timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00097 
00098         std::string normLabel = label_ + ": Ortho (Norm)";
00099         timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00100 
00101         std::string ipLabel = label_ + ": Ortho (Inner Product)";
00102         timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel); 
00103     };
00104 
00106     ~ICGSOrthoManager() {}
00108 
00109 
00111 
00112 
00114     void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
00115 
00117     void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
00118 
00120     MagnitudeType getBlkTol() const { return blk_tol_; } 
00121 
00123     MagnitudeType getSingTol() const { return sing_tol_; } 
00124 
00126 
00127 
00129 
00130 
00158     void project ( MV &X, Teuchos::RCP<MV> MX, 
00159                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00160                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00161 
00162 
00165     void project ( MV &X, 
00166                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00167                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00168       project(X,Teuchos::null,C,Q);
00169     }
00170 
00171 
00172  
00197     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00198                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00199 
00200 
00203     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00204       return normalize(X,Teuchos::null,B);
00205     }
00206 
00207 
00240     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00241                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00242                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00243                               Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00244 
00247     int projectAndNormalize ( MV &X, 
00248                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00249                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00250                               Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00251       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00252     }
00253 
00255 
00257 
00258 
00262     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00263     orthonormError(const MV &X) const {
00264       return orthonormError(X,Teuchos::null);
00265     }
00266 
00271     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00272     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00273 
00277     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00278     orthogError(const MV &X1, const MV &X2) const {
00279       return orthogError(X1,Teuchos::null,X2);
00280     }
00281 
00286     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00287     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00288 
00290 
00292 
00293 
00296     void setLabel(const std::string& label);
00297 
00300     const std::string& getLabel() const { return label_; }
00301 
00303 
00304   private:
00305     
00307     int max_ortho_steps_;
00308     MagnitudeType blk_tol_;
00309     MagnitudeType sing_tol_;
00310 
00312     std::string label_;
00313     Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, 
00314                                         timerNorm_, timerScale_, timerInnerProd_;
00315   
00317     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00318       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00319       bool completeBasis, int howMany = -1 ) const;
00320     
00322     bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
00323          Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00324          Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00325 
00327     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00328         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00329         Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00330 
00331     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00332            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00333            Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00334            Teuchos::Array<Teuchos::RCP<const MV> > Q) const;    
00335   };
00336 
00338   // Set the label for this orthogonalization manager and create new timers if it's changed
00339   template<class ScalarType, class MV, class OP>
00340   void ICGSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00341   {   
00342     if (label != label_) {
00343       label_ = label;
00344       std::string orthoLabel = label_ + ": Orthogonalization";
00345       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00346 
00347       std::string updateLabel = label_ + ": Ortho (Update)";
00348       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00349 
00350       std::string normLabel = label_ + ": Ortho (Norm)";
00351       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00352 
00353       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00354       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel);
00355     }
00356   } 
00357 
00359   // Compute the distance from orthonormality
00360   template<class ScalarType, class MV, class OP>
00361   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00362   ICGSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00363     const ScalarType ONE = SCT::one();
00364     int rank = MVT::GetNumberVecs(X);
00365     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00366     innerProd(X,X,MX,xTx);
00367     for (int i=0; i<rank; i++) {
00368       xTx(i,i) -= ONE;
00369     }
00370     return xTx.normFrobenius();
00371   }
00372 
00374   // Compute the distance from orthogonality
00375   template<class ScalarType, class MV, class OP>
00376   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00377   ICGSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00378     int r1 = MVT::GetNumberVecs(X1);
00379     int r2  = MVT::GetNumberVecs(X2);
00380     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00381     innerProd(X2,X1,MX1,xTx);
00382     return xTx.normFrobenius();
00383   }
00384 
00386   // Find an Op-orthonormal basis for span(X) - span(W)
00387   template<class ScalarType, class MV, class OP>
00388   int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00389                                     MV &X, Teuchos::RCP<MV> MX, 
00390                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00391                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00392                                     Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00393 
00394     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00395 
00396     ScalarType    ONE  = SCT::one();
00397     ScalarType    ZERO  = SCT::zero();
00398 
00399     int nq = Q.length();
00400     int xc = MVT::GetNumberVecs( X );
00401     int xr = MVT::GetVecLength( X );
00402     int rank = xc;
00403 
00404     /* if the user doesn't want to store the coefficienets, 
00405      * allocate some local memory for them 
00406      */
00407     if ( B == Teuchos::null ) {
00408       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00409     }
00410 
00411     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00412     if (this->_hasOp) {
00413       if (MX == Teuchos::null) {
00414         // we need to allocate space for MX
00415         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00416         OPT::Apply(*(this->_Op),X,*MX);
00417       }
00418     }
00419     else {
00420       // Op == I  -->  MX = X (ignore it if the user passed it in)
00421       MX = Teuchos::rcp( &X, false );
00422     }
00423 
00424     int mxc = MVT::GetNumberVecs( *MX );
00425     int mxr = MVT::GetVecLength( *MX );
00426 
00427     // short-circuit
00428     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::ICGSOrthoManager::projectAndNormalize(): X must be non-empty" );
00429 
00430     int numbas = 0;
00431     for (int i=0; i<nq; i++) {
00432       numbas += MVT::GetNumberVecs( *Q[i] );
00433     }
00434 
00435     // check size of B
00436     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00437                         "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00438     // check size of X and MX
00439     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00440                         "Belos::ICGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00441     // check size of X w.r.t. MX 
00442     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00443                         "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00444     // check feasibility
00445     //TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00446     //                    "Belos::ICGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00447 
00448     // Some flags for checking dependency returns from the internal orthogonalization methods
00449     bool dep_flg = false;
00450 
00451     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00452     Teuchos::RCP<MV> tmpX, tmpMX;
00453     tmpX = MVT::CloneCopy(X);
00454     if (this->_hasOp) {
00455       tmpMX = MVT::CloneCopy(*MX);
00456     }
00457 
00458     if (xc == 1) {
00459 
00460       // Use the cheaper block orthogonalization.
00461       // NOTE: Don't check for dependencies because the update has one std::vector.
00462       dep_flg = blkOrtho1( X, MX, C, Q );
00463 
00464       // Normalize the new block X
00465       {
00466         Teuchos::TimeMonitor normTimer( *timerNorm_ );
00467         if ( B == Teuchos::null ) {
00468           B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00469         }
00470         std::vector<ScalarType> diag(xc);
00471         MVT::MvDot( X, *MX, diag );
00472         (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
00473         rank = 1; 
00474         MVT::MvAddMv( ONE/(*B)(0,0), X, ZERO, X, X );
00475         if (this->_hasOp) {
00476           // Update MXj.
00477     MVT::MvAddMv( ONE/(*B)(0,0), *MX, ZERO, *MX, *MX );
00478         }
00479       }
00480 
00481     } 
00482     else {
00483 
00484       // Use the cheaper block orthogonalization.
00485       dep_flg = blkOrtho( X, MX, C, Q );
00486 
00487       // If a dependency has been detected in this block, then perform
00488       // the more expensive single-std::vector orthogonalization.
00489       if (dep_flg) {
00490         rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00491 
00492         // Copy tmpX back into X.
00493         MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00494         if (this->_hasOp) {
00495     MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00496         }
00497       } 
00498       else {
00499         // There is no dependency, so orthonormalize new block X
00500         rank = findBasis( X, MX, B, false );
00501         if (rank < xc) {
00502     // A dependency was found during orthonormalization of X,
00503     // rerun orthogonalization using more expensive single-std::vector orthogonalization.
00504     rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00505 
00506     // Copy tmpX back into X.
00507     MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00508     if (this->_hasOp) {
00509       MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00510     }
00511         }    
00512       }
00513     } // if (xc == 1) {
00514 
00515     // this should not raise an std::exception; but our post-conditions oblige us to check
00516     TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00517                         "Belos::ICGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00518 
00519     // Return the rank of X.
00520     return rank;
00521   }
00522 
00523 
00524 
00526   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00527   template<class ScalarType, class MV, class OP>
00528   int ICGSOrthoManager<ScalarType, MV, OP>::normalize(
00529                                 MV &X, Teuchos::RCP<MV> MX, 
00530                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00531 
00532     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00533 
00534     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00535     return findBasis(X, MX, B, true);
00536   }
00537 
00538 
00539 
00541   template<class ScalarType, class MV, class OP>
00542   void ICGSOrthoManager<ScalarType, MV, OP>::project(
00543                           MV &X, Teuchos::RCP<MV> MX, 
00544                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00545                           Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00546     // For the inner product defined by the operator Op or the identity (Op == 0)
00547     //   -> Orthogonalize X against each Q[i]
00548     // Modify MX accordingly
00549     //
00550     // Note that when Op is 0, MX is not referenced
00551     //
00552     // Parameter variables
00553     //
00554     // X  : Vectors to be transformed
00555     //
00556     // MX : Image of the block std::vector X by the mass matrix
00557     //
00558     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00559     //
00560     
00561     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00562     
00563     int xc = MVT::GetNumberVecs( X );
00564     int xr = MVT::GetVecLength( X );
00565     int nq = Q.length();
00566     std::vector<int> qcs(nq);
00567     // short-circuit
00568     if (nq == 0 || xc == 0 || xr == 0) {
00569       return;
00570     }
00571     int qr = MVT::GetVecLength ( *Q[0] );
00572     // if we don't have enough C, expand it with null references
00573     // if we have too many, resize to throw away the latter ones
00574     // if we have exactly as many as we have Q, this call has no effect
00575     C.resize(nq);
00576 
00577 
00578     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00579     if (this->_hasOp) {
00580       if (MX == Teuchos::null) {
00581         // we need to allocate space for MX
00582         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00583         OPT::Apply(*(this->_Op),X,*MX);
00584       }
00585     }
00586     else {
00587       // Op == I  -->  MX = X (ignore it if the user passed it in)
00588       MX = Teuchos::rcp( &X, false );
00589     }
00590     int mxc = MVT::GetNumberVecs( *MX );
00591     int mxr = MVT::GetVecLength( *MX );
00592 
00593     // check size of X and Q w.r.t. common sense
00594     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00595                         "Belos::ICGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00596     // check size of X w.r.t. MX and Q
00597     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00598                         "Belos::ICGSOrthoManager::project(): Size of X not consistant with MX,Q" );
00599 
00600     // tally up size of all Q and check/allocate C
00601     int baslen = 0;
00602     for (int i=0; i<nq; i++) {
00603       TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00604                           "Belos::ICGSOrthoManager::project(): Q lengths not mutually consistant" );
00605       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00606       TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00607                           "Belos::ICGSOrthoManager::project(): Q has less rows than columns" );
00608       baslen += qcs[i];
00609 
00610       // check size of C[i]
00611       if ( C[i] == Teuchos::null ) {
00612         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00613       }
00614       else {
00615         TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00616                            "Belos::ICGSOrthoManager::project(): Size of Q not consistant with size of C" );
00617       }
00618     }
00619 
00620     // Use the cheaper block orthogonalization, don't check for rank deficiency.
00621     blkOrtho( X, MX, C, Q );
00622 
00623   }  
00624 
00626   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00627   // the rank is numvectors(X)
00628   template<class ScalarType, class MV, class OP>
00629   int ICGSOrthoManager<ScalarType, MV, OP>::findBasis(
00630                   MV &X, Teuchos::RCP<MV> MX, 
00631                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00632                   bool completeBasis, int howMany ) const {
00633     // For the inner product defined by the operator Op or the identity (Op == 0)
00634     //   -> Orthonormalize X 
00635     // Modify MX accordingly
00636     //
00637     // Note that when Op is 0, MX is not referenced
00638     //
00639     // Parameter variables
00640     //
00641     // X  : Vectors to be orthonormalized
00642     //
00643     // MX : Image of the multivector X under the operator Op
00644     //
00645     // Op  : Pointer to the operator for the inner product
00646     //
00647     //
00648 
00649     Teuchos::TimeMonitor normTimer( *timerNorm_ );
00650 
00651     const ScalarType ONE  = SCT::one();
00652     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00653 
00654     int xc = MVT::GetNumberVecs( X );
00655     int xr = MVT::GetVecLength( X );
00656 
00657     if (howMany == -1) {
00658       howMany = xc;
00659     }
00660 
00661     /*******************************************************
00662      *  If _hasOp == false, we will not reference MX below *
00663      *******************************************************/
00664 
00665     // if Op==null, MX == X (via pointer)
00666     // Otherwise, either the user passed in MX or we will allocated and compute it
00667     if (this->_hasOp) {
00668       if (MX == Teuchos::null) {
00669         // we need to allocate space for MX
00670         MX = MVT::Clone(X,xc);
00671         OPT::Apply(*(this->_Op),X,*MX);
00672       }
00673     }
00674 
00675     /* if the user doesn't want to store the coefficienets, 
00676      * allocate some local memory for them 
00677      */
00678     if ( B == Teuchos::null ) {
00679       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00680     }
00681 
00682     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00683     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00684 
00685     // check size of C, B
00686     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00687                         "Belos::ICGSOrthoManager::findBasis(): X must be non-empty" );
00688     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00689                         "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
00690     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00691                         "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
00692     TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00693                         "Belos::ICGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
00694     TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
00695                         "Belos::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
00696 
00697     /* xstart is which column we are starting the process with, based on howMany
00698      * columns before xstart are assumed to be Op-orthonormal already
00699      */
00700     int xstart = xc - howMany;
00701 
00702     for (int j = xstart; j < xc; j++) {
00703 
00704       // numX is 
00705       // * number of currently orthonormal columns of X
00706       // * the index of the current column of X
00707       int numX = j;
00708       bool addVec = false;
00709 
00710       // Get a view of the std::vector currently being worked on.
00711       std::vector<int> index(1);
00712       index[0] = numX;
00713       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
00714       Teuchos::RCP<MV> MXj;
00715       if ((this->_hasOp)) {
00716         // MXj is a view of the current std::vector in MX
00717         MXj = MVT::CloneViewNonConst( *MX, index );
00718       }
00719       else {
00720         // MXj is a pointer to Xj, and MUST NOT be modified
00721         MXj = Xj;
00722       }
00723 
00724       // Get a view of the previous vectors.
00725       std::vector<int> prev_idx( numX );
00726       Teuchos::RCP<const MV> prevX, prevMX;
00727 
00728       if (numX > 0) {
00729         for (int i=0; i<numX; i++) {
00730           prev_idx[i] = i;
00731         }
00732         prevX = MVT::CloneView( X, prev_idx );
00733         if (this->_hasOp) {
00734           prevMX = MVT::CloneView( *MX, prev_idx );
00735         }
00736       } 
00737 
00738       // Make storage for these Gram-Schmidt iterations.
00739       Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
00740       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00741       //
00742       // Save old MXj std::vector and compute Op-norm
00743       //
00744       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00745       MVT::MvDot( *Xj, *MXj, oldDot );
00746       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
00747       TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
00748         "Belos::ICGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
00749 
00750       if (numX > 0) {
00751   
00752         Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
00753   
00754   for (int i=0; i<max_ortho_steps_; ++i) {
00755     
00756     // product <- prevX^T MXj
00757           {
00758             Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00759       innerProd(*prevX,*Xj,MXj,P2);
00760           }
00761     
00762     // Xj <- Xj - prevX prevX^T MXj   
00763     //     = Xj - prevX product
00764           {
00765             Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00766       MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00767           }
00768     
00769     // Update MXj
00770     if (this->_hasOp) {
00771       // MXj <- Op*Xj_new
00772       //      = Op*(Xj_old - prevX prevX^T MXj)
00773       //      = MXj - prevMX product
00774             Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00775       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00776     }
00777     
00778     // Set coefficients
00779     if ( i==0 )
00780       product = P2;
00781     else
00782       product += P2;
00783         }
00784   
00785       } // if (numX > 0)
00786 
00787       // Compute Op-norm with old MXj
00788       MVT::MvDot( *Xj, *oldMXj, newDot );
00789 
00790       // Check to see if the new std::vector is dependent.
00791       if (completeBasis) {
00792   //
00793   // We need a complete basis, so add random vectors if necessary
00794   //
00795   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
00796     
00797     // Add a random std::vector and orthogonalize it against previous vectors in block.
00798     addVec = true;
00799 #ifdef ORTHO_DEBUG
00800     std::cout << "Belos::ICGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
00801 #endif
00802     //
00803     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
00804     Teuchos::RCP<MV> tempMXj;
00805     MVT::MvRandom( *tempXj );
00806     if (this->_hasOp) {
00807       tempMXj = MVT::Clone( X, 1 );
00808       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
00809     } 
00810     else {
00811       tempMXj = tempXj;
00812     }
00813     MVT::MvDot( *tempXj, *tempMXj, oldDot );
00814     //
00815     for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
00816             {
00817               Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00818         innerProd(*prevX,*tempXj,tempMXj,product);
00819             }
00820             {
00821               Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00822         MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
00823             }
00824       if (this->_hasOp) {
00825               Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00826         MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
00827       }
00828     }
00829     // Compute new Op-norm
00830     MVT::MvDot( *tempXj, *tempMXj, newDot );
00831     //
00832     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
00833       // Copy std::vector into current column of _basisvecs
00834       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
00835       if (this->_hasOp) {
00836         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
00837       }
00838     }
00839     else {
00840       return numX;
00841     } 
00842   }
00843       }
00844       else {
00845   //
00846   // We only need to detect dependencies.
00847   //
00848   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
00849     return numX;
00850   }
00851       }
00852       
00853       // If we haven't left this method yet, then we can normalize the new std::vector Xj.
00854       // Normalize Xj.
00855       // Xj <- Xj / std::sqrt(newDot)
00856       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
00857       {
00858         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
00859         if (this->_hasOp) {
00860     // Update MXj.
00861     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
00862         }
00863       }
00864 
00865       // If we've added a random std::vector, enter a zero in the j'th diagonal element.
00866       if (addVec) {
00867   (*B)(j,j) = ZERO;
00868       }
00869       else {
00870   (*B)(j,j) = diag;
00871       }
00872 
00873       // Save the coefficients, if we are working on the original std::vector and not a randomly generated one
00874       if (!addVec) {
00875   for (int i=0; i<numX; i++) {
00876     (*B)(i,j) = product(i,0);
00877   }
00878       }
00879 
00880     } // for (j = 0; j < xc; ++j)
00881 
00882     return xc;
00883   }
00884 
00886   // Routine to compute the block orthogonalization
00887   template<class ScalarType, class MV, class OP>
00888   bool 
00889   ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
00890                 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00891                 Teuchos::Array<Teuchos::RCP<const MV> > Q) const
00892   {
00893     int nq = Q.length();
00894     int xc = MVT::GetNumberVecs( X );
00895     const ScalarType ONE  = SCT::one();
00896 
00897     std::vector<int> qcs( nq );
00898     for (int i=0; i<nq; i++) {
00899       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00900     }
00901 
00902     // Perform the Gram-Schmidt transformation for a block of vectors
00903 
00904     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
00905     // Define the product Q^T * (Op*X)
00906     for (int i=0; i<nq; i++) {
00907       // Multiply Q' with MX
00908       {
00909         Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00910         innerProd(*Q[i],X,MX,*C[i]);
00911       }
00912       // Multiply by Q and subtract the result in X
00913       {
00914         Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00915         MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
00916       }
00917 
00918       // Update MX, with the least number of applications of Op as possible
00919       if (this->_hasOp) {
00920         if (xc <= qcs[i]) {
00921           OPT::Apply( *(this->_Op), X, *MX);
00922         }
00923         else {
00924           // this will possibly be used again below; don't delete it
00925           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
00926           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
00927           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
00928         }
00929       }
00930     }
00931 
00932     // Do as many steps of classical Gram-Schmidt as required by max_ortho_steps_
00933     for (int j = 1; j < max_ortho_steps_; ++j) {
00934       
00935       for (int i=0; i<nq; i++) {
00936   Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
00937         
00938   // Apply another step of classical Gram-Schmidt
00939   {
00940           Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00941           innerProd(*Q[i],X,MX,C2);
00942         }
00943   *C[i] += C2;
00944         {
00945           Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00946     MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
00947         }
00948         
00949   // Update MX, with the least number of applications of Op as possible
00950   if (this->_hasOp) {
00951     if (MQ[i].get()) {
00952       // MQ was allocated and computed above; use it
00953       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
00954     }
00955     else if (xc <= qcs[i]) {
00956       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
00957       OPT::Apply( *(this->_Op), X, *MX);
00958     }
00959   }
00960       } // for (int i=0; i<nq; i++)
00961     } // for (int j = 0; j < max_ortho_steps; ++j)
00962   
00963     return false;
00964   }
00965 
00967   // Routine to compute the block orthogonalization
00968   template<class ScalarType, class MV, class OP>
00969   bool 
00970   ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00971                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00972                Teuchos::Array<Teuchos::RCP<const MV> > Q) const
00973   {
00974     int nq = Q.length();
00975     int xc = MVT::GetNumberVecs( X );
00976     bool dep_flg = false;
00977     const ScalarType ONE  = SCT::one();
00978 
00979     std::vector<int> qcs( nq );
00980     for (int i=0; i<nq; i++) {
00981       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00982     }
00983 
00984     // Perform the Gram-Schmidt transformation for a block of vectors
00985     
00986     // Compute the initial Op-norms
00987     std::vector<ScalarType> oldDot( xc );
00988     MVT::MvDot( X, *MX, oldDot );
00989 
00990     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
00991     // Define the product Q^T * (Op*X)
00992     for (int i=0; i<nq; i++) {
00993       // Multiply Q' with MX
00994       {
00995         Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00996         innerProd(*Q[i],X,MX,*C[i]);
00997       }
00998       // Multiply by Q and subtract the result in X
00999       {
01000         Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01001         MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
01002       }
01003       // Update MX, with the least number of applications of Op as possible
01004       if (this->_hasOp) {
01005         if (xc <= qcs[i]) {
01006           OPT::Apply( *(this->_Op), X, *MX);
01007         }
01008         else {
01009           // this will possibly be used again below; don't delete it
01010           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01011           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01012           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01013         }
01014       }
01015     }
01016 
01017     // Do as many steps of classical Gram-Schmidt as required by max_ortho_steps_
01018     for (int j = 1; j < max_ortho_steps_; ++j) {
01019       
01020       for (int i=0; i<nq; i++) {
01021   Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
01022         
01023   // Apply another step of classical Gram-Schmidt
01024         {
01025           Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01026     innerProd(*Q[i],X,MX,C2);
01027         }
01028   *C[i] += C2;
01029         {
01030           Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01031     MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
01032         }
01033         
01034   // Update MX, with the least number of applications of Op as possible
01035   if (this->_hasOp) {
01036     if (MQ[i].get()) {
01037       // MQ was allocated and computed above; use it
01038       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01039     }
01040     else if (xc <= qcs[i]) {
01041       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01042       OPT::Apply( *(this->_Op), X, *MX);
01043     }
01044   }
01045       } // for (int i=0; i<nq; i++)
01046     } // for (int j = 0; j < max_ortho_steps; ++j)
01047   
01048     // Compute new Op-norms
01049     std::vector<ScalarType> newDot(xc);
01050     MVT::MvDot( X, *MX, newDot );
01051  
01052     // Check to make sure the new block of vectors are not dependent on previous vectors
01053     for (int i=0; i<xc; i++){
01054       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01055   dep_flg = true;
01056   break;
01057       }
01058     } // end for (i=0;...)
01059 
01060     return dep_flg;
01061   }
01062   
01064   // Routine to compute the block orthogonalization using single-std::vector orthogonalization
01065   template<class ScalarType, class MV, class OP>
01066   int
01067   ICGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01068                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01069                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01070                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const
01071   {
01072     const ScalarType ONE  = SCT::one();
01073     const ScalarType ZERO  = SCT::zero();
01074     
01075     int nq = Q.length();
01076     int xc = MVT::GetNumberVecs( X );
01077     std::vector<int> indX( 1 );
01078     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01079 
01080     std::vector<int> qcs( nq );
01081     for (int i=0; i<nq; i++) {
01082       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01083     }
01084 
01085     // Create pointers for the previous vectors of X that have already been orthonormalized.
01086     Teuchos::RCP<const MV> lastQ;
01087     Teuchos::RCP<MV> Xj, MXj;
01088     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01089 
01090     // Perform the Gram-Schmidt transformation for each std::vector in the block of vectors.
01091     for (int j=0; j<xc; j++) {
01092       
01093       bool dep_flg = false;
01094       
01095       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01096       if (j > 0) {
01097   std::vector<int> index( j );
01098   for (int ind=0; ind<j; ind++) {
01099     index[ind] = ind;
01100   }
01101   lastQ = MVT::CloneView( X, index );
01102 
01103   // Add these views to the Q and C arrays.
01104   Q.push_back( lastQ );
01105   C.push_back( B );
01106   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01107       }
01108       
01109       // Get a view of the current std::vector in X to orthogonalize.
01110       indX[0] = j;
01111       Xj = MVT::CloneViewNonConst( X, indX );
01112       if (this->_hasOp) {
01113   MXj = MVT::CloneViewNonConst( *MX, indX );
01114       }
01115       else {
01116   MXj = Xj;
01117       }
01118       
01119       // Compute the initial Op-norms
01120       MVT::MvDot( *Xj, *MXj, oldDot );
01121       
01122       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.length());
01123       // Define the product Q^T * (Op*X)
01124       for (int i=0; i<Q.length(); i++) {
01125 
01126   // Get a view of the current serial dense matrix
01127   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01128 
01129   // Multiply Q' with MX
01130   innerProd(*Q[i],*Xj,MXj,tempC);
01131   // Multiply by Q and subtract the result in Xj
01132   MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
01133   
01134   // Update MXj, with the least number of applications of Op as possible
01135   if (this->_hasOp) {
01136     if (xc <= qcs[i]) {
01137       OPT::Apply( *(this->_Op), *Xj, *MXj);
01138     }
01139     else {
01140       // this will possibly be used again below; don't delete it
01141       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01142       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01143       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01144     }
01145   }
01146       }
01147      
01148       // Do any additional steps of classical Gram-Schmidt orthogonalization 
01149       for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
01150   
01151   for (int i=0; i<Q.length(); i++) {
01152     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01153     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01154     
01155     // Apply another step of classical Gram-Schmidt
01156     innerProd(*Q[i],*Xj,MXj,C2);
01157     tempC += C2;
01158     MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
01159     
01160     // Update MXj, with the least number of applications of Op as possible
01161     if (this->_hasOp) {
01162       if (MQ[i].get()) {
01163         // MQ was allocated and computed above; use it
01164         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01165       }
01166       else if (xc <= qcs[i]) {
01167         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01168         OPT::Apply( *(this->_Op), *Xj, *MXj);
01169       }
01170     }
01171   } // for (int i=0; i<Q.length(); i++)
01172   
01173       } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
01174       
01175       // Check for linear dependence.
01176       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01177   dep_flg = true;
01178       }
01179       
01180       // Normalize the new std::vector if it's not dependent
01181       if (!dep_flg) {
01182   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01183   
01184   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01185   if (this->_hasOp) {
01186     // Update MXj.
01187     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01188   }
01189   
01190   // Enter value on diagonal of B.
01191   (*B)(j,j) = diag;
01192       }
01193       else {
01194   // Create a random std::vector and orthogonalize it against all previous columns of Q.
01195   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01196   Teuchos::RCP<MV> tempMXj;
01197   MVT::MvRandom( *tempXj );
01198   if (this->_hasOp) {
01199     tempMXj = MVT::Clone( X, 1 );
01200     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01201   } 
01202   else {
01203     tempMXj = tempXj;
01204   }
01205   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01206   //
01207   for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
01208     
01209     for (int i=0; i<Q.length(); i++) {
01210       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01211       
01212       // Apply another step of classical Gram-Schmidt
01213       innerProd(*Q[i],*tempXj,tempMXj,product);
01214       MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
01215       
01216       // Update MXj, with the least number of applications of Op as possible
01217       if (this->_hasOp) {
01218         if (MQ[i].get()) {
01219     // MQ was allocated and computed above; use it
01220     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01221         }
01222         else if (xc <= qcs[i]) {
01223     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01224     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01225         }
01226       }
01227     } // for (int i=0; i<nq; i++)
01228     
01229   }
01230   
01231   // Compute the Op-norms after the correction step.
01232   MVT::MvDot( *tempXj, *tempMXj, newDot );
01233   
01234   // Copy std::vector into current column of Xj
01235   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01236     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01237     
01238     // Enter value on diagonal of B.
01239     (*B)(j,j) = ZERO;
01240 
01241     // Copy std::vector into current column of _basisvecs
01242     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01243     if (this->_hasOp) {
01244       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01245     }
01246   }
01247   else {
01248     return j;
01249   } 
01250       } // if (!dep_flg)
01251 
01252       // Remove the vectors from array
01253       if (j > 0) {
01254   Q.resize( nq );
01255   C.resize( nq );
01256   qcs.resize( nq );
01257       }
01258 
01259     } // for (int j=0; j<xc; j++)
01260 
01261     return xc;
01262   }
01263 
01264 } // namespace Belos
01265 
01266 #endif // BELOS_ICGS_ORTHOMANAGER_HPP
01267 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines