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