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 terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 
00034 #ifndef BELOS_IMGS_ORTHOMANAGER_HPP
00035 #define BELOS_IMGS_ORTHOMANAGER_HPP
00036 
00044 // #define ORTHO_DEBUG
00045 
00046 #include "BelosConfigDefs.hpp"
00047 #include "BelosMultiVecTraits.hpp"
00048 #include "BelosOperatorTraits.hpp"
00049 #include "BelosMatOrthoManager.hpp"
00050 
00051 namespace Belos {
00052 
00053   template<class ScalarType, class MV, class OP>
00054   class IMGSOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00055 
00056   private:
00057     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00058     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00059     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00060     typedef MultiVecTraits<ScalarType,MV>      MVT;
00061     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00062 
00063   public:
00064     
00066 
00067 
00068     IMGSOrthoManager( const std::string& label = "Belos",
00069                       Teuchos::RCP<const OP> Op = Teuchos::null,
00070           const int max_ortho_steps = 1,
00071           const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00072           const MagnitudeType sing_tol = 10*MGT::eps() )
00073       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00074   max_ortho_steps_( max_ortho_steps ),
00075   blk_tol_( blk_tol ),
00076   sing_tol_( sing_tol ),
00077         label_( label )
00078     {
00079         std::string orthoLabel = label_ + ": Orthogonalization";
00080         timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00081 
00082         std::string updateLabel = label_ + ": Ortho (Update)";
00083         timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00084 
00085         std::string normLabel = label_ + ": Ortho (Norm)";
00086         timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00087 
00088         std::string ipLabel = label_ + ": Ortho (Inner Product)";
00089         timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel); 
00090     };
00091 
00093     ~IMGSOrthoManager() {}
00095 
00096 
00098 
00099 
00101     void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
00102 
00104     void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
00105 
00107     MagnitudeType getBlkTol() const { return blk_tol_; } 
00108 
00110     MagnitudeType getSingTol() const { return sing_tol_; } 
00111 
00113 
00114 
00116 
00117 
00145     void project ( MV &X, Teuchos::RCP<MV> MX, 
00146                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00147                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00148 
00149 
00152     void project ( MV &X, 
00153                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00154                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00155       project(X,Teuchos::null,C,Q);
00156     }
00157 
00158 
00159  
00184     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00185                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00186 
00187 
00190     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00191       return normalize(X,Teuchos::null,B);
00192     }
00193 
00194 
00227     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00228                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00229                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00230                               Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00231 
00234     int projectAndNormalize ( MV &X, 
00235                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00236                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00237                               Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00238       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00239     }
00240 
00242 
00244 
00245 
00249     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00250     orthonormError(const MV &X) const {
00251       return orthonormError(X,Teuchos::null);
00252     }
00253 
00258     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00259     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00260 
00264     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00265     orthogError(const MV &X1, const MV &X2) const {
00266       return orthogError(X1,Teuchos::null,X2);
00267     }
00268 
00273     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00274     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00275 
00277 
00279 
00280 
00283     void setLabel(const std::string& label);
00284 
00287     const std::string& getLabel() const { return label_; }
00288 
00290 
00291   private:
00292     
00294     int max_ortho_steps_;
00295     MagnitudeType blk_tol_;
00296     MagnitudeType sing_tol_;
00297 
00299     std::string label_;
00300     Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, 
00301                                         timerNorm_, timerScale_, timerInnerProd_;
00302   
00304     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00305       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00306       bool completeBasis, int howMany = -1 ) const;
00307     
00309     bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
00310          Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00311          Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00312 
00314     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00315         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00316         Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00317 
00318     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00319            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00320            Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00321            Teuchos::Array<Teuchos::RCP<const MV> > Q) const;    
00322   };
00323 
00325   // Set the label for this orthogonalization manager and create new timers if it's changed
00326   template<class ScalarType, class MV, class OP>
00327   void IMGSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00328   {   
00329     if (label != label_) {
00330       label_ = label;
00331       std::string orthoLabel = label_ + ": Orthogonalization";
00332       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00333 
00334       std::string updateLabel = label_ + ": Ortho (Update)";
00335       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00336 
00337       std::string normLabel = label_ + ": Ortho (Norm)";
00338       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00339 
00340       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00341       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel);
00342     }
00343   } 
00344 
00346   // Compute the distance from orthonormality
00347   template<class ScalarType, class MV, class OP>
00348   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00349   IMGSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00350     const ScalarType ONE = SCT::one();
00351     int rank = MVT::GetNumberVecs(X);
00352     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00353     innerProd(X,X,MX,xTx);
00354     for (int i=0; i<rank; i++) {
00355       xTx(i,i) -= ONE;
00356     }
00357     return xTx.normFrobenius();
00358   }
00359 
00361   // Compute the distance from orthogonality
00362   template<class ScalarType, class MV, class OP>
00363   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00364   IMGSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00365     int r1 = MVT::GetNumberVecs(X1);
00366     int r2  = MVT::GetNumberVecs(X2);
00367     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00368     innerProd(X2,X1,MX1,xTx);
00369     return xTx.normFrobenius();
00370   }
00371 
00373   // Find an Op-orthonormal basis for span(X) - span(W)
00374   template<class ScalarType, class MV, class OP>
00375   int IMGSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00376                                     MV &X, Teuchos::RCP<MV> MX, 
00377                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00378                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00379                                     Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00380 
00381     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00382 
00383     ScalarType    ONE  = SCT::one();
00384     ScalarType    ZERO  = SCT::zero();
00385 
00386     int nq = Q.length();
00387     int xc = MVT::GetNumberVecs( X );
00388     int xr = MVT::GetVecLength( X );
00389     int rank = xc;
00390 
00391     /* if the user doesn't want to store the coefficienets, 
00392      * allocate some local memory for them 
00393      */
00394     if ( B == Teuchos::null ) {
00395       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00396     }
00397 
00398     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00399     if (this->_hasOp) {
00400       if (MX == Teuchos::null) {
00401         // we need to allocate space for MX
00402         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00403         OPT::Apply(*(this->_Op),X,*MX);
00404       }
00405     }
00406     else {
00407       // Op == I  -->  MX = X (ignore it if the user passed it in)
00408       MX = Teuchos::rcp( &X, false );
00409     }
00410 
00411     int mxc = MVT::GetNumberVecs( *MX );
00412     int mxr = MVT::GetVecLength( *MX );
00413 
00414     // short-circuit
00415     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
00416 
00417     int numbas = 0;
00418     for (int i=0; i<nq; i++) {
00419       numbas += MVT::GetNumberVecs( *Q[i] );
00420     }
00421 
00422     // check size of B
00423     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00424                         "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00425     // check size of X and MX
00426     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00427                         "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00428     // check size of X w.r.t. MX 
00429     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00430                         "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00431     // check feasibility
00432     //TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00433     //                    "Belos::IMGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00434 
00435     // Some flags for checking dependency returns from the internal orthogonalization methods
00436     bool dep_flg = false;
00437 
00438     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00439     Teuchos::RCP<MV> tmpX, tmpMX;
00440     tmpX = MVT::CloneCopy(X);
00441     if (this->_hasOp) {
00442       tmpMX = MVT::CloneCopy(*MX);
00443     }
00444 
00445     if (xc == 1) {
00446 
00447       // Use the cheaper block orthogonalization.
00448       // NOTE: Don't check for dependencies because the update has one std::vector.
00449       dep_flg = blkOrtho1( X, MX, C, Q );
00450 
00451       // Normalize the new block X
00452       {
00453         Teuchos::TimeMonitor normTimer( *timerNorm_ );
00454         if ( B == Teuchos::null ) {
00455           B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00456         }
00457         std::vector<ScalarType> diag(xc);
00458         MVT::MvDot( X, *MX, diag );
00459         (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
00460         rank = 1; 
00461         MVT::MvAddMv( ONE/(*B)(0,0), X, ZERO, X, X );
00462         if (this->_hasOp) {
00463           // Update MXj.
00464     MVT::MvAddMv( ONE/(*B)(0,0), *MX, ZERO, *MX, *MX );
00465         }
00466       }
00467 
00468     } 
00469     else {
00470 
00471       // Use the cheaper block orthogonalization.
00472       dep_flg = blkOrtho( X, MX, C, Q );
00473 
00474       // If a dependency has been detected in this block, then perform
00475       // the more expensive single-std::vector orthogonalization.
00476       if (dep_flg) {
00477         rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00478 
00479         // Copy tmpX back into X.
00480         MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00481         if (this->_hasOp) {
00482     MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00483         }
00484       } 
00485       else {
00486         // There is no dependency, so orthonormalize new block X
00487         rank = findBasis( X, MX, B, false );
00488         if (rank < xc) {
00489     // A dependency was found during orthonormalization of X,
00490     // rerun orthogonalization using more expensive single-std::vector orthogonalization.
00491     rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00492 
00493     // Copy tmpX back into X.
00494     MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00495     if (this->_hasOp) {
00496       MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00497     }
00498         }    
00499       }
00500     } // if (xc == 1) {
00501 
00502     // this should not raise an std::exception; but our post-conditions oblige us to check
00503     TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00504                         "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00505 
00506     // Return the rank of X.
00507     return rank;
00508   }
00509 
00510 
00511 
00513   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00514   template<class ScalarType, class MV, class OP>
00515   int IMGSOrthoManager<ScalarType, MV, OP>::normalize(
00516                                 MV &X, Teuchos::RCP<MV> MX, 
00517                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00518 
00519     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00520 
00521     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00522     return findBasis(X, MX, B, true);
00523   }
00524 
00525 
00526 
00528   template<class ScalarType, class MV, class OP>
00529   void IMGSOrthoManager<ScalarType, MV, OP>::project(
00530                           MV &X, Teuchos::RCP<MV> MX, 
00531                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00532                           Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00533     // For the inner product defined by the operator Op or the identity (Op == 0)
00534     //   -> Orthogonalize X against each Q[i]
00535     // Modify MX accordingly
00536     //
00537     // Note that when Op is 0, MX is not referenced
00538     //
00539     // Parameter variables
00540     //
00541     // X  : Vectors to be transformed
00542     //
00543     // MX : Image of the block std::vector X by the mass matrix
00544     //
00545     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00546     //
00547     
00548     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00549     
00550     int xc = MVT::GetNumberVecs( X );
00551     int xr = MVT::GetVecLength( X );
00552     int nq = Q.length();
00553     std::vector<int> qcs(nq);
00554     // short-circuit
00555     if (nq == 0 || xc == 0 || xr == 0) {
00556       return;
00557     }
00558     int qr = MVT::GetVecLength ( *Q[0] );
00559     // if we don't have enough C, expand it with null references
00560     // if we have too many, resize to throw away the latter ones
00561     // if we have exactly as many as we have Q, this call has no effect
00562     C.resize(nq);
00563 
00564 
00565     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00566     if (this->_hasOp) {
00567       if (MX == Teuchos::null) {
00568         // we need to allocate space for MX
00569         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00570         OPT::Apply(*(this->_Op),X,*MX);
00571       }
00572     }
00573     else {
00574       // Op == I  -->  MX = X (ignore it if the user passed it in)
00575       MX = Teuchos::rcp( &X, false );
00576     }
00577     int mxc = MVT::GetNumberVecs( *MX );
00578     int mxr = MVT::GetVecLength( *MX );
00579 
00580     // check size of X and Q w.r.t. common sense
00581     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00582                         "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00583     // check size of X w.r.t. MX and Q
00584     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00585                         "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
00586 
00587     // tally up size of all Q and check/allocate C
00588     int baslen = 0;
00589     for (int i=0; i<nq; i++) {
00590       TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00591                           "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
00592       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00593       TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00594                           "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
00595       baslen += qcs[i];
00596 
00597       // check size of C[i]
00598       if ( C[i] == Teuchos::null ) {
00599         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00600       }
00601       else {
00602         TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00603                            "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
00604       }
00605     }
00606 
00607     // Use the cheaper block orthogonalization, don't check for rank deficiency.
00608     blkOrtho( X, MX, C, Q );
00609 
00610   }  
00611 
00613   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00614   // the rank is numvectors(X)
00615   template<class ScalarType, class MV, class OP>
00616   int IMGSOrthoManager<ScalarType, MV, OP>::findBasis(
00617                   MV &X, Teuchos::RCP<MV> MX, 
00618                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00619                   bool completeBasis, int howMany ) const {
00620     // For the inner product defined by the operator Op or the identity (Op == 0)
00621     //   -> Orthonormalize X 
00622     // Modify MX accordingly
00623     //
00624     // Note that when Op is 0, MX is not referenced
00625     //
00626     // Parameter variables
00627     //
00628     // X  : Vectors to be orthonormalized
00629     //
00630     // MX : Image of the multivector X under the operator Op
00631     //
00632     // Op  : Pointer to the operator for the inner product
00633     //
00634     //
00635 
00636     Teuchos::TimeMonitor normTimer( *timerNorm_ );
00637 
00638     const ScalarType ONE  = SCT::one();
00639     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00640 
00641     int xc = MVT::GetNumberVecs( X );
00642     int xr = MVT::GetVecLength( X );
00643 
00644     if (howMany == -1) {
00645       howMany = xc;
00646     }
00647 
00648     /*******************************************************
00649      *  If _hasOp == false, we will not reference MX below *
00650      *******************************************************/
00651 
00652     // if Op==null, MX == X (via pointer)
00653     // Otherwise, either the user passed in MX or we will allocated and compute it
00654     if (this->_hasOp) {
00655       if (MX == Teuchos::null) {
00656         // we need to allocate space for MX
00657         MX = MVT::Clone(X,xc);
00658         OPT::Apply(*(this->_Op),X,*MX);
00659       }
00660     }
00661 
00662     /* if the user doesn't want to store the coefficienets, 
00663      * allocate some local memory for them 
00664      */
00665     if ( B == Teuchos::null ) {
00666       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00667     }
00668 
00669     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00670     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00671 
00672     // check size of C, B
00673     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00674                         "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
00675     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00676                         "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
00677     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00678                         "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
00679     TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00680                         "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
00681     TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
00682                         "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
00683 
00684     /* xstart is which column we are starting the process with, based on howMany
00685      * columns before xstart are assumed to be Op-orthonormal already
00686      */
00687     int xstart = xc - howMany;
00688 
00689     for (int j = xstart; j < xc; j++) {
00690 
00691       // numX is 
00692       // * number of currently orthonormal columns of X
00693       // * the index of the current column of X
00694       int numX = j;
00695       bool addVec = false;
00696 
00697       // Get a view of the std::vector currently being worked on.
00698       std::vector<int> index(1);
00699       index[0] = numX;
00700       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
00701       Teuchos::RCP<MV> MXj;
00702       if ((this->_hasOp)) {
00703         // MXj is a view of the current std::vector in MX
00704         MXj = MVT::CloneViewNonConst( *MX, index );
00705       }
00706       else {
00707         // MXj is a pointer to Xj, and MUST NOT be modified
00708         MXj = Xj;
00709       }
00710 
00711       // Make storage for these Gram-Schmidt iterations.
00712       Teuchos::SerialDenseVector<int,ScalarType> product(numX);
00713       Teuchos::SerialDenseVector<int,ScalarType> P2(1);
00714       Teuchos::RCP<const MV> prevX, prevMX;
00715   
00716       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00717       //
00718       // Save old MXj std::vector and compute Op-norm
00719       //
00720       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00721       MVT::MvDot( *Xj, *MXj, oldDot );
00722       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
00723       TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
00724         "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
00725 
00726       // Perform MGS one vector at a time
00727       for (int ii=0; ii<numX; ii++) {
00728   
00729   index[0] = ii;
00730   prevX = MVT::CloneView( X, index );
00731   if (this->_hasOp) {
00732     prevMX = MVT::CloneView( *MX, index );
00733   }
00734   
00735   for (int i=0; i<max_ortho_steps_; ++i) {
00736     
00737     // product <- prevX^T MXj
00738     {
00739       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00740       innerProd(*prevX,*Xj,MXj,P2);
00741     }
00742     
00743     // Xj <- Xj - prevX prevX^T MXj   
00744     //     = Xj - prevX product
00745     {
00746       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00747       MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00748     }
00749     
00750     // Update MXj
00751     if (this->_hasOp) {
00752       // MXj <- Op*Xj_new
00753       //      = Op*(Xj_old - prevX prevX^T MXj)
00754       //      = MXj - prevMX product
00755       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00756       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00757     }
00758     
00759     // Set coefficients
00760     if ( i==0 )
00761       product[ii] = P2[0];
00762     else
00763       product[ii] += P2[0];
00764     
00765   } // for (int i=0; i<max_ortho_steps_; ++i)
00766 
00767       } // for (int ii=0; ii<numX; ++ii)  
00768   
00769   // Compute Op-norm with old MXj
00770       MVT::MvDot( *Xj, *oldMXj, newDot );
00771       
00772       // Check to see if the new std::vector is dependent.
00773       if (completeBasis) {
00774   //
00775   // We need a complete basis, so add random vectors if necessary
00776   //
00777   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
00778     
00779     // Add a random std::vector and orthogonalize it against previous vectors in block.
00780     addVec = true;
00781 #ifdef ORTHO_DEBUG
00782     std::cout << "Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
00783 #endif
00784     //
00785     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
00786     Teuchos::RCP<MV> tempMXj;
00787     MVT::MvRandom( *tempXj );
00788     if (this->_hasOp) {
00789       tempMXj = MVT::Clone( X, 1 );
00790       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
00791     } 
00792     else {
00793       tempMXj = tempXj;
00794     }
00795     MVT::MvDot( *tempXj, *tempMXj, oldDot );
00796     //
00797     // Perform MGS one vector at a time
00798     for (int ii=0; ii<numX; ii++) {
00799       
00800       index[0] = ii;
00801       prevX = MVT::CloneView( X, index );
00802       if (this->_hasOp) {
00803         prevMX = MVT::CloneView( *MX, index );
00804       }
00805       
00806       for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
00807         {
00808     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00809     innerProd(*prevX,*tempXj,tempMXj,P2);
00810         }
00811         {
00812     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00813     MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
00814         }
00815         if (this->_hasOp) {
00816     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00817     MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
00818         }
00819 
00820         // Set coefficients
00821         if ( num_orth==0 )
00822     product[ii] = P2[0];
00823         else
00824     product[ii] += P2[0];       
00825       }
00826     }
00827       
00828     // Compute new Op-norm
00829     MVT::MvDot( *tempXj, *tempMXj, newDot );
00830     //
00831     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
00832       // Copy std::vector into current column of _basisvecs
00833       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
00834       if (this->_hasOp) {
00835         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
00836       }
00837     }
00838     else {
00839       return numX;
00840     } 
00841   }
00842       }
00843       else {
00844   //
00845   // We only need to detect dependencies.
00846   //
00847   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
00848     return numX;
00849   }
00850       }
00851       
00852       
00853       // If we haven't left this method yet, then we can normalize the new 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);
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   IMGSOrthoManager<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     std::vector<int> index(1);
00904     Teuchos::RCP<const MV> tempQ;
00905 
00906     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
00907     // Define the product Q^T * (Op*X)
00908     for (int i=0; i<nq; i++) {
00909 
00910       // Perform MGS one vector at a time
00911       for (int ii=0; ii<qcs[i]; ii++) {
00912   
00913   index[0] = ii;
00914   tempQ = MVT::CloneView( *Q[i], index );
00915   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
00916 
00917   // Multiply Q' with MX
00918   {
00919     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00920     innerProd(*tempQ,X,MX,tempC);
00921   }
00922   // Multiply by Q and subtract the result in X
00923   {
00924     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00925     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
00926   }
00927       }
00928       
00929       // Update MX, with the least number of applications of Op as possible
00930       if (this->_hasOp) {
00931   if (xc <= qcs[i]) {
00932     OPT::Apply( *(this->_Op), X, *MX);
00933   }
00934         else {
00935           // this will possibly be used again below; don't delete it
00936           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
00937           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
00938           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
00939         }
00940       }
00941     }
00942     
00943     // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
00944     for (int j = 1; j < max_ortho_steps_; ++j) {
00945       
00946       for (int i=0; i<nq; i++) {
00947 
00948   Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
00949         
00950   // Perform MGS one vector at a time
00951   for (int ii=0; ii<qcs[i]; ii++) {
00952     
00953     index[0] = ii;
00954     tempQ = MVT::CloneView( *Q[i], index );
00955     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
00956     Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
00957     
00958     // Apply another step of modified Gram-Schmidt
00959     {
00960       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00961       innerProd( *tempQ, X, MX, tempC2 );
00962     }
00963     tempC += tempC2;
00964     {
00965       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00966       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
00967     }
00968     
00969   }
00970 
00971   // Update MX, with the least number of applications of Op as possible
00972   if (this->_hasOp) {
00973     if (MQ[i].get()) {
00974       // MQ was allocated and computed above; use it
00975       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
00976     }
00977     else if (xc <= qcs[i]) {
00978       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
00979       OPT::Apply( *(this->_Op), X, *MX);
00980     }
00981   }
00982       } // for (int i=0; i<nq; i++)
00983     } // for (int j = 0; j < max_ortho_steps; ++j)
00984   
00985     return false;
00986   }
00987 
00989   // Routine to compute the block orthogonalization
00990   template<class ScalarType, class MV, class OP>
00991   bool 
00992   IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00993                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00994                Teuchos::Array<Teuchos::RCP<const MV> > Q) const
00995   {
00996     int nq = Q.length();
00997     int xc = MVT::GetNumberVecs( X );
00998     bool dep_flg = false;
00999     const ScalarType ONE  = SCT::one();
01000 
01001     std::vector<int> qcs( nq );
01002     for (int i=0; i<nq; i++) {
01003       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01004     }
01005 
01006     // Perform the Gram-Schmidt transformation for a block of vectors
01007     
01008     // Compute the initial Op-norms
01009     std::vector<ScalarType> oldDot( xc );
01010     MVT::MvDot( X, *MX, oldDot );
01011 
01012     std::vector<int> index(1);
01013     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01014     Teuchos::RCP<const MV> tempQ;
01015 
01016     // Define the product Q^T * (Op*X)
01017     for (int i=0; i<nq; i++) {
01018 
01019       // Perform MGS one vector at a time
01020       for (int ii=0; ii<qcs[i]; ii++) {
01021   
01022   index[0] = ii;
01023   tempQ = MVT::CloneView( *Q[i], index );
01024   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
01025 
01026   // Multiply Q' with MX
01027   {
01028     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01029     innerProd( *tempQ, X, MX, tempC);
01030   }
01031   // Multiply by Q and subtract the result in X
01032   {
01033     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01034     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
01035   }
01036       }
01037 
01038       // Update MX, with the least number of applications of Op as possible
01039       if (this->_hasOp) {
01040         if (xc <= qcs[i]) {
01041           OPT::Apply( *(this->_Op), X, *MX);
01042         }
01043         else {
01044           // this will possibly be used again below; don't delete it
01045           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01046           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01047           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01048         }
01049       }
01050     }
01051 
01052     // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
01053     for (int j = 1; j < max_ortho_steps_; ++j) {
01054       
01055       for (int i=0; i<nq; i++) {
01056   Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
01057         
01058   // Perform MGS one vector at a time
01059   for (int ii=0; ii<qcs[i]; ii++) {
01060     
01061     index[0] = ii;
01062     tempQ = MVT::CloneView( *Q[i], index );
01063     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
01064     Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
01065     
01066     // Apply another step of modified Gram-Schmidt
01067     {
01068       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01069       innerProd( *tempQ, X, MX, tempC2 );
01070     }
01071     tempC += tempC2;
01072     {
01073       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01074       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
01075     }
01076   }
01077 
01078   // Update MX, with the least number of applications of Op as possible
01079   if (this->_hasOp) {
01080     if (MQ[i].get()) {
01081       // MQ was allocated and computed above; use it
01082       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01083     }
01084     else if (xc <= qcs[i]) {
01085       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01086       OPT::Apply( *(this->_Op), X, *MX);
01087     }
01088   }
01089       } // for (int i=0; i<nq; i++)
01090     } // for (int j = 0; j < max_ortho_steps; ++j)
01091   
01092     // Compute new Op-norms
01093     std::vector<ScalarType> newDot(xc);
01094     MVT::MvDot( X, *MX, newDot );
01095     
01096     // Check to make sure the new block of vectors are not dependent on previous vectors
01097     for (int i=0; i<xc; i++){
01098       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01099   dep_flg = true;
01100   break;
01101       }
01102     } // end for (i=0;...)
01103     
01104     return dep_flg;
01105   }
01106   
01108   // Routine to compute the block orthogonalization using single-std::vector orthogonalization
01109   template<class ScalarType, class MV, class OP>
01110   int
01111   IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01112                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01113                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01114                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const
01115   {
01116     const ScalarType ONE  = SCT::one();
01117     const ScalarType ZERO  = SCT::zero();
01118     
01119     int nq = Q.length();
01120     int xc = MVT::GetNumberVecs( X );
01121     std::vector<int> indX( 1 );
01122     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01123 
01124     std::vector<int> qcs( nq );
01125     for (int i=0; i<nq; i++) {
01126       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01127     }
01128 
01129     // Create pointers for the previous vectors of X that have already been orthonormalized.
01130     Teuchos::RCP<const MV> lastQ;
01131     Teuchos::RCP<MV> Xj, MXj;
01132     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01133 
01134     // Perform the Gram-Schmidt transformation for each std::vector in the block of vectors.
01135     for (int j=0; j<xc; j++) {
01136       
01137       bool dep_flg = false;
01138       
01139       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01140       if (j > 0) {
01141   std::vector<int> index( j );
01142   for (int ind=0; ind<j; ind++) {
01143     index[ind] = ind;
01144   }
01145   lastQ = MVT::CloneView( X, index );
01146 
01147   // Add these views to the Q and C arrays.
01148   Q.push_back( lastQ );
01149   C.push_back( B );
01150   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01151       }
01152       
01153       // Get a view of the current std::vector in X to orthogonalize.
01154       indX[0] = j;
01155       Xj = MVT::CloneViewNonConst( X, indX );
01156       if (this->_hasOp) {
01157   MXj = MVT::CloneViewNonConst( *MX, indX );
01158       }
01159       else {
01160   MXj = Xj;
01161       }
01162       
01163       // Compute the initial Op-norms
01164       MVT::MvDot( *Xj, *MXj, oldDot );
01165       
01166       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.length());
01167       Teuchos::RCP<const MV> tempQ;
01168 
01169       // Define the product Q^T * (Op*X)
01170       for (int i=0; i<Q.length(); i++) {
01171 
01172   // Perform MGS one vector at a time
01173   for (int ii=0; ii<qcs[i]; ii++) {
01174     
01175     indX[0] = ii;
01176     tempQ = MVT::CloneView( *Q[i], indX );
01177     // Get a view of the current serial dense matrix
01178     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
01179     
01180     // Multiply Q' with MX
01181     innerProd(*tempQ,*Xj,MXj,tempC);
01182 
01183     // Multiply by Q and subtract the result in Xj
01184     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
01185   }
01186  
01187   // Update MXj, with the least number of applications of Op as possible
01188   if (this->_hasOp) {
01189     if (xc <= qcs[i]) {
01190       OPT::Apply( *(this->_Op), *Xj, *MXj);
01191     }
01192     else {
01193       // this will possibly be used again below; don't delete it
01194       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01195       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01196       Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01197       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01198     }
01199   }
01200       }
01201       
01202       // Do any additional steps of modified Gram-Schmidt orthogonalization 
01203       for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
01204   
01205   for (int i=0; i<Q.length(); i++) {
01206     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01207     
01208     // Perform MGS one vector at a time
01209     for (int ii=0; ii<qcs[i]; ii++) {
01210       
01211       indX[0] = ii;
01212       tempQ = MVT::CloneView( *Q[i], indX );
01213       // Get a view of the current serial dense matrix
01214       Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
01215       
01216       // Apply another step of modified Gram-Schmidt
01217       innerProd( *tempQ, *Xj, MXj, tempC2);
01218       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj ); 
01219     }
01220     
01221     // Add the coefficients into C[i]
01222     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01223     tempC += C2;
01224     
01225     // Update MXj, with the least number of applications of Op as possible
01226     if (this->_hasOp) {
01227       if (MQ[i].get()) {
01228         // MQ was allocated and computed above; use it
01229         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01230       }
01231       else if (xc <= qcs[i]) {
01232         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01233         OPT::Apply( *(this->_Op), *Xj, *MXj);
01234       }
01235     }
01236   } // for (int i=0; i<Q.length(); i++)
01237   
01238       } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
01239       
01240       // Check for linear dependence.
01241       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01242   dep_flg = true;
01243       }
01244       
01245       // Normalize the new std::vector if it's not dependent
01246       if (!dep_flg) {
01247   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01248   
01249   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01250   if (this->_hasOp) {
01251     // Update MXj.
01252     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01253   }
01254   
01255   // Enter value on diagonal of B.
01256   (*B)(j,j) = diag;
01257       }
01258       else {
01259   // Create a random std::vector and orthogonalize it against all previous columns of Q.
01260   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01261   Teuchos::RCP<MV> tempMXj;
01262   MVT::MvRandom( *tempXj );
01263   if (this->_hasOp) {
01264     tempMXj = MVT::Clone( X, 1 );
01265     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01266   } 
01267   else {
01268     tempMXj = tempXj;
01269   }
01270   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01271   //
01272   for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
01273     
01274     for (int i=0; i<Q.length(); i++) {
01275       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01276 
01277       // Perform MGS one vector at a time
01278       for (int ii=0; ii<qcs[i]; ii++) {
01279         
01280         indX[0] = ii;
01281         tempQ = MVT::CloneView( *Q[i], indX );
01282         Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
01283       
01284         // Apply another step of modified Gram-Schmidt
01285         innerProd( *tempQ, *tempXj, tempMXj, tempC );
01286         MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
01287         
01288       }
01289 
01290       // Update MXj, with the least number of applications of Op as possible
01291       if (this->_hasOp) {
01292         if (MQ[i].get()) {
01293     // MQ was allocated and computed above; use it
01294     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01295         }
01296         else if (xc <= qcs[i]) {
01297     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01298     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01299         }
01300       }
01301     } // for (int i=0; i<nq; i++)   
01302   } // for (int num_orth=0; num_orth<max_orth_steps_; num_orth++)
01303   
01304   // Compute the Op-norms after the correction step.
01305   MVT::MvDot( *tempXj, *tempMXj, newDot );
01306   
01307   // Copy std::vector into current column of Xj
01308   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01309     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01310     
01311     // Enter value on diagonal of B.
01312     (*B)(j,j) = ZERO;
01313     
01314     // Copy std::vector into current column of _basisvecs
01315     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01316     if (this->_hasOp) {
01317       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01318     }
01319   }
01320   else {
01321     return j;
01322   } 
01323       } // if (!dep_flg)
01324       
01325       // Remove the vectors from array
01326       if (j > 0) {
01327   Q.resize( nq );
01328   C.resize( nq );
01329   qcs.resize( nq );
01330       }
01331       
01332     } // for (int j=0; j<xc; j++)
01333     
01334     return xc;
01335   }
01336   
01337 } // namespace Belos
01338 
01339 #endif // BELOS_IMGS_ORTHOMANAGER_HPP
01340 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:06 2011 for Belos by  doxygen 1.6.3