Belos Version of the Day
BelosIMGSOrthoManager.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00047 #ifndef BELOS_IMGS_ORTHOMANAGER_HPP
00048 #define BELOS_IMGS_ORTHOMANAGER_HPP
00049 
00057 // #define ORTHO_DEBUG
00058 
00059 #include "BelosConfigDefs.hpp"
00060 #include "BelosMultiVecTraits.hpp"
00061 #include "BelosOperatorTraits.hpp"
00062 #include "BelosMatOrthoManager.hpp"
00063 #include "Teuchos_SerialDenseMatrix.hpp"
00064 #include "Teuchos_SerialDenseVector.hpp"
00065 
00066 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00067 #include "Teuchos_TimeMonitor.hpp"
00068 #endif // BELOS_TEUCHOS_TIME_MONITOR
00069 
00070 namespace Belos {
00071 
00075   template<class ScalarType>
00076   Teuchos::RCP<const Teuchos::ParameterList> 
00077   getDefaultImgsParameters()
00078   {
00079     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00080     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00081 
00082     // This part makes this class method non-reentrant.
00083     static Teuchos::RCP<Teuchos::ParameterList> params;
00084     if (! params.is_null())
00085       return params;
00086 
00087     // Default parameter values for IMGS orthogonalization.
00088     // Documentation will be embedded in the parameter list.
00089     const int defaultMaxNumOrthogPasses = 2;
00090     const magnitude_type eps = STM::eps();
00091     const magnitude_type defaultBlkTol = magnitude_type(10) * STM::squareroot(eps);
00092     const magnitude_type defaultSingTol = magnitude_type(10) * eps;
00093 
00094     params = Teuchos::parameterList();
00095     params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00096      "Maximum number of orthogonalization passes "
00097      "(includes the first).  Default is 2, since "
00098      "\"twice is enough\" for Krylov methods.");
00099     params->set ("blkTol", defaultBlkTol, 
00100      "Block reorthogonalization threshhold.");
00101     params->set ("singTol", defaultSingTol, 
00102      "Singular block detection threshold.");
00103     return params;
00104   }
00105 
00109   template<class ScalarType>
00110   Teuchos::RCP<const Teuchos::ParameterList> 
00111   getFastImgsParameters()
00112   {
00113     using Teuchos::ParameterList;
00114     using Teuchos::RCP;
00115     using Teuchos::rcp;
00116     using Teuchos::ScalarTraits;
00117     typedef typename ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00118     typedef ScalarTraits<magnitude_type> STM;
00119 
00120     // This part makes this class method non-reentrant.
00121     static RCP<ParameterList> params;
00122     if (params.is_null())
00123       {
00124   RCP<const ParameterList> defaultParams = getDefaultImgsParameters<ScalarType>();
00125   // Start with a clone of the default parameters
00126   params = rcp (new ParameterList (*defaultParams));
00127 
00128   const int maxBlkOrtho = 1;
00129   params->set ("maxNumOrthogPasses", maxBlkOrtho);
00130 
00131   const magnitude_type blkTol = STM::zero();
00132   params->set ("blkTol", blkTol);
00133 
00134   const magnitude_type singTol = STM::zero();
00135   params->set ("singTol", singTol);
00136       }
00137     return params;
00138   }
00139 
00145   template<class ScalarType>
00146   void
00147   readImgsParameters (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00148           int& maxNumOrthogPasses,
00149           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& blkTol,
00150           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& singTol)
00151   {
00152     using Teuchos::ParameterList;
00153     using Teuchos::RCP;
00154     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00155     const magnitude_type zero = Teuchos::ScalarTraits<magnitude_type>::zero();
00156     RCP<const ParameterList> defaultParams = getDefaultImgsParameters<ScalarType>();
00157 
00158     // Using temporary variables and fetching all values before
00159     // setting the output arguments ensures the strong exception
00160     // guarantee for this function: if an exception is thrown, no
00161     // externally visible side effects (in this case, setting the
00162     // output arguments) have taken place.
00163     int _maxNumOrthogPasses;
00164     magnitude_type _blkTol, _singTol;
00165     if (params.is_null())
00166       {
00167   _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00168   _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00169   _singTol = defaultParams->get<magnitude_type> ("singTol");
00170       }
00171     else
00172       {
00173   try {
00174     _maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00175     if (_maxNumOrthogPasses < 1)
00176       _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00177   } catch (Teuchos::Exceptions::InvalidParameter&) {
00178     _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00179   }
00180 
00181   try {
00182     _blkTol = params->get<magnitude_type> ("blkTol");
00183     if (_blkTol < zero)
00184       _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00185   } catch (Teuchos::Exceptions::InvalidParameter&) {
00186     _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00187   }
00188 
00189   try {
00190     _singTol = params->get<magnitude_type> ("singTol");
00191     if (_singTol < zero)
00192       _singTol = defaultParams->get<magnitude_type> ("singTol");
00193   } catch (Teuchos::Exceptions::InvalidParameter&) {
00194     _singTol = defaultParams->get<magnitude_type> ("singTol");
00195   }
00196       }
00197     maxNumOrthogPasses = _maxNumOrthogPasses;
00198     blkTol = _blkTol;
00199     singTol = _singTol;
00200   }
00201 
00202 
00203   template<class ScalarType, class MV, class OP>
00204   class IMGSOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00205 
00206   private:
00207     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00208     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00209     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00210     typedef MultiVecTraits<ScalarType,MV>      MVT;
00211     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00212 
00213   public:
00214     
00216 
00217 
00218     IMGSOrthoManager( const std::string& label = "Belos",
00219                       Teuchos::RCP<const OP> Op = Teuchos::null,
00220           const int max_ortho_steps = 1,
00221           const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00222           const MagnitudeType sing_tol = 10*MGT::eps() )
00223       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00224   max_ortho_steps_( max_ortho_steps ),
00225   blk_tol_( blk_tol ),
00226   sing_tol_( sing_tol ),
00227         label_( label )
00228     {
00229         std::string orthoLabel = label_ + ": Orthogonalization";
00230 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00231         timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00232 #endif
00233 
00234         std::string updateLabel = label_ + ": Ortho (Update)";
00235 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00236         timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00237 #endif
00238 
00239         std::string normLabel = label_ + ": Ortho (Norm)";
00240 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00241         timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00242 #endif
00243 
00244         std::string ipLabel = label_ + ": Ortho (Inner Product)";
00245 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00246         timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel); 
00247 #endif
00248     };
00249 
00251     ~IMGSOrthoManager() {}
00253 
00254 
00256 
00257 
00259     void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
00260 
00262     void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
00263 
00265     MagnitudeType getBlkTol() const { return blk_tol_; } 
00266 
00268     MagnitudeType getSingTol() const { return sing_tol_; } 
00269 
00271 
00272 
00274 
00275 
00303     void project ( MV &X, Teuchos::RCP<MV> MX, 
00304                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00305                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00306 
00307 
00310     void project ( MV &X, 
00311                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00312                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00313       project(X,Teuchos::null,C,Q);
00314     }
00315 
00316 
00317  
00342     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00343                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00344 
00345 
00348     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00349       return normalize(X,Teuchos::null,B);
00350     }
00351 
00352 
00394     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00395                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00396                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00397                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00398 
00401     int projectAndNormalize ( MV &X, 
00402                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00403                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00404                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const {
00405       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00406     }
00407 
00409 
00411 
00412 
00416     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00417     orthonormError(const MV &X) const {
00418       return orthonormError(X,Teuchos::null);
00419     }
00420 
00425     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00426     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00427 
00431     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00432     orthogError(const MV &X1, const MV &X2) const {
00433       return orthogError(X1,Teuchos::null,X2);
00434     }
00435 
00440     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00441     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00442 
00444 
00446 
00447 
00450     void setLabel(const std::string& label);
00451 
00454     const std::string& getLabel() const { return label_; }
00455 
00457 
00458   private:
00459     
00461     int max_ortho_steps_;
00462     MagnitudeType blk_tol_;
00463     MagnitudeType sing_tol_;
00464 
00466     std::string label_;
00467 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00468     Teuchos::RCP< Teuchos::Time > timerOrtho_, timerUpdate_, 
00469       timerNorm_, timerScale_, timerInnerProd_;
00470 #endif // BELOS_TEUCHOS_TIME_MONITOR
00471   
00473     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00474       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00475       bool completeBasis, int howMany = -1 ) const;
00476     
00478     bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
00479          Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00480          Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00481 
00483     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00484         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00485         Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00486 
00500     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00501            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00502            Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00503            Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;    
00504   };
00505 
00507   // Set the label for this orthogonalization manager and create new timers if it's changed
00508   template<class ScalarType, class MV, class OP>
00509   void IMGSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00510   {   
00511     if (label != label_) {
00512       label_ = label;
00513       std::string orthoLabel = label_ + ": Orthogonalization";
00514 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00515       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00516 #endif
00517 
00518       std::string updateLabel = label_ + ": Ortho (Update)";
00519 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00520       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00521 #endif
00522 
00523       std::string normLabel = label_ + ": Ortho (Norm)";
00524 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00525       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00526 #endif
00527 
00528       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00529 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00530       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel);
00531 #endif
00532     }
00533   } 
00534 
00536   // Compute the distance from orthonormality
00537   template<class ScalarType, class MV, class OP>
00538   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00539   IMGSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00540     const ScalarType ONE = SCT::one();
00541     int rank = MVT::GetNumberVecs(X);
00542     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00543     MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
00544     for (int i=0; i<rank; i++) {
00545       xTx(i,i) -= ONE;
00546     }
00547     return xTx.normFrobenius();
00548   }
00549 
00551   // Compute the distance from orthogonality
00552   template<class ScalarType, class MV, class OP>
00553   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00554   IMGSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00555     int r1 = MVT::GetNumberVecs(X1);
00556     int r2  = MVT::GetNumberVecs(X2);
00557     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00558     MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
00559     return xTx.normFrobenius();
00560   }
00561 
00563   // Find an Op-orthonormal basis for span(X) - span(W)
00564   template<class ScalarType, class MV, class OP>
00565   int IMGSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00566                                     MV &X, Teuchos::RCP<MV> MX, 
00567                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00568                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00569                                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const 
00570   {
00571     using Teuchos::Array;
00572     using Teuchos::null;
00573     using Teuchos::is_null;
00574     using Teuchos::RCP;
00575     using Teuchos::rcp;
00576     using Teuchos::SerialDenseMatrix;
00577     typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
00578     typedef typename Array< RCP< const MV > >::size_type size_type;
00579 
00580 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00581     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00582 #endif
00583 
00584     ScalarType    ONE  = SCT::one();
00585     ScalarType    ZERO  = SCT::zero();
00586 
00587     int nq = Q.size();
00588     int xc = MVT::GetNumberVecs( X );
00589     int xr = MVT::GetVecLength( X );
00590     int rank = xc;
00591 
00592     // If the user doesn't want to store the normalization
00593     // coefficients, allocate some local memory for them.  This will
00594     // go away at the end of this method.
00595     if (is_null (B)) {
00596       B = rcp (new serial_dense_matrix_type (xc, xc));
00597     }
00598     // Likewise, if the user doesn't want to store the projection
00599     // coefficients, allocate some local memory for them.  Also make
00600     // sure that all the entries of C are the right size.  We're going
00601     // to overwrite them anyway, so we don't have to worry about the
00602     // contents (other than to resize them if they are the wrong
00603     // size).
00604     if (C.size() < nq)
00605       C.resize (nq);
00606     for (size_type k = 0; k < nq; ++k)
00607       {
00608   const int numRows = MVT::GetNumberVecs (*Q[k]);
00609   const int numCols = xc; // Number of vectors in X
00610   
00611   if (is_null (C[k]))
00612     C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
00613   else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
00614     {
00615       int err = C[k]->reshape (numRows, numCols);
00616       TEST_FOR_EXCEPTION(err != 0, std::runtime_error, 
00617              "IMGS orthogonalization: failed to reshape "
00618              "C[" << k << "] (the array of block "
00619              "coefficients resulting from projecting X "
00620              "against Q[1:" << nq << "]).");
00621     }
00622       }
00623 
00624     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00625     if (this->_hasOp) {
00626       if (MX == Teuchos::null) {
00627         // we need to allocate space for MX
00628         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00629         OPT::Apply(*(this->_Op),X,*MX);
00630       }
00631     }
00632     else {
00633       // Op == I  -->  MX = X (ignore it if the user passed it in)
00634       MX = Teuchos::rcp( &X, false );
00635     }
00636 
00637     int mxc = MVT::GetNumberVecs( *MX );
00638     int mxr = MVT::GetVecLength( *MX );
00639 
00640     // short-circuit
00641     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
00642 
00643     int numbas = 0;
00644     for (int i=0; i<nq; i++) {
00645       numbas += MVT::GetNumberVecs( *Q[i] );
00646     }
00647 
00648     // check size of B
00649     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00650                         "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00651     // check size of X and MX
00652     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00653                         "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00654     // check size of X w.r.t. MX 
00655     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00656                         "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00657     // check feasibility
00658     //TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00659     //                    "Belos::IMGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00660 
00661     // Some flags for checking dependency returns from the internal orthogonalization methods
00662     bool dep_flg = false;
00663 
00664     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00665     Teuchos::RCP<MV> tmpX, tmpMX;
00666     tmpX = MVT::CloneCopy(X);
00667     if (this->_hasOp) {
00668       tmpMX = MVT::CloneCopy(*MX);
00669     }
00670 
00671     if (xc == 1) {
00672 
00673       // Use the cheaper block orthogonalization.
00674       // NOTE: Don't check for dependencies because the update has one vector.
00675       dep_flg = blkOrtho1( X, MX, C, Q );
00676 
00677       // Normalize the new block X
00678       {
00679 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00680         Teuchos::TimeMonitor normTimer( *timerNorm_ );
00681 #endif
00682         if ( B == Teuchos::null ) {
00683           B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00684         }
00685         std::vector<ScalarType> diag(xc);
00686         MVT::MvDot( X, *MX, diag );
00687         (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
00688         rank = 1; 
00689         MVT::MvAddMv( ONE/(*B)(0,0), X, ZERO, X, X );
00690         if (this->_hasOp) {
00691           // Update MXj.
00692     MVT::MvAddMv( ONE/(*B)(0,0), *MX, ZERO, *MX, *MX );
00693         }
00694       }
00695 
00696     } 
00697     else {
00698 
00699       // Use the cheaper block orthogonalization.
00700       dep_flg = blkOrtho( X, MX, C, Q );
00701 
00702       // If a dependency has been detected in this block, then perform
00703       // the more expensive nonblock (single vector at a time)
00704       // orthogonalization.
00705       if (dep_flg) {
00706         rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00707 
00708         // Copy tmpX back into X.
00709         MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00710         if (this->_hasOp) {
00711     MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00712         }
00713       } 
00714       else {
00715         // There is no dependency, so orthonormalize new block X
00716         rank = findBasis( X, MX, B, false );
00717         if (rank < xc) {
00718     // A dependency was found during orthonormalization of X,
00719     // rerun orthogonalization using more expensive single-
00720     // vector orthogonalization.
00721     rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00722 
00723     // Copy tmpX back into X.
00724     MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00725     if (this->_hasOp) {
00726       MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00727     }
00728         }    
00729       }
00730     } // if (xc == 1) {
00731 
00732     // this should not raise an std::exception; but our post-conditions oblige us to check
00733     TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00734                         "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00735 
00736     // Return the rank of X.
00737     return rank;
00738   }
00739 
00740 
00741 
00743   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00744   template<class ScalarType, class MV, class OP>
00745   int IMGSOrthoManager<ScalarType, MV, OP>::normalize(
00746                                 MV &X, Teuchos::RCP<MV> MX, 
00747                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00748 
00749 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00750     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00751 #endif
00752 
00753     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00754     return findBasis(X, MX, B, true);
00755   }
00756 
00757 
00758 
00760   template<class ScalarType, class MV, class OP>
00761   void IMGSOrthoManager<ScalarType, MV, OP>::project(
00762                           MV &X, Teuchos::RCP<MV> MX, 
00763                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00764                           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00765     // For the inner product defined by the operator Op or the identity (Op == 0)
00766     //   -> Orthogonalize X against each Q[i]
00767     // Modify MX accordingly
00768     //
00769     // Note that when Op is 0, MX is not referenced
00770     //
00771     // Parameter variables
00772     //
00773     // X  : Vectors to be transformed
00774     //
00775     // MX : Image of the block of vectors X by the mass matrix
00776     //
00777     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00778     //
00779     
00780 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00781     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00782 #endif
00783     
00784     int xc = MVT::GetNumberVecs( X );
00785     int xr = MVT::GetVecLength( X );
00786     int nq = Q.size();
00787     std::vector<int> qcs(nq);
00788     // short-circuit
00789     if (nq == 0 || xc == 0 || xr == 0) {
00790       return;
00791     }
00792     int qr = MVT::GetVecLength ( *Q[0] );
00793     // if we don't have enough C, expand it with null references
00794     // if we have too many, resize to throw away the latter ones
00795     // if we have exactly as many as we have Q, this call has no effect
00796     C.resize(nq);
00797 
00798 
00799     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00800     if (this->_hasOp) {
00801       if (MX == Teuchos::null) {
00802         // we need to allocate space for MX
00803         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00804         OPT::Apply(*(this->_Op),X,*MX);
00805       }
00806     }
00807     else {
00808       // Op == I  -->  MX = X (ignore it if the user passed it in)
00809       MX = Teuchos::rcp( &X, false );
00810     }
00811     int mxc = MVT::GetNumberVecs( *MX );
00812     int mxr = MVT::GetVecLength( *MX );
00813 
00814     // check size of X and Q w.r.t. common sense
00815     TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00816                         "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00817     // check size of X w.r.t. MX and Q
00818     TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00819                         "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
00820 
00821     // tally up size of all Q and check/allocate C
00822     int baslen = 0;
00823     for (int i=0; i<nq; i++) {
00824       TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00825                           "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
00826       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00827       TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00828                           "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
00829       baslen += qcs[i];
00830 
00831       // check size of C[i]
00832       if ( C[i] == Teuchos::null ) {
00833         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00834       }
00835       else {
00836         TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00837                            "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
00838       }
00839     }
00840 
00841     // Use the cheaper block orthogonalization, don't check for rank deficiency.
00842     blkOrtho( X, MX, C, Q );
00843 
00844   }  
00845 
00847   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00848   // the rank is numvectors(X)
00849   template<class ScalarType, class MV, class OP>
00850   int IMGSOrthoManager<ScalarType, MV, OP>::findBasis(
00851                   MV &X, Teuchos::RCP<MV> MX, 
00852                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00853                   bool completeBasis, int howMany ) const {
00854     // For the inner product defined by the operator Op or the identity (Op == 0)
00855     //   -> Orthonormalize X 
00856     // Modify MX accordingly
00857     //
00858     // Note that when Op is 0, MX is not referenced
00859     //
00860     // Parameter variables
00861     //
00862     // X  : Vectors to be orthonormalized
00863     //
00864     // MX : Image of the multivector X under the operator Op
00865     //
00866     // Op  : Pointer to the operator for the inner product
00867     //
00868     //
00869 
00870 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00871     Teuchos::TimeMonitor normTimer( *timerNorm_ );
00872 #endif
00873 
00874     const ScalarType ONE  = SCT::one();
00875     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00876 
00877     int xc = MVT::GetNumberVecs( X );
00878     int xr = MVT::GetVecLength( X );
00879 
00880     if (howMany == -1) {
00881       howMany = xc;
00882     }
00883 
00884     /*******************************************************
00885      *  If _hasOp == false, we will not reference MX below *
00886      *******************************************************/
00887 
00888     // if Op==null, MX == X (via pointer)
00889     // Otherwise, either the user passed in MX or we will allocated and compute it
00890     if (this->_hasOp) {
00891       if (MX == Teuchos::null) {
00892         // we need to allocate space for MX
00893         MX = MVT::Clone(X,xc);
00894         OPT::Apply(*(this->_Op),X,*MX);
00895       }
00896     }
00897 
00898     /* if the user doesn't want to store the coefficienets, 
00899      * allocate some local memory for them 
00900      */
00901     if ( B == Teuchos::null ) {
00902       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00903     }
00904 
00905     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00906     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
00907 
00908     // check size of C, B
00909     TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
00910                         "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
00911     TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00912                         "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
00913     TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
00914                         "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
00915     TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
00916                         "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
00917     TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
00918                         "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
00919 
00920     /* xstart is which column we are starting the process with, based on howMany
00921      * columns before xstart are assumed to be Op-orthonormal already
00922      */
00923     int xstart = xc - howMany;
00924 
00925     for (int j = xstart; j < xc; j++) {
00926 
00927       // numX is 
00928       // * number of currently orthonormal columns of X
00929       // * the index of the current column of X
00930       int numX = j;
00931       bool addVec = false;
00932 
00933       // Get a view of the vector currently being worked on.
00934       std::vector<int> index(1);
00935       index[0] = numX;
00936       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
00937       Teuchos::RCP<MV> MXj;
00938       if ((this->_hasOp)) {
00939         // MXj is a view of the current vector in MX
00940         MXj = MVT::CloneViewNonConst( *MX, index );
00941       }
00942       else {
00943         // MXj is a pointer to Xj, and MUST NOT be modified
00944         MXj = Xj;
00945       }
00946 
00947       // Make storage for these Gram-Schmidt iterations.
00948       Teuchos::SerialDenseVector<int,ScalarType> product(numX);
00949       Teuchos::SerialDenseVector<int,ScalarType> P2(1);
00950       Teuchos::RCP<const MV> prevX, prevMX;
00951   
00952       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
00953       //
00954       // Save old MXj vector and compute Op-norm
00955       //
00956       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
00957       MVT::MvDot( *Xj, *MXj, oldDot );
00958       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
00959       TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
00960         "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
00961 
00962       // Perform MGS one vector at a time
00963       for (int ii=0; ii<numX; ii++) {
00964   
00965   index[0] = ii;
00966   prevX = MVT::CloneView( X, index );
00967   if (this->_hasOp) {
00968     prevMX = MVT::CloneView( *MX, index );
00969   }
00970   
00971   for (int i=0; i<max_ortho_steps_; ++i) {
00972     
00973     // product <- prevX^T MXj
00974     {
00975 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00976       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
00977 #endif
00978       MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
00979     }
00980     
00981     // Xj <- Xj - prevX prevX^T MXj   
00982     //     = Xj - prevX product
00983     {
00984 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00985       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00986 #endif
00987       MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00988     }
00989     
00990     // Update MXj
00991     if (this->_hasOp) {
00992       // MXj <- Op*Xj_new
00993       //      = Op*(Xj_old - prevX prevX^T MXj)
00994       //      = MXj - prevMX product
00995 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00996       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
00997 #endif
00998       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00999     }
01000     
01001     // Set coefficients
01002     if ( i==0 )
01003       product[ii] = P2[0];
01004     else
01005       product[ii] += P2[0];
01006     
01007   } // for (int i=0; i<max_ortho_steps_; ++i)
01008 
01009       } // for (int ii=0; ii<numX; ++ii)  
01010   
01011   // Compute Op-norm with old MXj
01012       MVT::MvDot( *Xj, *oldMXj, newDot );
01013       
01014       // Check to see if the new vector is dependent.
01015       if (completeBasis) {
01016   //
01017   // We need a complete basis, so add random vectors if necessary
01018   //
01019   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
01020     
01021     // Add a random vector and orthogonalize it against previous vectors in block.
01022     addVec = true;
01023 #ifdef ORTHO_DEBUG
01024     std::cout << "Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
01025 #endif
01026     //
01027     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01028     Teuchos::RCP<MV> tempMXj;
01029     MVT::MvRandom( *tempXj );
01030     if (this->_hasOp) {
01031       tempMXj = MVT::Clone( X, 1 );
01032       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01033     } 
01034     else {
01035       tempMXj = tempXj;
01036     }
01037     MVT::MvDot( *tempXj, *tempMXj, oldDot );
01038     //
01039     // Perform MGS one vector at a time
01040     for (int ii=0; ii<numX; ii++) {
01041       
01042       index[0] = ii;
01043       prevX = MVT::CloneView( X, index );
01044       if (this->_hasOp) {
01045         prevMX = MVT::CloneView( *MX, index );
01046       }
01047       
01048       for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
01049         {
01050 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01051     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01052 #endif
01053     MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,P2);
01054         }
01055         {
01056 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01057     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01058 #endif
01059     MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
01060         }
01061         if (this->_hasOp) {
01062 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01063     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01064 #endif
01065     MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
01066         }
01067 
01068         // Set coefficients
01069         if ( num_orth==0 )
01070     product[ii] = P2[0];
01071         else
01072     product[ii] += P2[0];       
01073       }
01074     }
01075       
01076     // Compute new Op-norm
01077     MVT::MvDot( *tempXj, *tempMXj, newDot );
01078     //
01079     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01080       // Copy vector into current column of _basisvecs
01081       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
01082       if (this->_hasOp) {
01083         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
01084       }
01085     }
01086     else {
01087       return numX;
01088     } 
01089   }
01090       }
01091       else {
01092   //
01093   // We only need to detect dependencies.
01094   //
01095   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
01096     return numX;
01097   }
01098       }
01099       
01100       
01101       // If we haven't left this method yet, then we can normalize the new vector Xj.
01102       // Normalize Xj.
01103       // Xj <- Xj / std::sqrt(newDot)
01104       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01105       {
01106         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01107         if (this->_hasOp) {
01108     // Update MXj.
01109     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01110         }
01111       }
01112       
01113       // If we've added a random vector, enter a zero in the j'th diagonal element.
01114       if (addVec) {
01115   (*B)(j,j) = ZERO;
01116       }
01117       else {
01118   (*B)(j,j) = diag;
01119       }
01120       
01121       // Save the coefficients, if we are working on the original vector and not a randomly generated one
01122       if (!addVec) {
01123   for (int i=0; i<numX; i++) {
01124     (*B)(i,j) = product(i);
01125   }
01126       }
01127       
01128     } // for (j = 0; j < xc; ++j)
01129     
01130     return xc;
01131   }
01132   
01134   // Routine to compute the block orthogonalization
01135   template<class ScalarType, class MV, class OP>
01136   bool 
01137   IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
01138                 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01139                 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01140   {
01141     int nq = Q.size();
01142     int xc = MVT::GetNumberVecs( X );
01143     const ScalarType ONE  = SCT::one();
01144 
01145     std::vector<int> qcs( nq );
01146     for (int i=0; i<nq; i++) {
01147       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01148     }
01149 
01150     // Perform the Gram-Schmidt transformation for a block of vectors
01151     std::vector<int> index(1);
01152     Teuchos::RCP<const MV> tempQ;
01153 
01154     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01155     // Define the product Q^T * (Op*X)
01156     for (int i=0; i<nq; i++) {
01157 
01158       // Perform MGS one vector at a time
01159       for (int ii=0; ii<qcs[i]; ii++) {
01160   
01161   index[0] = ii;
01162   tempQ = MVT::CloneView( *Q[i], index );
01163   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
01164 
01165   // Multiply Q' with MX
01166   {
01167 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01168     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01169 #endif
01170     MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,X,MX,tempC);
01171   }
01172   // Multiply by Q and subtract the result in X
01173   {
01174 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01175     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01176 #endif
01177     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
01178   }
01179       }
01180       
01181       // Update MX, with the least number of applications of Op as possible
01182       if (this->_hasOp) {
01183   if (xc <= qcs[i]) {
01184     OPT::Apply( *(this->_Op), X, *MX);
01185   }
01186         else {
01187           // this will possibly be used again below; don't delete it
01188           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01189           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01190           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01191         }
01192       }
01193     }
01194     
01195     // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
01196     for (int j = 1; j < max_ortho_steps_; ++j) {
01197       
01198       for (int i=0; i<nq; i++) {
01199 
01200   Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
01201         
01202   // Perform MGS one vector at a time
01203   for (int ii=0; ii<qcs[i]; ii++) {
01204     
01205     index[0] = ii;
01206     tempQ = MVT::CloneView( *Q[i], index );
01207     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
01208     Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
01209     
01210     // Apply another step of modified Gram-Schmidt
01211     {
01212 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01213       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01214 #endif
01215       MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
01216     }
01217     tempC += tempC2;
01218     {
01219 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01220       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01221 #endif
01222       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
01223     }
01224     
01225   }
01226 
01227   // Update MX, with the least number of applications of Op as possible
01228   if (this->_hasOp) {
01229     if (MQ[i].get()) {
01230       // MQ was allocated and computed above; use it
01231       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01232     }
01233     else if (xc <= qcs[i]) {
01234       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01235       OPT::Apply( *(this->_Op), X, *MX);
01236     }
01237   }
01238       } // for (int i=0; i<nq; i++)
01239     } // for (int j = 0; j < max_ortho_steps; ++j)
01240   
01241     return false;
01242   }
01243 
01245   // Routine to compute the block orthogonalization
01246   template<class ScalarType, class MV, class OP>
01247   bool 
01248   IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
01249                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01250                Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01251   {
01252     int nq = Q.size();
01253     int xc = MVT::GetNumberVecs( X );
01254     bool dep_flg = false;
01255     const ScalarType ONE  = SCT::one();
01256 
01257     std::vector<int> qcs( nq );
01258     for (int i=0; i<nq; i++) {
01259       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01260     }
01261 
01262     // Perform the Gram-Schmidt transformation for a block of vectors
01263     
01264     // Compute the initial Op-norms
01265     std::vector<ScalarType> oldDot( xc );
01266     MVT::MvDot( X, *MX, oldDot );
01267 
01268     std::vector<int> index(1);
01269     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01270     Teuchos::RCP<const MV> tempQ;
01271 
01272     // Define the product Q^T * (Op*X)
01273     for (int i=0; i<nq; i++) {
01274 
01275       // Perform MGS one vector at a time
01276       for (int ii=0; ii<qcs[i]; ii++) {
01277   
01278   index[0] = ii;
01279   tempQ = MVT::CloneView( *Q[i], index );
01280   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
01281 
01282   // Multiply Q' with MX
01283   {
01284 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01285     Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01286 #endif
01287     MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC);
01288   }
01289   // Multiply by Q and subtract the result in X
01290   {
01291 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01292     Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01293 #endif
01294     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
01295   }
01296       }
01297 
01298       // Update MX, with the least number of applications of Op as possible
01299       if (this->_hasOp) {
01300         if (xc <= qcs[i]) {
01301           OPT::Apply( *(this->_Op), X, *MX);
01302         }
01303         else {
01304           // this will possibly be used again below; don't delete it
01305           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01306           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01307           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01308         }
01309       }
01310     }
01311 
01312     // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
01313     for (int j = 1; j < max_ortho_steps_; ++j) {
01314       
01315       for (int i=0; i<nq; i++) {
01316   Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
01317         
01318   // Perform MGS one vector at a time
01319   for (int ii=0; ii<qcs[i]; ii++) {
01320     
01321     index[0] = ii;
01322     tempQ = MVT::CloneView( *Q[i], index );
01323     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
01324     Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
01325     
01326     // Apply another step of modified Gram-Schmidt
01327     {
01328 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01329       Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01330 #endif
01331       MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
01332     }
01333     tempC += tempC2;
01334     {
01335 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01336       Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01337 #endif
01338       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
01339     }
01340   }
01341 
01342   // Update MX, with the least number of applications of Op as possible
01343   if (this->_hasOp) {
01344     if (MQ[i].get()) {
01345       // MQ was allocated and computed above; use it
01346       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01347     }
01348     else if (xc <= qcs[i]) {
01349       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01350       OPT::Apply( *(this->_Op), X, *MX);
01351     }
01352   }
01353       } // for (int i=0; i<nq; i++)
01354     } // for (int j = 0; j < max_ortho_steps; ++j)
01355   
01356     // Compute new Op-norms
01357     std::vector<ScalarType> newDot(xc);
01358     MVT::MvDot( X, *MX, newDot );
01359     
01360     // Check to make sure the new block of vectors are not dependent on previous vectors
01361     for (int i=0; i<xc; i++){
01362       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01363   dep_flg = true;
01364   break;
01365       }
01366     } // end for (i=0;...)
01367     
01368     return dep_flg;
01369   }
01370   
01371   template<class ScalarType, class MV, class OP>
01372   int
01373   IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01374                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01375                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01376                    Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
01377   {
01378     Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
01379 
01380     const ScalarType ONE  = SCT::one();
01381     const ScalarType ZERO  = SCT::zero();
01382     
01383     int nq = Q.size();
01384     int xc = MVT::GetNumberVecs( X );
01385     std::vector<int> indX( 1 );
01386     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01387 
01388     std::vector<int> qcs( nq );
01389     for (int i=0; i<nq; i++) {
01390       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01391     }
01392 
01393     // Create pointers for the previous vectors of X that have already been orthonormalized.
01394     Teuchos::RCP<const MV> lastQ;
01395     Teuchos::RCP<MV> Xj, MXj;
01396     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01397 
01398     // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
01399     for (int j=0; j<xc; j++) {
01400       
01401       bool dep_flg = false;
01402       
01403       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01404       if (j > 0) {
01405   std::vector<int> index( j );
01406   for (int ind=0; ind<j; ind++) {
01407     index[ind] = ind;
01408   }
01409   lastQ = MVT::CloneView( X, index );
01410 
01411   // Add these views to the Q and C arrays.
01412   Q.push_back( lastQ );
01413   C.push_back( B );
01414   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01415       }
01416       
01417       // Get a view of the current vector in X to orthogonalize.
01418       indX[0] = j;
01419       Xj = MVT::CloneViewNonConst( X, indX );
01420       if (this->_hasOp) {
01421   MXj = MVT::CloneViewNonConst( *MX, indX );
01422       }
01423       else {
01424   MXj = Xj;
01425       }
01426       
01427       // Compute the initial Op-norms
01428       MVT::MvDot( *Xj, *MXj, oldDot );
01429       
01430       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
01431       Teuchos::RCP<const MV> tempQ;
01432 
01433       // Define the product Q^T * (Op*X)
01434       for (int i=0; i<Q.size(); i++) {
01435 
01436   // Perform MGS one vector at a time
01437   for (int ii=0; ii<qcs[i]; ii++) {
01438     
01439     indX[0] = ii;
01440     tempQ = MVT::CloneView( *Q[i], indX );
01441     // Get a view of the current serial dense matrix
01442     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
01443     
01444     // Multiply Q' with MX
01445     MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,*Xj,MXj,tempC);
01446 
01447     // Multiply by Q and subtract the result in Xj
01448     MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
01449   }
01450  
01451   // Update MXj, with the least number of applications of Op as possible
01452   if (this->_hasOp) {
01453     if (xc <= qcs[i]) {
01454       OPT::Apply( *(this->_Op), *Xj, *MXj);
01455     }
01456     else {
01457       // this will possibly be used again below; don't delete it
01458       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01459       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01460       Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01461       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01462     }
01463   }
01464       }
01465       
01466       // Do any additional steps of modified Gram-Schmidt orthogonalization 
01467       for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
01468   
01469   for (int i=0; i<Q.size(); i++) {
01470     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01471     
01472     // Perform MGS one vector at a time
01473     for (int ii=0; ii<qcs[i]; ii++) {
01474       
01475       indX[0] = ii;
01476       tempQ = MVT::CloneView( *Q[i], indX );
01477       // Get a view of the current serial dense matrix
01478       Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
01479       
01480       // Apply another step of modified Gram-Schmidt
01481       MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *Xj, MXj, tempC2);
01482       MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj ); 
01483     }
01484     
01485     // Add the coefficients into C[i]
01486     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01487     tempC += C2;
01488     
01489     // Update MXj, with the least number of applications of Op as possible
01490     if (this->_hasOp) {
01491       if (MQ[i].get()) {
01492         // MQ was allocated and computed above; use it
01493         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01494       }
01495       else if (xc <= qcs[i]) {
01496         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01497         OPT::Apply( *(this->_Op), *Xj, *MXj);
01498       }
01499     }
01500   } // for (int i=0; i<Q.size(); i++)
01501   
01502       } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
01503       
01504       // Check for linear dependence.
01505       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01506   dep_flg = true;
01507       }
01508       
01509       // Normalize the new vector if it's not dependent
01510       if (!dep_flg) {
01511   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01512   
01513   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01514   if (this->_hasOp) {
01515     // Update MXj.
01516     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01517   }
01518   
01519   // Enter value on diagonal of B.
01520   (*B)(j,j) = diag;
01521       }
01522       else {
01523   // Create a random vector and orthogonalize it against all previous columns of Q.
01524   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01525   Teuchos::RCP<MV> tempMXj;
01526   MVT::MvRandom( *tempXj );
01527   if (this->_hasOp) {
01528     tempMXj = MVT::Clone( X, 1 );
01529     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01530   } 
01531   else {
01532     tempMXj = tempXj;
01533   }
01534   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01535   //
01536   for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
01537     
01538     for (int i=0; i<Q.size(); i++) {
01539       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01540 
01541       // Perform MGS one vector at a time
01542       for (int ii=0; ii<qcs[i]; ii++) {
01543         
01544         indX[0] = ii;
01545         tempQ = MVT::CloneView( *Q[i], indX );
01546         Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
01547       
01548         // Apply another step of modified Gram-Schmidt
01549         MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *tempXj, tempMXj, tempC );
01550         MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
01551         
01552       }
01553 
01554       // Update MXj, with the least number of applications of Op as possible
01555       if (this->_hasOp) {
01556         if (MQ[i].get()) {
01557     // MQ was allocated and computed above; use it
01558     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01559         }
01560         else if (xc <= qcs[i]) {
01561     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01562     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01563         }
01564       }
01565     } // for (int i=0; i<nq; i++)   
01566   } // for (int num_orth=0; num_orth<max_orth_steps_; num_orth++)
01567   
01568   // Compute the Op-norms after the correction step.
01569   MVT::MvDot( *tempXj, *tempMXj, newDot );
01570   
01571   // Copy vector into current column of Xj
01572   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01573     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01574     
01575     // Enter value on diagonal of B.
01576     (*B)(j,j) = ZERO;
01577     
01578     // Copy vector into current column of _basisvecs
01579     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01580     if (this->_hasOp) {
01581       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01582     }
01583   }
01584   else {
01585     return j;
01586   } 
01587       } // if (!dep_flg)
01588       
01589       // Remove the vectors from array
01590       if (j > 0) {
01591   Q.resize( nq );
01592   C.resize( nq );
01593   qcs.resize( nq );
01594       }
01595       
01596     } // for (int j=0; j<xc; j++)
01597     
01598     return xc;
01599   }
01600   
01601 } // namespace Belos
01602 
01603 #endif // BELOS_IMGS_ORTHOMANAGER_HPP
01604 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines