BelosDGKSOrthoManager.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under 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_DGKS_ORTHOMANAGER_HPP
00035 #define BELOS_DGKS_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 DGKSOrthoManager : 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     DGKSOrthoManager( const std::string& label = "Belos",
00069                       Teuchos::RCP<const OP> Op = Teuchos::null,
00070           const int max_blk_ortho = 2,
00071           const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00072           const MagnitudeType dep_tol = MGT::one()/MGT::squareroot( 2*MGT::one() ),
00073           const MagnitudeType sing_tol = 10*MGT::eps() )
00074       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00075   max_blk_ortho_( max_blk_ortho ),
00076   blk_tol_( blk_tol ),
00077   dep_tol_( dep_tol ),
00078   sing_tol_( sing_tol ),
00079   label_( label )
00080     {
00081         std::string orthoLabel = label_ + ": Orthogonalization";
00082         timerOrtho_ = Teuchos::TimeMonitor::getNewTimer( orthoLabel );
00083     }    
00084 
00086     ~DGKSOrthoManager() {}
00088 
00089 
00091 
00092 
00094     void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
00095 
00097     void setDepTol( const MagnitudeType dep_tol ) { dep_tol_ = dep_tol; }
00098 
00100     void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
00101 
00103     MagnitudeType getBlkTol() const { return blk_tol_; } 
00104 
00106     MagnitudeType getDepTol() const { return dep_tol_; } 
00107 
00109     MagnitudeType getSingTol() const { return sing_tol_; } 
00110 
00112 
00113 
00115 
00116 
00144     void project ( MV &X, Teuchos::RCP<MV> MX, 
00145                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00146                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00147 
00148 
00151     void project ( MV &X, 
00152                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00153                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00154       project(X,Teuchos::null,C,Q);
00155     }
00156 
00157 
00158  
00183     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00184                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00185 
00186 
00189     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00190       return normalize(X,Teuchos::null,B);
00191     }
00192 
00193 
00226     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00227                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00228                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00229                               Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00230 
00233     int projectAndNormalize ( MV &X, 
00234                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00235                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00236                               Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00237       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00238     }
00239 
00241 
00243 
00244 
00248     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00249     orthonormError(const MV &X) const {
00250       return orthonormError(X,Teuchos::null);
00251     }
00252 
00257     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00258     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00259 
00263     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00264     orthogError(const MV &X1, const MV &X2) const {
00265       return orthogError(X1,Teuchos::null,X2);
00266     }
00267 
00272     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00273     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00274 
00276 
00278 
00279 
00282     void setLabel(const std::string& label);
00283 
00286     const std::string& getLabel() const { return label_; }
00287 
00289 
00290   private:
00291     
00293     int max_blk_ortho_;
00294     MagnitudeType blk_tol_;
00295     MagnitudeType dep_tol_;
00296     MagnitudeType sing_tol_;
00297 
00299     std::string label_;
00300     Teuchos::RCP<Teuchos::Time> timerOrtho_;
00301 
00303     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00304       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00305       bool completeBasis, int howMany = -1 ) const;
00306     
00308     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00309         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00310         Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00311 
00312     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00313            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00314            Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00315            Teuchos::Array<Teuchos::RCP<const MV> > Q) const;    
00316   };
00317 
00319   // Set the label for this orthogonalization manager and create new timers if it's changed
00320   template<class ScalarType, class MV, class OP>
00321   void DGKSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00322   {
00323     if (label != label_) {
00324       label_ = label;
00325       std::string orthoLabel = label_ + ": Orthogonalization";
00326       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00327     }
00328   }
00329 
00331   // Compute the distance from orthonormality
00332   template<class ScalarType, class MV, class OP>
00333   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00334   DGKSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00335     const ScalarType ONE = SCT::one();
00336     int rank = MVT::GetNumberVecs(X);
00337     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00338     innerProd(X,X,MX,xTx);
00339     for (int i=0; i<rank; i++) {
00340       xTx(i,i) -= ONE;
00341     }
00342     return xTx.normFrobenius();
00343   }
00344 
00346   // Compute the distance from orthogonality
00347   template<class ScalarType, class MV, class OP>
00348   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00349   DGKSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00350     int r1 = MVT::GetNumberVecs(X1);
00351     int r2  = MVT::GetNumberVecs(X2);
00352     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00353     innerProd(X2,X1,MX1,xTx);
00354     return xTx.normFrobenius();
00355   }
00356 
00358   // Find an Op-orthonormal basis for span(X) - span(W)
00359   template<class ScalarType, class MV, class OP>
00360   int DGKSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00361                                     MV &X, Teuchos::RCP<MV> MX, 
00362                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00363                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00364                                     Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00365     
00366     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00367 
00368     ScalarType    ONE  = SCT::one();
00369     ScalarType    ZERO  = SCT::zero();
00370 
00371     int nq = Q.length();
00372     int xc = MVT::GetNumberVecs( X );
00373     int xr = MVT::GetVecLength( X );
00374     int rank = xc;
00375 
00376     /* if the user doesn't want to store the coefficienets, 
00377      * allocate some local memory for them 
00378      */
00379     if ( B == Teuchos::null ) {
00380       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00381     }
00382 
00383     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00384     if (this->_hasOp) {
00385       if (MX == Teuchos::null) {
00386         // we need to allocate space for MX
00387         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00388         OPT::Apply(*(this->_Op),X,*MX);
00389       }
00390     }
00391     else {
00392       // Op == I  -->  MX = X (ignore it if the user passed it in)
00393       MX = Teuchos::rcp( &X, false );
00394     }
00395 
00396     int mxc = MVT::GetNumberVecs( *MX );
00397     int mxr = MVT::GetVecLength( *MX );
00398 
00399     // short-circuit
00400     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
00401 
00402     int numbas = 0;
00403     for (int i=0; i<nq; i++) {
00404       numbas += MVT::GetNumberVecs( *Q[i] );
00405     }
00406 
00407     // check size of B
00408     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00409                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00410     // check size of X and MX
00411     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00412                         "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00413     // check size of X w.r.t. MX 
00414     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00415                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00416     // check feasibility
00417     //TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00418     //                    "Belos::DGKSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00419 
00420     // Some flags for checking dependency returns from the internal orthogonalization methods
00421     bool dep_flg = false;
00422 
00423     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00424     Teuchos::RCP<MV> tmpX, tmpMX;
00425     tmpX = MVT::CloneCopy(X);
00426     if (this->_hasOp) {
00427       tmpMX = MVT::CloneCopy(*MX);
00428     }
00429 
00430     // Use the cheaper block orthogonalization.
00431     dep_flg = blkOrtho( X, MX, C, Q );
00432 
00433     // If a dependency has been detected in this block, then perform
00434     // the more expensive single-std::vector orthogonalization.
00435     if (dep_flg) {
00436       rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00437 
00438       // Copy tmpX back into X.
00439       MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00440       if (this->_hasOp) {
00441   MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00442       }
00443     } 
00444     else {
00445       // There is no dependency, so orthonormalize new block X
00446       rank = findBasis( X, MX, B, false );
00447       if (rank < xc) {
00448   // A dependency was found during orthonormalization of X,
00449   // rerun orthogonalization using more expensive single-std::vector orthogonalization.
00450   rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00451 
00452   // Copy tmpX back into X.
00453   MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00454   if (this->_hasOp) {
00455     MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00456   }
00457       }    
00458     }
00459 
00460     // this should not raise an std::exception; but our post-conditions oblige us to check
00461     TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00462                         "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00463 
00464     // Return the rank of X.
00465     return rank;
00466   }
00467 
00468 
00469 
00471   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00472   template<class ScalarType, class MV, class OP>
00473   int DGKSOrthoManager<ScalarType, MV, OP>::normalize(
00474                                 MV &X, Teuchos::RCP<MV> MX, 
00475                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00476 
00477     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00478 
00479     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00480     return findBasis(X, MX, B, true);
00481   }
00482 
00483 
00484 
00486   template<class ScalarType, class MV, class OP>
00487   void DGKSOrthoManager<ScalarType, MV, OP>::project(
00488                           MV &X, Teuchos::RCP<MV> MX, 
00489                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00490                           Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00491     // For the inner product defined by the operator Op or the identity (Op == 0)
00492     //   -> Orthogonalize X against each Q[i]
00493     // Modify MX accordingly
00494     //
00495     // Note that when Op is 0, MX is not referenced
00496     //
00497     // Parameter variables
00498     //
00499     // X  : Vectors to be transformed
00500     //
00501     // MX : Image of the block std::vector X by the mass matrix
00502     //
00503     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00504     //
00505 
00506     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00507 
00508     int xc = MVT::GetNumberVecs( X );
00509     int xr = MVT::GetVecLength( X );
00510     int nq = Q.length();
00511     std::vector<int> qcs(nq);
00512     // short-circuit
00513     if (nq == 0 || xc == 0 || xr == 0) {
00514       return;
00515     }
00516     int qr = MVT::GetVecLength ( *Q[0] );
00517     // if we don't have enough C, expand it with null references
00518     // if we have too many, resize to throw away the latter ones
00519     // if we have exactly as many as we have Q, this call has no effect
00520     C.resize(nq);
00521 
00522 
00523     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00524     if (this->_hasOp) {
00525       if (MX == Teuchos::null) {
00526         // we need to allocate space for MX
00527         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00528         OPT::Apply(*(this->_Op),X,*MX);
00529       }
00530     }
00531     else {
00532       // Op == I  -->  MX = X (ignore it if the user passed it in)
00533       MX = Teuchos::rcp( &X, false );
00534     }
00535     int mxc = MVT::GetNumberVecs( *MX );
00536     int mxr = MVT::GetVecLength( *MX );
00537 
00538     // check size of X and Q w.r.t. common sense
00539     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00540                         "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00541     // check size of X w.r.t. MX and Q
00542     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00543                         "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
00544 
00545     // tally up size of all Q and check/allocate C
00546     int baslen = 0;
00547     for (int i=0; i<nq; i++) {
00548       TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00549                           "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
00550       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00551       TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00552                           "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
00553       baslen += qcs[i];
00554 
00555       // check size of C[i]
00556       if ( C[i] == Teuchos::null ) {
00557         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00558       }
00559       else {
00560         TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00561                            "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
00562       }
00563     }
00564 
00565     // Use the cheaper block orthogonalization, don't check for rank deficiency.
00566     blkOrtho( X, MX, C, Q );
00567 
00568   }  
00569 
00571   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00572   // the rank is numvectors(X)
00573   template<class ScalarType, class MV, class OP>
00574   int DGKSOrthoManager<ScalarType, MV, OP>::findBasis(
00575                   MV &X, Teuchos::RCP<MV> MX, 
00576                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00577                   bool completeBasis, int howMany ) const {
00578     // For the inner product defined by the operator Op or the identity (Op == 0)
00579     //   -> Orthonormalize X 
00580     // Modify MX accordingly
00581     //
00582     // Note that when Op is 0, MX is not referenced
00583     //
00584     // Parameter variables
00585     //
00586     // X  : Vectors to be orthonormalized
00587     //
00588     // MX : Image of the multivector X under the operator Op
00589     //
00590     // Op  : Pointer to the operator for the inner product
00591     //
00592     //
00593 
00594     const ScalarType ONE  = SCT::one();
00595     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00596 
00597     int xc = MVT::GetNumberVecs( X );
00598     int xr = MVT::GetVecLength( X );
00599 
00600     if (howMany == -1) {
00601       howMany = xc;
00602     }
00603 
00604     /*******************************************************
00605      *  If _hasOp == false, we will not reference MX below *
00606      *******************************************************/
00607 
00608     // if Op==null, MX == X (via pointer)
00609     // Otherwise, either the user passed in MX or we will allocated and compute it
00610     if (this->_hasOp) {
00611       if (MX == Teuchos::null) {
00612         // we need to allocate space for MX
00613         MX = MVT::Clone(X,xc);
00614         OPT::Apply(*(this->_Op),X,*MX);
00615       }
00616     }
00617 
00618     /* if the user doesn't want to store the coefficienets, 
00619      * allocate some local memory for them 
00620      */
00621     if ( B == Teuchos::null ) {
00622       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00623     }
00624 
00625     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00626     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00627 
00628     // check size of C, B
00629     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00630                         "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
00631     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00632                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
00633     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00634                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
00635     TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00636                         "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
00637     TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
00638                         "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
00639 
00640     /* xstart is which column we are starting the process with, based on howMany
00641      * columns before xstart are assumed to be Op-orthonormal already
00642      */
00643     int xstart = xc - howMany;
00644 
00645     for (int j = xstart; j < xc; j++) {
00646 
00647       // numX is 
00648       // * number of currently orthonormal columns of X
00649       // * the index of the current column of X
00650       int numX = j;
00651       bool addVec = false;
00652 
00653       // Get a view of the std::vector currently being worked on.
00654       std::vector<int> index(1);
00655       index[0] = numX;
00656       Teuchos::RCP<MV> Xj = MVT::CloneView( X, index );
00657       Teuchos::RCP<MV> MXj;
00658       if ((this->_hasOp)) {
00659         // MXj is a view of the current std::vector in MX
00660         MXj = MVT::CloneView( *MX, index );
00661       }
00662       else {
00663         // MXj is a pointer to Xj, and MUST NOT be modified
00664         MXj = Xj;
00665       }
00666 
00667       // Get a view of the previous vectors.
00668       std::vector<int> prev_idx( numX );
00669       Teuchos::RCP<const MV> prevX, prevMX;
00670 
00671       if (numX > 0) {
00672         for (int i=0; i<numX; i++) {
00673           prev_idx[i] = i;
00674         }
00675         prevX = MVT::CloneView( X, prev_idx );
00676         if (this->_hasOp) {
00677           prevMX = MVT::CloneView( *MX, prev_idx );
00678         }
00679       } 
00680 
00681       // Make storage for these Gram-Schmidt iterations.
00682       Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
00683       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00684       //
00685       // Save old MXj std::vector and compute Op-norm
00686       //
00687       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00688       MVT::MvDot( *Xj, *MXj, oldDot );
00689       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
00690       TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
00691         "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
00692 
00693       if (numX > 0) {
00694   // Apply the first step of Gram-Schmidt
00695   
00696   // product <- prevX^T MXj
00697   innerProd(*prevX,*Xj,MXj,product);
00698   
00699   // Xj <- Xj - prevX prevX^T MXj   
00700   //     = Xj - prevX product
00701   MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
00702   
00703   // Update MXj
00704   if (this->_hasOp) {
00705     // MXj <- Op*Xj_new
00706     //      = Op*(Xj_old - prevX prevX^T MXj)
00707     //      = MXj - prevMX product
00708     MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
00709   }
00710   
00711   // Compute new Op-norm
00712   MVT::MvDot( *Xj, *MXj, newDot );
00713   
00714   // Check if a correction is needed.
00715   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(dep_tol_*oldDot[0]) ) {
00716     // Apply the second step of Gram-Schmidt
00717     // This is the same as above
00718     Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
00719     
00720     innerProd(*prevX,*Xj,MXj,P2);
00721     product += P2;
00722     MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00723     if ((this->_hasOp)) {
00724       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00725     }
00726   } // if (newDot[0] < dep_tol_*oldDot[0])
00727   
00728       } // if (numX > 0)
00729 
00730         // Compute Op-norm with old MXj
00731       MVT::MvDot( *Xj, *oldMXj, newDot );
00732 
00733       // Check to see if the new std::vector is dependent.
00734       if (completeBasis) {
00735   //
00736   // We need a complete basis, so add random vectors if necessary
00737   //
00738   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
00739     
00740     // Add a random std::vector and orthogonalize it against previous vectors in block.
00741     addVec = true;
00742 #ifdef ORTHO_DEBUG
00743     std::cout << "Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
00744 #endif
00745     //
00746     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
00747     Teuchos::RCP<MV> tempMXj;
00748     MVT::MvRandom( *tempXj );
00749     if (this->_hasOp) {
00750       tempMXj = MVT::Clone( X, 1 );
00751       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
00752     } 
00753     else {
00754       tempMXj = tempXj;
00755     }
00756     MVT::MvDot( *tempXj, *tempMXj, oldDot );
00757     //
00758     for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
00759       innerProd(*prevX,*tempXj,tempMXj,product);
00760       MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
00761       if (this->_hasOp) {
00762         MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
00763       }
00764     }
00765     // Compute new Op-norm
00766     MVT::MvDot( *tempXj, *tempMXj, newDot );
00767     //
00768     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
00769       // Copy std::vector into current column of _basisvecs
00770       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
00771       if (this->_hasOp) {
00772         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
00773       }
00774     }
00775     else {
00776       return numX;
00777     } 
00778   }
00779       }
00780       else {
00781   //
00782   // We only need to detect dependencies.
00783   //
00784   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
00785     return numX;
00786   }
00787       }
00788       
00789       // If we haven't left this method yet, then we can normalize the new std::vector Xj.
00790       // Normalize Xj.
00791       // Xj <- Xj / std::sqrt(newDot)
00792       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
00793 
00794       if (SCT::magnitude(diag) > ZERO) {      
00795         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
00796         if (this->_hasOp) {
00797     // Update MXj.
00798     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
00799         }
00800       }
00801 
00802       // If we've added a random std::vector, enter a zero in the j'th diagonal element.
00803       if (addVec) {
00804   (*B)(j,j) = ZERO;
00805       }
00806       else {
00807   (*B)(j,j) = diag;
00808       }
00809 
00810       // Save the coefficients, if we are working on the original std::vector and not a randomly generated one
00811       if (!addVec) {
00812   for (int i=0; i<numX; i++) {
00813     (*B)(i,j) = product(i,0);
00814   }
00815       }
00816 
00817     } // for (j = 0; j < xc; ++j)
00818 
00819     return xc;
00820   }
00821 
00823   // Routine to compute the block orthogonalization
00824   template<class ScalarType, class MV, class OP>
00825   bool 
00826   DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00827                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00828                Teuchos::Array<Teuchos::RCP<const MV> > Q) const
00829   {
00830     int nq = Q.length();
00831     int xc = MVT::GetNumberVecs( X );
00832     bool dep_flg = false;
00833     const ScalarType ONE  = SCT::one();
00834 
00835     std::vector<int> qcs( nq );
00836     for (int i=0; i<nq; i++) {
00837       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00838     }
00839 
00840     // Perform the Gram-Schmidt transformation for a block of vectors
00841     
00842     // Compute the initial Op-norms
00843     std::vector<ScalarType> oldDot( xc );
00844     MVT::MvDot( X, *MX, oldDot );
00845 
00846     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
00847     // Define the product Q^T * (Op*X)
00848     for (int i=0; i<nq; i++) {
00849       // Multiply Q' with MX
00850       innerProd(*Q[i],X,MX,*C[i]);
00851       // Multiply by Q and subtract the result in X
00852       MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
00853 
00854       // Update MX, with the least number of applications of Op as possible
00855       if (this->_hasOp) {
00856         if (xc <= qcs[i]) {
00857           OPT::Apply( *(this->_Op), X, *MX);
00858         }
00859         else {
00860           // this will possibly be used again below; don't delete it
00861           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
00862           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
00863           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
00864         }
00865       }
00866     }
00867 
00868     // Do as many steps of classical Gram-Schmidt as required by max_blk_ortho_
00869     for (int j = 1; j < max_blk_ortho_; ++j) {
00870       
00871       for (int i=0; i<nq; i++) {
00872   Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
00873         
00874   // Apply another step of classical Gram-Schmidt
00875   innerProd(*Q[i],X,MX,C2);
00876   *C[i] += C2;
00877   MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
00878         
00879   // Update MX, with the least number of applications of Op as possible
00880   if (this->_hasOp) {
00881     if (MQ[i].get()) {
00882       // MQ was allocated and computed above; use it
00883       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
00884     }
00885     else if (xc <= qcs[i]) {
00886       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
00887       OPT::Apply( *(this->_Op), X, *MX);
00888     }
00889   }
00890       } // for (int i=0; i<nq; i++)
00891     } // for (int j = 0; j < max_blk_ortho; ++j)
00892   
00893     // Compute new Op-norms
00894     std::vector<ScalarType> newDot(xc);
00895     MVT::MvDot( X, *MX, newDot );
00896  
00897     // Check to make sure the new block of vectors are not dependent on previous vectors
00898     for (int i=0; i<xc; i++){
00899       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
00900   dep_flg = true;
00901   break;
00902       }
00903     } // end for (i=0;...)
00904 
00905     return dep_flg;
00906   }
00907   
00909   // Routine to compute the block orthogonalization using single-std::vector orthogonalization
00910   template<class ScalarType, class MV, class OP>
00911   int
00912   DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00913                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00914                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00915                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const
00916   {
00917     const ScalarType ONE  = SCT::one();
00918     const ScalarType ZERO  = SCT::zero();
00919     
00920     int nq = Q.length();
00921     int xc = MVT::GetNumberVecs( X );
00922     std::vector<int> indX( 1 );
00923     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00924 
00925     std::vector<int> qcs( nq );
00926     for (int i=0; i<nq; i++) {
00927       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00928     }
00929 
00930     // Create pointers for the previous vectors of X that have already been orthonormalized.
00931     Teuchos::RCP<const MV> lastQ;
00932     Teuchos::RCP<MV> Xj, MXj;
00933     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
00934 
00935     // Perform the Gram-Schmidt transformation for each std::vector in the block of vectors.
00936     for (int j=0; j<xc; j++) {
00937       
00938       bool dep_flg = false;
00939       
00940       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
00941       if (j > 0) {
00942   std::vector<int> index( j );
00943   for (int ind=0; ind<j; ind++) {
00944     index[ind] = ind;
00945   }
00946   lastQ = MVT::CloneView( X, index );
00947 
00948   // Add these views to the Q and C arrays.
00949   Q.push_back( lastQ );
00950   C.push_back( B );
00951   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
00952       }
00953       
00954       // Get a view of the current std::vector in X to orthogonalize.
00955       indX[0] = j;
00956       Xj = MVT::CloneView( X, indX );
00957       if (this->_hasOp) {
00958   MXj = MVT::CloneView( *MX, indX );
00959       }
00960       else {
00961   MXj = Xj;
00962       }
00963       
00964       // Compute the initial Op-norms
00965       MVT::MvDot( *Xj, *MXj, oldDot );
00966       
00967       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.length());
00968       // Define the product Q^T * (Op*X)
00969       for (int i=0; i<Q.length(); i++) {
00970 
00971   // Get a view of the current serial dense matrix
00972   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
00973 
00974   // Multiply Q' with MX
00975   innerProd(*Q[i],*Xj,MXj,tempC);
00976   // Multiply by Q and subtract the result in Xj
00977   MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
00978   
00979   // Update MXj, with the least number of applications of Op as possible
00980   if (this->_hasOp) {
00981     if (xc <= qcs[i]) {
00982       OPT::Apply( *(this->_Op), *Xj, *MXj);
00983     }
00984     else {
00985       // this will possibly be used again below; don't delete it
00986       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
00987       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
00988       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
00989     }
00990   }
00991       }
00992       
00993       // Compute the Op-norms
00994       MVT::MvDot( *Xj, *MXj, newDot );
00995       
00996       // Do one step of classical Gram-Schmidt orthogonalization 
00997       // with a second correction step if needed.
00998       
00999       if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
01000   
01001   for (int i=0; i<Q.length(); i++) {
01002     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01003     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01004     
01005     // Apply another step of classical Gram-Schmidt
01006     innerProd(*Q[i],*Xj,MXj,C2);
01007     tempC += C2;
01008     MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
01009     
01010     // Update MXj, with the least number of applications of Op as possible
01011     if (this->_hasOp) {
01012       if (MQ[i].get()) {
01013         // MQ was allocated and computed above; use it
01014         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01015       }
01016       else if (xc <= qcs[i]) {
01017         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01018         OPT::Apply( *(this->_Op), *Xj, *MXj);
01019       }
01020     }
01021   } // for (int i=0; i<Q.length(); i++)
01022   
01023   // Compute the Op-norms after the correction step.
01024   MVT::MvDot( *Xj, *MXj, newDot );
01025   
01026       } // if ()
01027       
01028       // Check for linear dependence.
01029       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01030   dep_flg = true;
01031       }
01032       
01033       // Normalize the new std::vector if it's not dependent
01034       if (!dep_flg) {
01035   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01036   
01037   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01038   if (this->_hasOp) {
01039     // Update MXj.
01040     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01041   }
01042   
01043   // Enter value on diagonal of B.
01044   (*B)(j,j) = diag;
01045       }
01046       else {
01047   // Create a random std::vector and orthogonalize it against all previous columns of Q.
01048   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01049   Teuchos::RCP<MV> tempMXj;
01050   MVT::MvRandom( *tempXj );
01051   if (this->_hasOp) {
01052     tempMXj = MVT::Clone( X, 1 );
01053     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01054   } 
01055   else {
01056     tempMXj = tempXj;
01057   }
01058   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01059   //
01060   for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
01061     
01062     for (int i=0; i<Q.length(); i++) {
01063       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01064       
01065       // Apply another step of classical Gram-Schmidt
01066       innerProd(*Q[i],*tempXj,tempMXj,product);
01067       MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
01068       
01069       // Update MXj, with the least number of applications of Op as possible
01070       if (this->_hasOp) {
01071         if (MQ[i].get()) {
01072     // MQ was allocated and computed above; use it
01073     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01074         }
01075         else if (xc <= qcs[i]) {
01076     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01077     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01078         }
01079       }
01080     } // for (int i=0; i<nq; i++)
01081     
01082   }
01083   
01084   // Compute the Op-norms after the correction step.
01085   MVT::MvDot( *tempXj, *tempMXj, newDot );
01086   
01087   // Copy std::vector into current column of Xj
01088   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01089     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01090     
01091     // Enter value on diagonal of B.
01092     (*B)(j,j) = ZERO;
01093 
01094     // Copy std::vector into current column of _basisvecs
01095     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01096     if (this->_hasOp) {
01097       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01098     }
01099   }
01100   else {
01101     return j;
01102   } 
01103       } // if (!dep_flg)
01104 
01105       // Remove the vectors from array
01106       if (j > 0) {
01107   Q.resize( nq );
01108   C.resize( nq );
01109   qcs.resize( nq );
01110       }
01111 
01112     } // for (int j=0; j<xc; j++)
01113 
01114     return xc;
01115   }
01116 
01117 } // namespace Belos
01118 
01119 #endif // BELOS_DGKS_ORTHOMANAGER_HPP
01120 

Generated on Tue Oct 20 12:48:34 2009 for Belos by doxygen 1.4.7