Belos Package Browser (Single Doxygen Collection) Development
BelosIMGSOrthoManager.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_IMGS_ORTHOMANAGER_HPP
00048 #define BELOS_IMGS_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 IMGSOrthoManager : 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     IMGSOrthoManager( const std::string& label = "Belos",
00082                       Teuchos::RCP<const OP> Op = Teuchos::null,
00083           const int max_ortho_steps = 1,
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     ~IMGSOrthoManager() {}
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 IMGSOrthoManager<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   IMGSOrthoManager<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   IMGSOrthoManager<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 IMGSOrthoManager<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::IMGSOrthoManager::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::IMGSOrthoManager::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::IMGSOrthoManager::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::IMGSOrthoManager::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::IMGSOrthoManager::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::IMGSOrthoManager::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 IMGSOrthoManager<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 IMGSOrthoManager<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::IMGSOrthoManager::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::IMGSOrthoManager::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::IMGSOrthoManager::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::IMGSOrthoManager::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::IMGSOrthoManager::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 IMGSOrthoManager<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::IMGSOrthoManager::findBasis(): X must be non-empty" );
00688     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00689                         "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
00690     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00691                         "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
00692     TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00693                         "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
00694     TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
00695                         "Belos::IMGSOrthoManager::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       // Make storage for these Gram-Schmidt iterations.
00725       Teuchos::SerialDenseVector<int,ScalarType> product(numX);
00726       Teuchos::SerialDenseVector<int,ScalarType> P2(1);
00727       Teuchos::RCP<const MV> prevX, prevMX;
00728   
00729       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00730       //
00731       // Save old MXj std::vector and compute Op-norm
00732       //
00733       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00734       MVT::MvDot( *Xj, *MXj, oldDot );
00735       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
00736       TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
00737         "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
00738 
00739       // Perform MGS one vector at a time
00740       for (int ii=0; ii<numX; ii++) {
00741   
00742   index[0] = ii;
00743   prevX = MVT::CloneView( X, index );
00744   if (this->_hasOp) {
00745     prevMX = MVT::CloneView( *MX, index );
00746   }
00747   
00748   for (int i=0; i<max_ortho_steps_; ++i) {
00749     
00750     // product <- prevX^T MXj
00751     {
00752       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00753       innerProd(*prevX,*Xj,MXj,P2);
00754     }
00755     
00756     // Xj <- Xj - prevX prevX^T MXj   
00757     //     = Xj - prevX product
00758     {
00759       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00760       MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00761     }
00762     
00763     // Update MXj
00764     if (this->_hasOp) {
00765       // MXj <- Op*Xj_new
00766       //      = Op*(Xj_old - prevX prevX^T MXj)
00767       //      = MXj - prevMX product
00768       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00769       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00770     }
00771     
00772     // Set coefficients
00773     if ( i==0 )
00774       product[ii] = P2[0];
00775     else
00776       product[ii] += P2[0];
00777     
00778   } // for (int i=0; i<max_ortho_steps_; ++i)
00779 
00780       } // for (int ii=0; ii<numX; ++ii)  
00781   
00782   // Compute Op-norm with old MXj
00783       MVT::MvDot( *Xj, *oldMXj, newDot );
00784       
00785       // Check to see if the new std::vector is dependent.
00786       if (completeBasis) {
00787   //
00788   // We need a complete basis, so add random vectors if necessary
00789   //
00790   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
00791     
00792     // Add a random std::vector and orthogonalize it against previous vectors in block.
00793     addVec = true;
00794 #ifdef ORTHO_DEBUG
00795     std::cout << "Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
00796 #endif
00797     //
00798     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
00799     Teuchos::RCP<MV> tempMXj;
00800     MVT::MvRandom( *tempXj );
00801     if (this->_hasOp) {
00802       tempMXj = MVT::Clone( X, 1 );
00803       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
00804     } 
00805     else {
00806       tempMXj = tempXj;
00807     }
00808     MVT::MvDot( *tempXj, *tempMXj, oldDot );
00809     //
00810     // Perform MGS one vector at a time
00811     for (int ii=0; ii<numX; ii++) {
00812       
00813       index[0] = ii;
00814       prevX = MVT::CloneView( X, index );
00815       if (this->_hasOp) {
00816         prevMX = MVT::CloneView( *MX, index );
00817       }
00818       
00819       for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
00820         {
00821     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00822     innerProd(*prevX,*tempXj,tempMXj,P2);
00823         }
00824         {
00825     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00826     MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
00827         }
00828         if (this->_hasOp) {
00829     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00830     MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
00831         }
00832 
00833         // Set coefficients
00834         if ( num_orth==0 )
00835     product[ii] = P2[0];
00836         else
00837     product[ii] += P2[0];       
00838       }
00839     }
00840       
00841     // Compute new Op-norm
00842     MVT::MvDot( *tempXj, *tempMXj, newDot );
00843     //
00844     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
00845       // Copy std::vector into current column of _basisvecs
00846       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
00847       if (this->_hasOp) {
00848         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
00849       }
00850     }
00851     else {
00852       return numX;
00853     } 
00854   }
00855       }
00856       else {
00857   //
00858   // We only need to detect dependencies.
00859   //
00860   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
00861     return numX;
00862   }
00863       }
00864       
00865       
00866       // If we haven't left this method yet, then we can normalize the new vector Xj.
00867       // Normalize Xj.
00868       // Xj <- Xj / std::sqrt(newDot)
00869       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
00870       {
00871         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
00872         if (this->_hasOp) {
00873     // Update MXj.
00874     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
00875         }
00876       }
00877       
00878       // If we've added a random std::vector, enter a zero in the j'th diagonal element.
00879       if (addVec) {
00880   (*B)(j,j) = ZERO;
00881       }
00882       else {
00883   (*B)(j,j) = diag;
00884       }
00885       
00886       // Save the coefficients, if we are working on the original std::vector and not a randomly generated one
00887       if (!addVec) {
00888   for (int i=0; i<numX; i++) {
00889     (*B)(i,j) = product(i);
00890   }
00891       }
00892       
00893     } // for (j = 0; j < xc; ++j)
00894     
00895     return xc;
00896   }
00897   
00899   // Routine to compute the block orthogonalization
00900   template<class ScalarType, class MV, class OP>
00901   bool 
00902   IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
00903                 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00904                 Teuchos::Array<Teuchos::RCP<const MV> > Q) const
00905   {
00906     int nq = Q.length();
00907     int xc = MVT::GetNumberVecs( X );
00908     const ScalarType ONE  = SCT::one();
00909 
00910     std::vector<int> qcs( nq );
00911     for (int i=0; i<nq; i++) {
00912       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00913     }
00914 
00915     // Perform the Gram-Schmidt transformation for a block of vectors
00916     std::vector<int> index(1);
00917     Teuchos::RCP<const MV> tempQ;
00918 
00919     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
00920     // Define the product Q^T * (Op*X)
00921     for (int i=0; i<nq; i++) {
00922 
00923       // Perform MGS one vector at a time
00924       for (int ii=0; ii<qcs[i]; ii++) {
00925   
00926   index[0] = ii;
00927   tempQ = MVT::CloneView( *Q[i], index );
00928   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
00929 
00930   // Multiply Q' with MX
00931   {
00932     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00933     innerProd(*tempQ,X,MX,tempC);
00934   }
00935   // Multiply by Q and subtract the result in X
00936   {
00937     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00938     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
00939   }
00940       }
00941       
00942       // Update MX, with the least number of applications of Op as possible
00943       if (this->_hasOp) {
00944   if (xc <= qcs[i]) {
00945     OPT::Apply( *(this->_Op), X, *MX);
00946   }
00947         else {
00948           // this will possibly be used again below; don't delete it
00949           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
00950           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
00951           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
00952         }
00953       }
00954     }
00955     
00956     // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
00957     for (int j = 1; j < max_ortho_steps_; ++j) {
00958       
00959       for (int i=0; i<nq; i++) {
00960 
00961   Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
00962         
00963   // Perform MGS one vector at a time
00964   for (int ii=0; ii<qcs[i]; ii++) {
00965     
00966     index[0] = ii;
00967     tempQ = MVT::CloneView( *Q[i], index );
00968     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
00969     Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
00970     
00971     // Apply another step of modified Gram-Schmidt
00972     {
00973       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00974       innerProd( *tempQ, X, MX, tempC2 );
00975     }
00976     tempC += tempC2;
00977     {
00978       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00979       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
00980     }
00981     
00982   }
00983 
00984   // Update MX, with the least number of applications of Op as possible
00985   if (this->_hasOp) {
00986     if (MQ[i].get()) {
00987       // MQ was allocated and computed above; use it
00988       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
00989     }
00990     else if (xc <= qcs[i]) {
00991       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
00992       OPT::Apply( *(this->_Op), X, *MX);
00993     }
00994   }
00995       } // for (int i=0; i<nq; i++)
00996     } // for (int j = 0; j < max_ortho_steps; ++j)
00997   
00998     return false;
00999   }
01000 
01002   // Routine to compute the block orthogonalization
01003   template<class ScalarType, class MV, class OP>
01004   bool 
01005   IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
01006                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01007                Teuchos::Array<Teuchos::RCP<const MV> > Q) const
01008   {
01009     int nq = Q.length();
01010     int xc = MVT::GetNumberVecs( X );
01011     bool dep_flg = false;
01012     const ScalarType ONE  = SCT::one();
01013 
01014     std::vector<int> qcs( nq );
01015     for (int i=0; i<nq; i++) {
01016       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01017     }
01018 
01019     // Perform the Gram-Schmidt transformation for a block of vectors
01020     
01021     // Compute the initial Op-norms
01022     std::vector<ScalarType> oldDot( xc );
01023     MVT::MvDot( X, *MX, oldDot );
01024 
01025     std::vector<int> index(1);
01026     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01027     Teuchos::RCP<const MV> tempQ;
01028 
01029     // Define the product Q^T * (Op*X)
01030     for (int i=0; i<nq; i++) {
01031 
01032       // Perform MGS one vector at a time
01033       for (int ii=0; ii<qcs[i]; ii++) {
01034   
01035   index[0] = ii;
01036   tempQ = MVT::CloneView( *Q[i], index );
01037   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
01038 
01039   // Multiply Q' with MX
01040   {
01041     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01042     innerProd( *tempQ, X, MX, tempC);
01043   }
01044   // Multiply by Q and subtract the result in X
01045   {
01046     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01047     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
01048   }
01049       }
01050 
01051       // Update MX, with the least number of applications of Op as possible
01052       if (this->_hasOp) {
01053         if (xc <= qcs[i]) {
01054           OPT::Apply( *(this->_Op), X, *MX);
01055         }
01056         else {
01057           // this will possibly be used again below; don't delete it
01058           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01059           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01060           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01061         }
01062       }
01063     }
01064 
01065     // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
01066     for (int j = 1; j < max_ortho_steps_; ++j) {
01067       
01068       for (int i=0; i<nq; i++) {
01069   Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
01070         
01071   // Perform MGS one vector at a time
01072   for (int ii=0; ii<qcs[i]; ii++) {
01073     
01074     index[0] = ii;
01075     tempQ = MVT::CloneView( *Q[i], index );
01076     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
01077     Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
01078     
01079     // Apply another step of modified Gram-Schmidt
01080     {
01081       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01082       innerProd( *tempQ, X, MX, tempC2 );
01083     }
01084     tempC += tempC2;
01085     {
01086       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01087       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
01088     }
01089   }
01090 
01091   // Update MX, with the least number of applications of Op as possible
01092   if (this->_hasOp) {
01093     if (MQ[i].get()) {
01094       // MQ was allocated and computed above; use it
01095       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01096     }
01097     else if (xc <= qcs[i]) {
01098       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01099       OPT::Apply( *(this->_Op), X, *MX);
01100     }
01101   }
01102       } // for (int i=0; i<nq; i++)
01103     } // for (int j = 0; j < max_ortho_steps; ++j)
01104   
01105     // Compute new Op-norms
01106     std::vector<ScalarType> newDot(xc);
01107     MVT::MvDot( X, *MX, newDot );
01108     
01109     // Check to make sure the new block of vectors are not dependent on previous vectors
01110     for (int i=0; i<xc; i++){
01111       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01112   dep_flg = true;
01113   break;
01114       }
01115     } // end for (i=0;...)
01116     
01117     return dep_flg;
01118   }
01119   
01121   // Routine to compute the block orthogonalization using single-std::vector orthogonalization
01122   template<class ScalarType, class MV, class OP>
01123   int
01124   IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01125                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01126                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01127                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const
01128   {
01129     const ScalarType ONE  = SCT::one();
01130     const ScalarType ZERO  = SCT::zero();
01131     
01132     int nq = Q.length();
01133     int xc = MVT::GetNumberVecs( X );
01134     std::vector<int> indX( 1 );
01135     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01136 
01137     std::vector<int> qcs( nq );
01138     for (int i=0; i<nq; i++) {
01139       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01140     }
01141 
01142     // Create pointers for the previous vectors of X that have already been orthonormalized.
01143     Teuchos::RCP<const MV> lastQ;
01144     Teuchos::RCP<MV> Xj, MXj;
01145     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01146 
01147     // Perform the Gram-Schmidt transformation for each std::vector in the block of vectors.
01148     for (int j=0; j<xc; j++) {
01149       
01150       bool dep_flg = false;
01151       
01152       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01153       if (j > 0) {
01154   std::vector<int> index( j );
01155   for (int ind=0; ind<j; ind++) {
01156     index[ind] = ind;
01157   }
01158   lastQ = MVT::CloneView( X, index );
01159 
01160   // Add these views to the Q and C arrays.
01161   Q.push_back( lastQ );
01162   C.push_back( B );
01163   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01164       }
01165       
01166       // Get a view of the current std::vector in X to orthogonalize.
01167       indX[0] = j;
01168       Xj = MVT::CloneViewNonConst( X, indX );
01169       if (this->_hasOp) {
01170   MXj = MVT::CloneViewNonConst( *MX, indX );
01171       }
01172       else {
01173   MXj = Xj;
01174       }
01175       
01176       // Compute the initial Op-norms
01177       MVT::MvDot( *Xj, *MXj, oldDot );
01178       
01179       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.length());
01180       Teuchos::RCP<const MV> tempQ;
01181 
01182       // Define the product Q^T * (Op*X)
01183       for (int i=0; i<Q.length(); i++) {
01184 
01185   // Perform MGS one vector at a time
01186   for (int ii=0; ii<qcs[i]; ii++) {
01187     
01188     indX[0] = ii;
01189     tempQ = MVT::CloneView( *Q[i], indX );
01190     // Get a view of the current serial dense matrix
01191     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
01192     
01193     // Multiply Q' with MX
01194     innerProd(*tempQ,*Xj,MXj,tempC);
01195 
01196     // Multiply by Q and subtract the result in Xj
01197     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
01198   }
01199  
01200   // Update MXj, with the least number of applications of Op as possible
01201   if (this->_hasOp) {
01202     if (xc <= qcs[i]) {
01203       OPT::Apply( *(this->_Op), *Xj, *MXj);
01204     }
01205     else {
01206       // this will possibly be used again below; don't delete it
01207       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01208       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01209       Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01210       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01211     }
01212   }
01213       }
01214       
01215       // Do any additional steps of modified Gram-Schmidt orthogonalization 
01216       for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
01217   
01218   for (int i=0; i<Q.length(); i++) {
01219     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01220     
01221     // Perform MGS one vector at a time
01222     for (int ii=0; ii<qcs[i]; ii++) {
01223       
01224       indX[0] = ii;
01225       tempQ = MVT::CloneView( *Q[i], indX );
01226       // Get a view of the current serial dense matrix
01227       Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
01228       
01229       // Apply another step of modified Gram-Schmidt
01230       innerProd( *tempQ, *Xj, MXj, tempC2);
01231       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj ); 
01232     }
01233     
01234     // Add the coefficients into C[i]
01235     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01236     tempC += C2;
01237     
01238     // Update MXj, with the least number of applications of Op as possible
01239     if (this->_hasOp) {
01240       if (MQ[i].get()) {
01241         // MQ was allocated and computed above; use it
01242         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01243       }
01244       else if (xc <= qcs[i]) {
01245         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01246         OPT::Apply( *(this->_Op), *Xj, *MXj);
01247       }
01248     }
01249   } // for (int i=0; i<Q.length(); i++)
01250   
01251       } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
01252       
01253       // Check for linear dependence.
01254       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01255   dep_flg = true;
01256       }
01257       
01258       // Normalize the new std::vector if it's not dependent
01259       if (!dep_flg) {
01260   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01261   
01262   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01263   if (this->_hasOp) {
01264     // Update MXj.
01265     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01266   }
01267   
01268   // Enter value on diagonal of B.
01269   (*B)(j,j) = diag;
01270       }
01271       else {
01272   // Create a random std::vector and orthogonalize it against all previous columns of Q.
01273   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01274   Teuchos::RCP<MV> tempMXj;
01275   MVT::MvRandom( *tempXj );
01276   if (this->_hasOp) {
01277     tempMXj = MVT::Clone( X, 1 );
01278     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01279   } 
01280   else {
01281     tempMXj = tempXj;
01282   }
01283   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01284   //
01285   for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
01286     
01287     for (int i=0; i<Q.length(); i++) {
01288       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01289 
01290       // Perform MGS one vector at a time
01291       for (int ii=0; ii<qcs[i]; ii++) {
01292         
01293         indX[0] = ii;
01294         tempQ = MVT::CloneView( *Q[i], indX );
01295         Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
01296       
01297         // Apply another step of modified Gram-Schmidt
01298         innerProd( *tempQ, *tempXj, tempMXj, tempC );
01299         MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
01300         
01301       }
01302 
01303       // Update MXj, with the least number of applications of Op as possible
01304       if (this->_hasOp) {
01305         if (MQ[i].get()) {
01306     // MQ was allocated and computed above; use it
01307     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01308         }
01309         else if (xc <= qcs[i]) {
01310     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01311     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01312         }
01313       }
01314     } // for (int i=0; i<nq; i++)   
01315   } // for (int num_orth=0; num_orth<max_orth_steps_; num_orth++)
01316   
01317   // Compute the Op-norms after the correction step.
01318   MVT::MvDot( *tempXj, *tempMXj, newDot );
01319   
01320   // Copy std::vector into current column of Xj
01321   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01322     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01323     
01324     // Enter value on diagonal of B.
01325     (*B)(j,j) = ZERO;
01326     
01327     // Copy std::vector into current column of _basisvecs
01328     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01329     if (this->_hasOp) {
01330       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01331     }
01332   }
01333   else {
01334     return j;
01335   } 
01336       } // if (!dep_flg)
01337       
01338       // Remove the vectors from array
01339       if (j > 0) {
01340   Q.resize( nq );
01341   C.resize( nq );
01342   qcs.resize( nq );
01343       }
01344       
01345     } // for (int j=0; j<xc; j++)
01346     
01347     return xc;
01348   }
01349   
01350 } // namespace Belos
01351 
01352 #endif // BELOS_IMGS_ORTHOMANAGER_HPP
01353 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines