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