Belos Version of the Day
BelosDGKSOrthoManager.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00047 #ifndef BELOS_DGKS_ORTHOMANAGER_HPP
00048 #define BELOS_DGKS_ORTHOMANAGER_HPP
00049 
00057 // #define ORTHO_DEBUG
00058 
00059 #include "BelosConfigDefs.hpp"
00060 #include "BelosMultiVecTraits.hpp"
00061 #include "BelosOperatorTraits.hpp"
00062 #include "BelosMatOrthoManager.hpp"
00063 
00064 namespace Belos {
00065 
00066   template<class ScalarType, class MV, class OP>
00067   class DGKSOrthoManager : 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     DGKSOrthoManager( const std::string& label = "Belos",
00082                       Teuchos::RCP<const OP> Op = Teuchos::null,
00083           const int max_blk_ortho = 2,
00084           const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00085           const MagnitudeType dep_tol = MGT::one()/MGT::squareroot( 2*MGT::one() ),
00086           const MagnitudeType sing_tol = 10*MGT::eps() )
00087       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00088   max_blk_ortho_( max_blk_ortho ),
00089   blk_tol_( blk_tol ),
00090   dep_tol_( dep_tol ),
00091   sing_tol_( sing_tol ),
00092   label_( label )
00093     {
00094         std::string orthoLabel = label_ + ": Orthogonalization";
00095         timerOrtho_ = Teuchos::TimeMonitor::getNewTimer( orthoLabel );
00096     }    
00097 
00099     ~DGKSOrthoManager() {}
00101 
00102 
00104 
00105 
00107     void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
00108 
00110     void setDepTol( const MagnitudeType dep_tol ) { dep_tol_ = dep_tol; }
00111 
00113     void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
00114 
00116     MagnitudeType getBlkTol() const { return blk_tol_; } 
00117 
00119     MagnitudeType getDepTol() const { return dep_tol_; } 
00120 
00122     MagnitudeType getSingTol() const { return sing_tol_; } 
00123 
00125 
00126 
00128 
00129 
00157     void project ( MV &X, Teuchos::RCP<MV> MX, 
00158                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00159                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00160 
00161 
00164     void project ( MV &X, 
00165                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00166                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00167       project(X,Teuchos::null,C,Q);
00168     }
00169 
00170 
00171  
00196     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00197                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00198 
00199 
00202     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00203       return normalize(X,Teuchos::null,B);
00204     }
00205 
00206 
00239     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00240                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00241                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00242                               Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00243 
00246     int projectAndNormalize ( MV &X, 
00247                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00248                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00249                               Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00250       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00251     }
00252 
00254 
00256 
00257 
00261     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00262     orthonormError(const MV &X) const {
00263       return orthonormError(X,Teuchos::null);
00264     }
00265 
00270     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00271     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00272 
00276     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00277     orthogError(const MV &X1, const MV &X2) const {
00278       return orthogError(X1,Teuchos::null,X2);
00279     }
00280 
00285     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00286     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00287 
00289 
00291 
00292 
00295     void setLabel(const std::string& label);
00296 
00299     const std::string& getLabel() const { return label_; }
00300 
00302 
00303   private:
00304     
00306     int max_blk_ortho_;
00307     MagnitudeType blk_tol_;
00308     MagnitudeType dep_tol_;
00309     MagnitudeType sing_tol_;
00310 
00312     std::string label_;
00313     Teuchos::RCP<Teuchos::Time> timerOrtho_;
00314 
00316     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00317       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00318       bool completeBasis, int howMany = -1 ) const;
00319     
00321     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00322         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00323         Teuchos::Array<Teuchos::RCP<const MV> > Q) const;
00324 
00325     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00326            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00327            Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00328            Teuchos::Array<Teuchos::RCP<const MV> > Q) const;    
00329   };
00330 
00332   // Set the label for this orthogonalization manager and create new timers if it's changed
00333   template<class ScalarType, class MV, class OP>
00334   void DGKSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00335   {
00336     if (label != label_) {
00337       label_ = label;
00338       std::string orthoLabel = label_ + ": Orthogonalization";
00339       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00340     }
00341   }
00342 
00344   // Compute the distance from orthonormality
00345   template<class ScalarType, class MV, class OP>
00346   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00347   DGKSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00348     const ScalarType ONE = SCT::one();
00349     int rank = MVT::GetNumberVecs(X);
00350     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00351     innerProd(X,X,MX,xTx);
00352     for (int i=0; i<rank; i++) {
00353       xTx(i,i) -= ONE;
00354     }
00355     return xTx.normFrobenius();
00356   }
00357 
00359   // Compute the distance from orthogonality
00360   template<class ScalarType, class MV, class OP>
00361   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00362   DGKSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00363     int r1 = MVT::GetNumberVecs(X1);
00364     int r2  = MVT::GetNumberVecs(X2);
00365     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00366     innerProd(X2,X1,MX1,xTx);
00367     return xTx.normFrobenius();
00368   }
00369 
00371   // Find an Op-orthonormal basis for span(X) - span(W)
00372   template<class ScalarType, class MV, class OP>
00373   int DGKSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00374                                     MV &X, Teuchos::RCP<MV> MX, 
00375                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00376                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00377                                     Teuchos::Array<Teuchos::RCP<const MV> > Q ) const {
00378     
00379     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00380 
00381     ScalarType    ONE  = SCT::one();
00382     ScalarType    ZERO  = SCT::zero();
00383 
00384     int nq = Q.length();
00385     int xc = MVT::GetNumberVecs( X );
00386     int xr = MVT::GetVecLength( X );
00387     int rank = xc;
00388 
00389     /* if the user doesn't want to store the coefficienets, 
00390      * allocate some local memory for them 
00391      */
00392     if ( B == Teuchos::null ) {
00393       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00394     }
00395 
00396     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00397     if (this->_hasOp) {
00398       if (MX == Teuchos::null) {
00399         // we need to allocate space for MX
00400         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00401         OPT::Apply(*(this->_Op),X,*MX);
00402       }
00403     }
00404     else {
00405       // Op == I  -->  MX = X (ignore it if the user passed it in)
00406       MX = Teuchos::rcp( &X, false );
00407     }
00408 
00409     int mxc = MVT::GetNumberVecs( *MX );
00410     int mxr = MVT::GetVecLength( *MX );
00411 
00412     // short-circuit
00413     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
00414 
00415     int numbas = 0;
00416     for (int i=0; i<nq; i++) {
00417       numbas += MVT::GetNumberVecs( *Q[i] );
00418     }
00419 
00420     // check size of B
00421     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00422                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00423     // check size of X and MX
00424     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00425                         "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00426     // check size of X w.r.t. MX 
00427     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00428                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00429     // check feasibility
00430     //TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00431     //                    "Belos::DGKSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00432 
00433     // Some flags for checking dependency returns from the internal orthogonalization methods
00434     bool dep_flg = false;
00435 
00436     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00437     Teuchos::RCP<MV> tmpX, tmpMX;
00438     tmpX = MVT::CloneCopy(X);
00439     if (this->_hasOp) {
00440       tmpMX = MVT::CloneCopy(*MX);
00441     }
00442 
00443     // Use the cheaper block orthogonalization.
00444     dep_flg = blkOrtho( X, MX, C, Q );
00445 
00446     // If a dependency has been detected in this block, then perform
00447     // the more expensive single-std::vector orthogonalization.
00448     if (dep_flg) {
00449       rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00450 
00451       // Copy tmpX back into X.
00452       MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00453       if (this->_hasOp) {
00454   MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00455       }
00456     } 
00457     else {
00458       // There is no dependency, so orthonormalize new block X
00459       rank = findBasis( X, MX, B, false );
00460       if (rank < xc) {
00461   // A dependency was found during orthonormalization of X,
00462   // rerun orthogonalization using more expensive single-std::vector orthogonalization.
00463   rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00464 
00465   // Copy tmpX back into X.
00466   MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00467   if (this->_hasOp) {
00468     MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00469   }
00470       }    
00471     }
00472 
00473     // this should not raise an std::exception; but our post-conditions oblige us to check
00474     TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00475                         "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00476 
00477     // Return the rank of X.
00478     return rank;
00479   }
00480 
00481 
00482 
00484   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00485   template<class ScalarType, class MV, class OP>
00486   int DGKSOrthoManager<ScalarType, MV, OP>::normalize(
00487                                 MV &X, Teuchos::RCP<MV> MX, 
00488                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00489 
00490     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00491 
00492     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00493     return findBasis(X, MX, B, true);
00494   }
00495 
00496 
00497 
00499   template<class ScalarType, class MV, class OP>
00500   void DGKSOrthoManager<ScalarType, MV, OP>::project(
00501                           MV &X, Teuchos::RCP<MV> MX, 
00502                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00503                           Teuchos::Array<Teuchos::RCP<const MV> > Q) const {
00504     // For the inner product defined by the operator Op or the identity (Op == 0)
00505     //   -> Orthogonalize X against each Q[i]
00506     // Modify MX accordingly
00507     //
00508     // Note that when Op is 0, MX is not referenced
00509     //
00510     // Parameter variables
00511     //
00512     // X  : Vectors to be transformed
00513     //
00514     // MX : Image of the block std::vector X by the mass matrix
00515     //
00516     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00517     //
00518 
00519     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00520 
00521     int xc = MVT::GetNumberVecs( X );
00522     int xr = MVT::GetVecLength( X );
00523     int nq = Q.length();
00524     std::vector<int> qcs(nq);
00525     // short-circuit
00526     if (nq == 0 || xc == 0 || xr == 0) {
00527       return;
00528     }
00529     int qr = MVT::GetVecLength ( *Q[0] );
00530     // if we don't have enough C, expand it with null references
00531     // if we have too many, resize to throw away the latter ones
00532     // if we have exactly as many as we have Q, this call has no effect
00533     C.resize(nq);
00534 
00535 
00536     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00537     if (this->_hasOp) {
00538       if (MX == Teuchos::null) {
00539         // we need to allocate space for MX
00540         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00541         OPT::Apply(*(this->_Op),X,*MX);
00542       }
00543     }
00544     else {
00545       // Op == I  -->  MX = X (ignore it if the user passed it in)
00546       MX = Teuchos::rcp( &X, false );
00547     }
00548     int mxc = MVT::GetNumberVecs( *MX );
00549     int mxr = MVT::GetVecLength( *MX );
00550 
00551     // check size of X and Q w.r.t. common sense
00552     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00553                         "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00554     // check size of X w.r.t. MX and Q
00555     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00556                         "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
00557 
00558     // tally up size of all Q and check/allocate C
00559     int baslen = 0;
00560     for (int i=0; i<nq; i++) {
00561       TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00562                           "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
00563       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00564       TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00565                           "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
00566       baslen += qcs[i];
00567 
00568       // check size of C[i]
00569       if ( C[i] == Teuchos::null ) {
00570         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00571       }
00572       else {
00573         TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00574                            "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
00575       }
00576     }
00577 
00578     // Use the cheaper block orthogonalization, don't check for rank deficiency.
00579     blkOrtho( X, MX, C, Q );
00580 
00581   }  
00582 
00584   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00585   // the rank is numvectors(X)
00586   template<class ScalarType, class MV, class OP>
00587   int DGKSOrthoManager<ScalarType, MV, OP>::findBasis(
00588                   MV &X, Teuchos::RCP<MV> MX, 
00589                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00590                   bool completeBasis, int howMany ) const {
00591     // For the inner product defined by the operator Op or the identity (Op == 0)
00592     //   -> Orthonormalize X 
00593     // Modify MX accordingly
00594     //
00595     // Note that when Op is 0, MX is not referenced
00596     //
00597     // Parameter variables
00598     //
00599     // X  : Vectors to be orthonormalized
00600     //
00601     // MX : Image of the multivector X under the operator Op
00602     //
00603     // Op  : Pointer to the operator for the inner product
00604     //
00605     //
00606 
00607     const ScalarType ONE  = SCT::one();
00608     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00609 
00610     int xc = MVT::GetNumberVecs( X );
00611     int xr = MVT::GetVecLength( X );
00612 
00613     if (howMany == -1) {
00614       howMany = xc;
00615     }
00616 
00617     /*******************************************************
00618      *  If _hasOp == false, we will not reference MX below *
00619      *******************************************************/
00620 
00621     // if Op==null, MX == X (via pointer)
00622     // Otherwise, either the user passed in MX or we will allocated and compute it
00623     if (this->_hasOp) {
00624       if (MX == Teuchos::null) {
00625         // we need to allocate space for MX
00626         MX = MVT::Clone(X,xc);
00627         OPT::Apply(*(this->_Op),X,*MX);
00628       }
00629     }
00630 
00631     /* if the user doesn't want to store the coefficienets, 
00632      * allocate some local memory for them 
00633      */
00634     if ( B == Teuchos::null ) {
00635       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00636     }
00637 
00638     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00639     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00640 
00641     // check size of C, B
00642     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00643                         "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
00644     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00645                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
00646     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00647                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
00648     TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00649                         "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
00650     TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
00651                         "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
00652 
00653     /* xstart is which column we are starting the process with, based on howMany
00654      * columns before xstart are assumed to be Op-orthonormal already
00655      */
00656     int xstart = xc - howMany;
00657 
00658     for (int j = xstart; j < xc; j++) {
00659 
00660       // numX is 
00661       // * number of currently orthonormal columns of X
00662       // * the index of the current column of X
00663       int numX = j;
00664       bool addVec = false;
00665 
00666       // Get a view of the std::vector currently being worked on.
00667       std::vector<int> index(1);
00668       index[0] = numX;
00669       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
00670       Teuchos::RCP<MV> MXj;
00671       if ((this->_hasOp)) {
00672         // MXj is a view of the current std::vector in MX
00673         MXj = MVT::CloneViewNonConst( *MX, index );
00674       }
00675       else {
00676         // MXj is a pointer to Xj, and MUST NOT be modified
00677         MXj = Xj;
00678       }
00679 
00680       // Get a view of the previous vectors.
00681       std::vector<int> prev_idx( numX );
00682       Teuchos::RCP<const MV> prevX, prevMX;
00683 
00684       if (numX > 0) {
00685         for (int i=0; i<numX; i++) {
00686           prev_idx[i] = i;
00687         }
00688         prevX = MVT::CloneView( X, prev_idx );
00689         if (this->_hasOp) {
00690           prevMX = MVT::CloneView( *MX, prev_idx );
00691         }
00692       } 
00693 
00694       // Make storage for these Gram-Schmidt iterations.
00695       Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
00696       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00697       //
00698       // Save old MXj std::vector and compute Op-norm
00699       //
00700       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00701       MVT::MvDot( *Xj, *MXj, oldDot );
00702       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
00703       TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
00704         "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
00705 
00706       if (numX > 0) {
00707   // Apply the first step of Gram-Schmidt
00708   
00709   // product <- prevX^T MXj
00710   innerProd(*prevX,*Xj,MXj,product);
00711   
00712   // Xj <- Xj - prevX prevX^T MXj   
00713   //     = Xj - prevX product
00714   MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
00715   
00716   // Update MXj
00717   if (this->_hasOp) {
00718     // MXj <- Op*Xj_new
00719     //      = Op*(Xj_old - prevX prevX^T MXj)
00720     //      = MXj - prevMX product
00721     MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
00722   }
00723   
00724   // Compute new Op-norm
00725   MVT::MvDot( *Xj, *MXj, newDot );
00726   
00727   // Check if a correction is needed.
00728   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(dep_tol_*oldDot[0]) ) {
00729     // Apply the second step of Gram-Schmidt
00730     // This is the same as above
00731     Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
00732     
00733     innerProd(*prevX,*Xj,MXj,P2);
00734     product += P2;
00735     MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00736     if ((this->_hasOp)) {
00737       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00738     }
00739   } // if (newDot[0] < dep_tol_*oldDot[0])
00740   
00741       } // if (numX > 0)
00742 
00743         // Compute Op-norm with old MXj
00744       MVT::MvDot( *Xj, *oldMXj, newDot );
00745 
00746       // Check to see if the new std::vector is dependent.
00747       if (completeBasis) {
00748   //
00749   // We need a complete basis, so add random vectors if necessary
00750   //
00751   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
00752     
00753     // Add a random std::vector and orthogonalize it against previous vectors in block.
00754     addVec = true;
00755 #ifdef ORTHO_DEBUG
00756     std::cout << "Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
00757 #endif
00758     //
00759     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
00760     Teuchos::RCP<MV> tempMXj;
00761     MVT::MvRandom( *tempXj );
00762     if (this->_hasOp) {
00763       tempMXj = MVT::Clone( X, 1 );
00764       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
00765     } 
00766     else {
00767       tempMXj = tempXj;
00768     }
00769     MVT::MvDot( *tempXj, *tempMXj, oldDot );
00770     //
00771     for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
00772       innerProd(*prevX,*tempXj,tempMXj,product);
00773       MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
00774       if (this->_hasOp) {
00775         MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
00776       }
00777     }
00778     // Compute new Op-norm
00779     MVT::MvDot( *tempXj, *tempMXj, newDot );
00780     //
00781     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
00782       // Copy std::vector into current column of _basisvecs
00783       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
00784       if (this->_hasOp) {
00785         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
00786       }
00787     }
00788     else {
00789       return numX;
00790     } 
00791   }
00792       }
00793       else {
00794   //
00795   // We only need to detect dependencies.
00796   //
00797   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
00798     return numX;
00799   }
00800       }
00801       
00802       // If we haven't left this method yet, then we can normalize the new std::vector Xj.
00803       // Normalize Xj.
00804       // Xj <- Xj / std::sqrt(newDot)
00805       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
00806 
00807       if (SCT::magnitude(diag) > ZERO) {      
00808         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
00809         if (this->_hasOp) {
00810     // Update MXj.
00811     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
00812         }
00813       }
00814 
00815       // If we've added a random std::vector, enter a zero in the j'th diagonal element.
00816       if (addVec) {
00817   (*B)(j,j) = ZERO;
00818       }
00819       else {
00820   (*B)(j,j) = diag;
00821       }
00822 
00823       // Save the coefficients, if we are working on the original std::vector and not a randomly generated one
00824       if (!addVec) {
00825   for (int i=0; i<numX; i++) {
00826     (*B)(i,j) = product(i,0);
00827   }
00828       }
00829 
00830     } // for (j = 0; j < xc; ++j)
00831 
00832     return xc;
00833   }
00834 
00836   // Routine to compute the block orthogonalization
00837   template<class ScalarType, class MV, class OP>
00838   bool 
00839   DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00840                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00841                Teuchos::Array<Teuchos::RCP<const MV> > Q) const
00842   {
00843     int nq = Q.length();
00844     int xc = MVT::GetNumberVecs( X );
00845     bool dep_flg = false;
00846     const ScalarType ONE  = SCT::one();
00847 
00848     std::vector<int> qcs( nq );
00849     for (int i=0; i<nq; i++) {
00850       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00851     }
00852 
00853     // Perform the Gram-Schmidt transformation for a block of vectors
00854     
00855     // Compute the initial Op-norms
00856     std::vector<ScalarType> oldDot( xc );
00857     MVT::MvDot( X, *MX, oldDot );
00858 
00859     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
00860     // Define the product Q^T * (Op*X)
00861     for (int i=0; i<nq; i++) {
00862       // Multiply Q' with MX
00863       innerProd(*Q[i],X,MX,*C[i]);
00864       // Multiply by Q and subtract the result in X
00865       MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
00866 
00867       // Update MX, with the least number of applications of Op as possible
00868       if (this->_hasOp) {
00869         if (xc <= qcs[i]) {
00870           OPT::Apply( *(this->_Op), X, *MX);
00871         }
00872         else {
00873           // this will possibly be used again below; don't delete it
00874           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
00875           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
00876           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
00877         }
00878       }
00879     }
00880 
00881     // Do as many steps of classical Gram-Schmidt as required by max_blk_ortho_
00882     for (int j = 1; j < max_blk_ortho_; ++j) {
00883       
00884       for (int i=0; i<nq; i++) {
00885   Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
00886         
00887   // Apply another step of classical Gram-Schmidt
00888   innerProd(*Q[i],X,MX,C2);
00889   *C[i] += C2;
00890   MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
00891         
00892   // Update MX, with the least number of applications of Op as possible
00893   if (this->_hasOp) {
00894     if (MQ[i].get()) {
00895       // MQ was allocated and computed above; use it
00896       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
00897     }
00898     else if (xc <= qcs[i]) {
00899       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
00900       OPT::Apply( *(this->_Op), X, *MX);
00901     }
00902   }
00903       } // for (int i=0; i<nq; i++)
00904     } // for (int j = 0; j < max_blk_ortho; ++j)
00905   
00906     // Compute new Op-norms
00907     std::vector<ScalarType> newDot(xc);
00908     MVT::MvDot( X, *MX, newDot );
00909  
00910     // Check to make sure the new block of vectors are not dependent on previous vectors
00911     for (int i=0; i<xc; i++){
00912       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
00913   dep_flg = true;
00914   break;
00915       }
00916     } // end for (i=0;...)
00917 
00918     return dep_flg;
00919   }
00920   
00922   // Routine to compute the block orthogonalization using single-std::vector orthogonalization
00923   template<class ScalarType, class MV, class OP>
00924   int
00925   DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00926                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00927                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00928                    Teuchos::Array<Teuchos::RCP<const MV> > Q) const
00929   {
00930     const ScalarType ONE  = SCT::one();
00931     const ScalarType ZERO  = SCT::zero();
00932     
00933     int nq = Q.length();
00934     int xc = MVT::GetNumberVecs( X );
00935     std::vector<int> indX( 1 );
00936     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00937 
00938     std::vector<int> qcs( nq );
00939     for (int i=0; i<nq; i++) {
00940       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00941     }
00942 
00943     // Create pointers for the previous vectors of X that have already been orthonormalized.
00944     Teuchos::RCP<const MV> lastQ;
00945     Teuchos::RCP<MV> Xj, MXj;
00946     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
00947 
00948     // Perform the Gram-Schmidt transformation for each std::vector in the block of vectors.
00949     for (int j=0; j<xc; j++) {
00950       
00951       bool dep_flg = false;
00952       
00953       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
00954       if (j > 0) {
00955   std::vector<int> index( j );
00956   for (int ind=0; ind<j; ind++) {
00957     index[ind] = ind;
00958   }
00959   lastQ = MVT::CloneView( X, index );
00960 
00961   // Add these views to the Q and C arrays.
00962   Q.push_back( lastQ );
00963   C.push_back( B );
00964   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
00965       }
00966       
00967       // Get a view of the current std::vector in X to orthogonalize.
00968       indX[0] = j;
00969       Xj = MVT::CloneViewNonConst( X, indX );
00970       if (this->_hasOp) {
00971   MXj = MVT::CloneViewNonConst( *MX, indX );
00972       }
00973       else {
00974   MXj = Xj;
00975       }
00976       
00977       // Compute the initial Op-norms
00978       MVT::MvDot( *Xj, *MXj, oldDot );
00979       
00980       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.length());
00981       // Define the product Q^T * (Op*X)
00982       for (int i=0; i<Q.length(); i++) {
00983 
00984   // Get a view of the current serial dense matrix
00985   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
00986 
00987   // Multiply Q' with MX
00988   innerProd(*Q[i],*Xj,MXj,tempC);
00989   // Multiply by Q and subtract the result in Xj
00990   MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
00991   
00992   // Update MXj, with the least number of applications of Op as possible
00993   if (this->_hasOp) {
00994     if (xc <= qcs[i]) {
00995       OPT::Apply( *(this->_Op), *Xj, *MXj);
00996     }
00997     else {
00998       // this will possibly be used again below; don't delete it
00999       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01000       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01001       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01002     }
01003   }
01004       }
01005       
01006       // Compute the Op-norms
01007       MVT::MvDot( *Xj, *MXj, newDot );
01008       
01009       // Do one step of classical Gram-Schmidt orthogonalization 
01010       // with a second correction step if needed.
01011       
01012       if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
01013   
01014   for (int i=0; i<Q.length(); i++) {
01015     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01016     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01017     
01018     // Apply another step of classical Gram-Schmidt
01019     innerProd(*Q[i],*Xj,MXj,C2);
01020     tempC += C2;
01021     MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
01022     
01023     // Update MXj, with the least number of applications of Op as possible
01024     if (this->_hasOp) {
01025       if (MQ[i].get()) {
01026         // MQ was allocated and computed above; use it
01027         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01028       }
01029       else if (xc <= qcs[i]) {
01030         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01031         OPT::Apply( *(this->_Op), *Xj, *MXj);
01032       }
01033     }
01034   } // for (int i=0; i<Q.length(); i++)
01035   
01036   // Compute the Op-norms after the correction step.
01037   MVT::MvDot( *Xj, *MXj, newDot );
01038   
01039       } // if ()
01040       
01041       // Check for linear dependence.
01042       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01043   dep_flg = true;
01044       }
01045       
01046       // Normalize the new std::vector if it's not dependent
01047       if (!dep_flg) {
01048   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01049   
01050   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01051   if (this->_hasOp) {
01052     // Update MXj.
01053     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01054   }
01055   
01056   // Enter value on diagonal of B.
01057   (*B)(j,j) = diag;
01058       }
01059       else {
01060   // Create a random std::vector and orthogonalize it against all previous columns of Q.
01061   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01062   Teuchos::RCP<MV> tempMXj;
01063   MVT::MvRandom( *tempXj );
01064   if (this->_hasOp) {
01065     tempMXj = MVT::Clone( X, 1 );
01066     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01067   } 
01068   else {
01069     tempMXj = tempXj;
01070   }
01071   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01072   //
01073   for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
01074     
01075     for (int i=0; i<Q.length(); i++) {
01076       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01077       
01078       // Apply another step of classical Gram-Schmidt
01079       innerProd(*Q[i],*tempXj,tempMXj,product);
01080       MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
01081       
01082       // Update MXj, with the least number of applications of Op as possible
01083       if (this->_hasOp) {
01084         if (MQ[i].get()) {
01085     // MQ was allocated and computed above; use it
01086     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01087         }
01088         else if (xc <= qcs[i]) {
01089     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01090     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01091         }
01092       }
01093     } // for (int i=0; i<nq; i++)
01094     
01095   }
01096   
01097   // Compute the Op-norms after the correction step.
01098   MVT::MvDot( *tempXj, *tempMXj, newDot );
01099   
01100   // Copy std::vector into current column of Xj
01101   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01102     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01103     
01104     // Enter value on diagonal of B.
01105     (*B)(j,j) = ZERO;
01106 
01107     // Copy std::vector into current column of _basisvecs
01108     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01109     if (this->_hasOp) {
01110       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01111     }
01112   }
01113   else {
01114     return j;
01115   } 
01116       } // if (!dep_flg)
01117 
01118       // Remove the vectors from array
01119       if (j > 0) {
01120   Q.resize( nq );
01121   C.resize( nq );
01122   qcs.resize( nq );
01123       }
01124 
01125     } // for (int j=0; j<xc; j++)
01126 
01127     return xc;
01128   }
01129 
01130 } // namespace Belos
01131 
01132 #endif // BELOS_DGKS_ORTHOMANAGER_HPP
01133 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines