Belos Version of the Day
BelosICGSOrthoManager.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_ICGS_ORTHOMANAGER_HPP
00048 #define BELOS_ICGS_ORTHOMANAGER_HPP
00049 
00057 // #define ORTHO_DEBUG
00058 
00059 #include "BelosConfigDefs.hpp"
00060 #include "BelosMultiVecTraits.hpp"
00061 #include "BelosOperatorTraits.hpp"
00062 #include "BelosMatOrthoManager.hpp"
00063 
00064 #include "Teuchos_as.hpp"
00065 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp"
00066 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00067 #include "Teuchos_TimeMonitor.hpp"
00068 #endif // BELOS_TEUCHOS_TIME_MONITOR
00069 
00070 namespace Belos {
00071 
00076   template<class ScalarType>
00077   Teuchos::RCP<const Teuchos::ParameterList> TEUCHOS_DEPRECATED
00078   getDefaultIcgsParameters()
00079   {
00080     using Teuchos::as;
00081     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00082     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00083 
00084     // Default parameter values for ICGS orthogonalization.
00085     // Documentation will be embedded in the parameter list.
00086     const int defaultMaxNumOrthogPasses = 2;
00087     const magnitude_type eps = STM::eps();
00088     const magnitude_type defaultBlkTol = 
00089       as<magnitude_type> (10) * STM::squareroot (eps);
00090     const magnitude_type defaultSingTol = as<magnitude_type> (10) * eps;
00091 
00092     Teuchos::RCP<Teuchos::ParameterList> params = 
00093       Teuchos::parameterList ("ICGS");
00094     params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00095                  "Maximum number of orthogonalization passes "
00096                  "(includes the first).  Default is 2, since "
00097                  "\"twice is enough\" for Krylov methods.");
00098     params->set ("blkTol", defaultBlkTol, 
00099                  "Block reorthogonalization threshhold.");
00100     params->set ("singTol", defaultSingTol, 
00101                  "Singular block detection threshold.");
00102     return params;
00103   }
00104 
00109   template<class ScalarType>
00110   Teuchos::RCP<const Teuchos::ParameterList> TEUCHOS_DEPRECATED
00111   getFastIcgsParameters()
00112   {
00113     using Teuchos::ParameterList;
00114     using Teuchos::RCP;
00115     using Teuchos::rcp;
00116     using Teuchos::ScalarTraits;
00117     typedef typename ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00118     typedef ScalarTraits<magnitude_type> STM;
00119 
00120     RCP<const ParameterList> defaultParams = getDefaultIcgsParameters<ScalarType>();
00121     // Start with a clone of the default parameters.
00122     RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
00123 
00124     const int maxBlkOrtho = 1;
00125     const magnitude_type blkTol = STM::zero();
00126     const magnitude_type singTol = STM::zero();
00127 
00128     params->set ("maxNumOrthogPasses", maxBlkOrtho);
00129     params->set ("blkTol", blkTol);
00130     params->set ("singTol", singTol);
00131 
00132     return params;
00133   }
00134 
00144   template<class ScalarType>
00145   void TEUCHOS_DEPRECATED
00146   readIcgsParameters (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00147                       int& maxNumOrthogPasses,
00148                       typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& blkTol,
00149                       typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& singTol)
00150   {
00151     using Teuchos::ParameterList;
00152     using Teuchos::RCP;
00153     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00154     const magnitude_type zero = Teuchos::ScalarTraits<magnitude_type>::zero();
00155     RCP<const ParameterList> defaultParams = getDefaultIcgsParameters<ScalarType>();
00156 
00157     // Using temporary variables and fetching all values before
00158     // setting the output arguments ensures the strong exception
00159     // guarantee for this function: if an exception is thrown, no
00160     // externally visible side effects (in this case, setting the
00161     // output arguments) have taken place.
00162     int _maxNumOrthogPasses;
00163     magnitude_type _blkTol, _singTol;
00164     if (params.is_null()) {
00165       _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00166       _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00167       _singTol = defaultParams->get<magnitude_type> ("singTol");
00168     } else {
00169       try {
00170         _maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00171         if (_maxNumOrthogPasses < 1)
00172           _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00173       } catch (Teuchos::Exceptions::InvalidParameter&) {
00174         _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00175       }
00176 
00177       try {
00178   _blkTol = params->get<magnitude_type> ("blkTol");
00179   if (_blkTol < zero)
00180     _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00181       } catch (Teuchos::Exceptions::InvalidParameter&) {
00182   try {
00183     // People may have used depTol instead of blkTol for this
00184     // parameter's name, by analogy with DGKS.
00185     _blkTol = params->get<magnitude_type> ("depTol");
00186   } catch (Teuchos::Exceptions::InvalidParameter&) {
00187     _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00188   }
00189       }
00190 
00191       try {
00192         _singTol = params->get<magnitude_type> ("singTol");
00193         if (_singTol < zero)
00194           _singTol = defaultParams->get<magnitude_type> ("singTol");
00195       } catch (Teuchos::Exceptions::InvalidParameter&) {
00196         _singTol = defaultParams->get<magnitude_type> ("singTol");
00197       }
00198     }
00199     maxNumOrthogPasses = _maxNumOrthogPasses;
00200     blkTol = _blkTol;
00201     singTol = _singTol;
00202   }
00203 
00204 
00205   template<class ScalarType, class MV, class OP>
00206   class ICGSOrthoManager : 
00207     public MatOrthoManager<ScalarType,MV,OP>, 
00208     public Teuchos::ParameterListAcceptorDefaultBase
00209   {
00210   private:
00211     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00212     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00213     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00214     typedef MultiVecTraits<ScalarType,MV>      MVT;
00215     typedef MultiVecTraitsExt<ScalarType,MV>   MVText;
00216     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00217 
00218   public:
00220 
00221 
00223     ICGSOrthoManager( const std::string& label = "Belos",
00224                       Teuchos::RCP<const OP> Op = Teuchos::null,
00225                       const int max_ortho_steps = 2,
00226                       const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00227                       const MagnitudeType sing_tol = 10*MGT::eps() )
00228       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00229       max_ortho_steps_( max_ortho_steps ),
00230       blk_tol_( blk_tol ),
00231       sing_tol_( sing_tol ),
00232       label_( label )
00233     {
00234 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00235       std::string orthoLabel = label_ + ": Orthogonalization";
00236       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00237 
00238       std::string updateLabel = label_ + ": Ortho (Update)";
00239       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00240 
00241       std::string normLabel = label_ + ": Ortho (Norm)";
00242       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00243 
00244       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00245       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel); 
00246 #endif
00247     }
00248 
00250     ICGSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
00251           const std::string& label = "Belos",
00252                       Teuchos::RCP<const OP> Op = Teuchos::null) :
00253       MatOrthoManager<ScalarType,MV,OP>(Op), 
00254       max_ortho_steps_ (2),
00255       blk_tol_ (10 * MGT::squareroot (MGT::eps())),
00256       sing_tol_ (10 * MGT::eps()),
00257       label_ (label)
00258     {
00259       setParameterList (plist);
00260 
00261 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00262       std::string orthoLabel = label_ + ": Orthogonalization";
00263       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00264 
00265       std::string updateLabel = label_ + ": Ortho (Update)";
00266       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00267 
00268       std::string normLabel = label_ + ": Ortho (Norm)";
00269       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00270 
00271       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00272       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel); 
00273 #endif
00274     }
00275 
00277     ~ICGSOrthoManager() {}
00279 
00281 
00282 
00283     void 
00284     setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00285     {
00286       using Teuchos::Exceptions::InvalidParameterName;
00287       using Teuchos::ParameterList;
00288       using Teuchos::parameterList;
00289       using Teuchos::RCP;
00290 
00291       RCP<const ParameterList> defaultParams = getValidParameters();
00292       RCP<ParameterList> params;
00293       if (plist.is_null()) {
00294   params = parameterList (*defaultParams);
00295       } else {
00296   params = plist;
00297   // Some users might want to specify "blkTol" as "depTol".  Due
00298   // to this case, we don't invoke
00299   // validateParametersAndSetDefaults on params.  Instead, we go
00300   // through the parameter list one parameter at a time and look
00301   // for alternatives.
00302       }
00303   
00304       // Using temporary variables and fetching all values before
00305       // setting the output arguments ensures the strong exception
00306       // guarantee for this function: if an exception is thrown, no
00307       // externally visible side effects (in this case, setting the
00308       // output arguments) have taken place.
00309       int maxNumOrthogPasses;
00310       MagnitudeType blkTol;
00311       MagnitudeType singTol;
00312 
00313       try {
00314   maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00315       } catch (InvalidParameterName&) {
00316   maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00317   params->set ("maxNumOrthogPasses", maxNumOrthogPasses);
00318       }
00319 
00320       // Handling of the "blkTol" parameter is a special case.  This
00321       // is because some users may prefer to call this parameter
00322       // "depTol" for consistency with DGKS.  However, our default
00323       // parameter list calls this "blkTol", and we don't want the
00324       // default list's value to override the user's value.  Thus, we
00325       // first check the user's parameter list for both names, and
00326       // only then access the default parameter list.
00327       try {
00328   blkTol = params->get<MagnitudeType> ("blkTol");
00329       } catch (InvalidParameterName&) {
00330   try {
00331     blkTol = params->get<MagnitudeType> ("depTol");
00332     // "depTol" is the wrong name, so remove it and replace with
00333     // "blkTol".  We'll set "blkTol" below.
00334     params->remove ("depTol");
00335   } catch (InvalidParameterName&) {
00336     blkTol = defaultParams->get<MagnitudeType> ("blkTol");
00337   }
00338   params->set ("blkTol", blkTol);
00339       }
00340 
00341       try {
00342   singTol = params->get<MagnitudeType> ("singTol");
00343       } catch (InvalidParameterName&) {
00344   singTol = defaultParams->get<MagnitudeType> ("singTol");
00345   params->set ("singTol", singTol);
00346       }
00347 
00348       max_ortho_steps_ = maxNumOrthogPasses;
00349       blk_tol_ = blkTol;
00350       sing_tol_ = singTol;
00351 
00352       setMyParamList (params);
00353     }
00354 
00355     Teuchos::RCP<const Teuchos::ParameterList> 
00356     getValidParameters () const
00357     {
00358       using Teuchos::as;
00359       using Teuchos::ParameterList;
00360       using Teuchos::parameterList;
00361       using Teuchos::RCP;
00362 
00363       if (defaultParams_.is_null()) {
00364   RCP<ParameterList> params = parameterList ("ICGS");
00365 
00366   // Default parameter values for ICGS orthogonalization.
00367   // Documentation will be embedded in the parameter list.
00368   const int defaultMaxNumOrthogPasses = 2;
00369   const MagnitudeType eps = MGT::eps();
00370   const MagnitudeType defaultBlkTol = 
00371     as<MagnitudeType> (10) * MGT::squareroot (eps);
00372   const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
00373 
00374   params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00375          "Maximum number of orthogonalization passes (includes the "
00376          "first).  Default is 2, since \"twice is enough\" for Krylov "
00377          "methods.");
00378   params->set ("blkTol", defaultBlkTol, "Block reorthogonalization "
00379          "threshhold.");
00380   params->set ("singTol", defaultSingTol, "Singular block detection "
00381          "threshold.");
00382   defaultParams_ = params;
00383       }
00384       return defaultParams_;
00385     }
00387 
00392     Teuchos::RCP<const Teuchos::ParameterList> 
00393     getFastParameters () const 
00394     {
00395       using Teuchos::as;
00396       using Teuchos::ParameterList;
00397       using Teuchos::parameterList;
00398       using Teuchos::RCP;
00399 
00400       RCP<const ParameterList> defaultParams = getValidParameters ();
00401       // Start with a clone of the default parameters.
00402       RCP<ParameterList> params = parameterList (*defaultParams);
00403 
00404       const int maxBlkOrtho = 1; // No block reorthogonalization
00405       const MagnitudeType blkTol = MGT::zero();
00406       const MagnitudeType singTol = MGT::zero();
00407 
00408       params->set ("maxNumOrthogPasses", maxBlkOrtho);
00409       params->set ("blkTol", blkTol);
00410       params->set ("singTol", singTol);
00411 
00412       return params;
00413     }
00414 
00416 
00417 
00419     void setBlkTol( const MagnitudeType blk_tol ) { 
00420       // Update the parameter list as well.
00421       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00422       if (! params.is_null()) {
00423   // If it's null, then we haven't called setParameterList()
00424   // yet.  It's entirely possible to construct the parameter
00425   // list on demand, so we don't try to create the parameter
00426   // list here.
00427   params->set ("blkTol", blk_tol);
00428       }
00429       blk_tol_ = blk_tol; 
00430     }
00431 
00433     void setSingTol( const MagnitudeType sing_tol ) { 
00434       // Update the parameter list as well.
00435       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00436       if (! params.is_null()) {
00437   // If it's null, then we haven't called setParameterList()
00438   // yet.  It's entirely possible to construct the parameter
00439   // list on demand, so we don't try to create the parameter
00440   // list here.
00441   params->set ("singTol", sing_tol);
00442       }
00443       sing_tol_ = sing_tol; 
00444     }
00445 
00447     MagnitudeType getBlkTol() const { return blk_tol_; } 
00448 
00450     MagnitudeType getSingTol() const { return sing_tol_; } 
00451 
00453 
00454 
00456 
00457 
00485     void project ( MV &X, Teuchos::RCP<MV> MX, 
00486                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00487                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00488 
00489 
00492     void project ( MV &X, 
00493                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00494                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00495       project(X,Teuchos::null,C,Q);
00496     }
00497 
00498 
00499  
00524     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00525                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00526 
00527 
00530     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00531       return normalize(X,Teuchos::null,B);
00532     }
00533 
00534 
00576     int projectAndNormalize ( MV &X, Teuchos::RCP<MV> MX, 
00577                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00578                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00579                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00580 
00583     int projectAndNormalize ( MV &X, 
00584                               Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00585                               Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00586                               Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const {
00587       return projectAndNormalize(X,Teuchos::null,C,B,Q);
00588     }
00589 
00591 
00593 
00594 
00598     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00599     orthonormError(const MV &X) const {
00600       return orthonormError(X,Teuchos::null);
00601     }
00602 
00607     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00608     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00609 
00613     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00614     orthogError(const MV &X1, const MV &X2) const {
00615       return orthogError(X1,Teuchos::null,X2);
00616     }
00617 
00622     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00623     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00624 
00626 
00628 
00629 
00632     void setLabel(const std::string& label);
00633 
00636     const std::string& getLabel() const { return label_; }
00637 
00639 
00640   private:
00642     int max_ortho_steps_;
00644     MagnitudeType blk_tol_;
00646     MagnitudeType sing_tol_;
00647 
00649     std::string label_;
00650 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00651     Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, 
00652       timerNorm_, timerScale_, timerInnerProd_;
00653 #endif // BELOS_TEUCHOS_TIME_MONITOR
00654 
00656     mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
00657   
00659     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00660                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00661                   bool completeBasis, int howMany = -1 ) const;
00662     
00664     bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
00665                      Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00666                      Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00667 
00669     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00670                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00671                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00672 
00686     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00687                        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00688                        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00689                        Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;    
00690   };
00691 
00693   // Set the label for this orthogonalization manager and create new timers if it's changed
00694   template<class ScalarType, class MV, class OP>
00695   void ICGSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00696   { 
00697     if (label != label_) {
00698       label_ = label;
00699       std::string orthoLabel = label_ + ": Orthogonalization";
00700 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00701       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00702 #endif
00703 
00704       std::string updateLabel = label_ + ": Ortho (Update)";
00705 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00706       timerUpdate_ = Teuchos::TimeMonitor::getNewTimer(updateLabel);
00707 #endif
00708 
00709       std::string normLabel = label_ + ": Ortho (Norm)";
00710 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00711       timerNorm_ = Teuchos::TimeMonitor::getNewTimer(normLabel);
00712 #endif
00713 
00714       std::string ipLabel = label_ + ": Ortho (Inner Product)";
00715 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00716       timerInnerProd_ = Teuchos::TimeMonitor::getNewTimer(ipLabel);
00717 #endif
00718     }
00719   } 
00720 
00722   // Compute the distance from orthonormality
00723   template<class ScalarType, class MV, class OP>
00724   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00725   ICGSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00726     const ScalarType ONE = SCT::one();
00727     int rank = MVT::GetNumberVecs(X);
00728     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00729     MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
00730     for (int i=0; i<rank; i++) {
00731       xTx(i,i) -= ONE;
00732     }
00733     return xTx.normFrobenius();
00734   }
00735 
00737   // Compute the distance from orthogonality
00738   template<class ScalarType, class MV, class OP>
00739   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00740   ICGSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00741     int r1 = MVT::GetNumberVecs(X1);
00742     int r2  = MVT::GetNumberVecs(X2);
00743     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00744     MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
00745     return xTx.normFrobenius();
00746   }
00747 
00749   // Find an Op-orthonormal basis for span(X) - span(W)
00750   template<class ScalarType, class MV, class OP>
00751   int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00752                                     MV &X, Teuchos::RCP<MV> MX, 
00753                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00754                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00755                                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const 
00756   {
00757     using Teuchos::Array;
00758     using Teuchos::null;
00759     using Teuchos::is_null;
00760     using Teuchos::RCP;
00761     using Teuchos::rcp;
00762     using Teuchos::SerialDenseMatrix;
00763     typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
00764     typedef typename Array< RCP< const MV > >::size_type size_type;
00765 
00766 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00767     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00768 #endif
00769 
00770     ScalarType    ONE  = SCT::one();
00771     ScalarType    ZERO  = SCT::zero();
00772 
00773     int nq = Q.size();
00774     int xc = MVT::GetNumberVecs( X );
00775     ptrdiff_t xr = MVText::GetGlobalLength( X );
00776     int rank = xc;
00777 
00778     // If the user doesn't want to store the normalization
00779     // coefficients, allocate some local memory for them.  This will
00780     // go away at the end of this method.
00781     if (is_null (B)) {
00782       B = rcp (new serial_dense_matrix_type (xc, xc));
00783     }
00784     // Likewise, if the user doesn't want to store the projection
00785     // coefficients, allocate some local memory for them.  Also make
00786     // sure that all the entries of C are the right size.  We're going
00787     // to overwrite them anyway, so we don't have to worry about the
00788     // contents (other than to resize them if they are the wrong
00789     // size).
00790     if (C.size() < nq)
00791       C.resize (nq);
00792     for (size_type k = 0; k < nq; ++k)
00793       {
00794         const int numRows = MVT::GetNumberVecs (*Q[k]);
00795         const int numCols = xc; // Number of vectors in X
00796 
00797         if (is_null (C[k]))
00798           C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
00799         else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
00800         {
00801           int err = C[k]->reshape (numRows, numCols);
00802           TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, 
00803               "IMGS orthogonalization: failed to reshape "
00804               "C[" << k << "] (the array of block "
00805               "coefficients resulting from projecting X "
00806               "against Q[1:" << nq << "]).");
00807         }
00808       }
00809 
00810     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00811     if (this->_hasOp) {
00812       if (MX == Teuchos::null) {
00813         // we need to allocate space for MX
00814         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00815         OPT::Apply(*(this->_Op),X,*MX);
00816       }
00817     }
00818     else {
00819       // Op == I  -->  MX = X (ignore it if the user passed it in)
00820       MX = Teuchos::rcp( &X, false );
00821     }
00822 
00823     int mxc = MVT::GetNumberVecs( *MX );
00824     ptrdiff_t mxr = MVText::GetGlobalLength( *MX );
00825 
00826     // short-circuit
00827     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::ICGSOrthoManager::projectAndNormalize(): X must be non-empty" );
00828 
00829     int numbas = 0;
00830     for (int i=0; i<nq; i++) {
00831       numbas += MVT::GetNumberVecs( *Q[i] );
00832     }
00833 
00834     // check size of B
00835     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00836                         "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00837     // check size of X and MX
00838     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00839                         "Belos::ICGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00840     // check size of X w.r.t. MX 
00841     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00842                         "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00843     // check feasibility
00844     //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00845     //                    "Belos::ICGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00846 
00847     // Some flags for checking dependency returns from the internal orthogonalization methods
00848     bool dep_flg = false;
00849 
00850     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00851     Teuchos::RCP<MV> tmpX, tmpMX;
00852     tmpX = MVT::CloneCopy(X);
00853     if (this->_hasOp) {
00854       tmpMX = MVT::CloneCopy(*MX);
00855     }
00856 
00857     if (xc == 1) {
00858 
00859       // Use the cheaper block orthogonalization.
00860       // NOTE: Don't check for dependencies because the update has one vector.
00861       dep_flg = blkOrtho1( X, MX, C, Q );
00862 
00863       // Normalize the new block X
00864       {
00865 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00866         Teuchos::TimeMonitor normTimer( *timerNorm_ );
00867 #endif
00868         if ( B == Teuchos::null ) {
00869           B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00870         }
00871         std::vector<ScalarType> diag(xc);
00872         MVT::MvDot( X, *MX, diag );
00873         (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
00874         rank = 1; 
00875         MVT::MvAddMv( ONE/(*B)(0,0), X, ZERO, X, X );
00876         if (this->_hasOp) {
00877           // Update MXj.
00878           MVT::MvAddMv( ONE/(*B)(0,0), *MX, ZERO, *MX, *MX );
00879         }
00880       }
00881 
00882     } 
00883     else {
00884 
00885       // Use the cheaper block orthogonalization.
00886       dep_flg = blkOrtho( X, MX, C, Q );
00887 
00888       // If a dependency has been detected in this block, then perform
00889       // the more expensive nonblock (single vector at a time)
00890       // orthogonalization.
00891       if (dep_flg) {
00892         rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00893 
00894         // Copy tmpX back into X.
00895         MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00896         if (this->_hasOp) {
00897           MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00898         }
00899       } 
00900       else {
00901         // There is no dependency, so orthonormalize new block X
00902         rank = findBasis( X, MX, B, false );
00903         if (rank < xc) {
00904           // A dependency was found during orthonormalization of X,
00905           // rerun orthogonalization using more expensive nonblock
00906           // orthogonalization.
00907           rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00908 
00909           // Copy tmpX back into X.
00910           MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00911           if (this->_hasOp) {
00912             MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00913           }
00914         }    
00915       }
00916     } // if (xc == 1) {
00917 
00918     // this should not raise an std::exception; but our post-conditions oblige us to check
00919     TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00920                         "Belos::ICGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00921 
00922     // Return the rank of X.
00923     return rank;
00924   }
00925 
00926 
00927 
00929   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00930   template<class ScalarType, class MV, class OP>
00931   int ICGSOrthoManager<ScalarType, MV, OP>::normalize(
00932                                 MV &X, Teuchos::RCP<MV> MX, 
00933                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00934 
00935 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00936     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00937 #endif
00938 
00939     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00940     return findBasis(X, MX, B, true);
00941   }
00942 
00943 
00944 
00946   template<class ScalarType, class MV, class OP>
00947   void ICGSOrthoManager<ScalarType, MV, OP>::project(
00948                           MV &X, Teuchos::RCP<MV> MX, 
00949                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00950                           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00951     // For the inner product defined by the operator Op or the identity (Op == 0)
00952     //   -> Orthogonalize X against each Q[i]
00953     // Modify MX accordingly
00954     //
00955     // Note that when Op is 0, MX is not referenced
00956     //
00957     // Parameter variables
00958     //
00959     // X  : Vectors to be transformed
00960     //
00961     // MX : Image of the block of vectors X by the mass matrix
00962     //
00963     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00964     //
00965     
00966 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00967     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00968 #endif
00969     
00970     int xc = MVT::GetNumberVecs( X );
00971     ptrdiff_t xr = MVText::GetGlobalLength( X );
00972     int nq = Q.size();
00973     std::vector<int> qcs(nq);
00974     // short-circuit
00975     if (nq == 0 || xc == 0 || xr == 0) {
00976       return;
00977     }
00978     ptrdiff_t qr = MVText::GetGlobalLength ( *Q[0] );
00979     // if we don't have enough C, expand it with null references
00980     // if we have too many, resize to throw away the latter ones
00981     // if we have exactly as many as we have Q, this call has no effect
00982     C.resize(nq);
00983 
00984 
00985     /******   DO NOT MODIFY *MX IF _hasOp == false   ******/
00986     if (this->_hasOp) {
00987       if (MX == Teuchos::null) {
00988         // we need to allocate space for MX
00989         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00990         OPT::Apply(*(this->_Op),X,*MX);
00991       }
00992     }
00993     else {
00994       // Op == I  -->  MX = X (ignore it if the user passed it in)
00995       MX = Teuchos::rcp( &X, false );
00996     }
00997     int mxc = MVT::GetNumberVecs( *MX );
00998     ptrdiff_t mxr = MVText::GetGlobalLength( *MX );
00999 
01000     // check size of X and Q w.r.t. common sense
01001     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
01002                         "Belos::ICGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
01003     // check size of X w.r.t. MX and Q
01004     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
01005                         "Belos::ICGSOrthoManager::project(): Size of X not consistant with MX,Q" );
01006 
01007     // tally up size of all Q and check/allocate C
01008     int baslen = 0;
01009     for (int i=0; i<nq; i++) {
01010       TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument, 
01011                           "Belos::ICGSOrthoManager::project(): Q lengths not mutually consistant" );
01012       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01013       TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
01014                           "Belos::ICGSOrthoManager::project(): Q has less rows than columns" );
01015       baslen += qcs[i];
01016 
01017       // check size of C[i]
01018       if ( C[i] == Teuchos::null ) {
01019         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
01020       }
01021       else {
01022         TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
01023                            "Belos::ICGSOrthoManager::project(): Size of Q not consistant with size of C" );
01024       }
01025     }
01026 
01027     // Use the cheaper block orthogonalization, don't check for rank deficiency.
01028     blkOrtho( X, MX, C, Q );
01029 
01030   }  
01031 
01033   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
01034   // the rank is numvectors(X)
01035   template<class ScalarType, class MV, class OP>
01036   int ICGSOrthoManager<ScalarType, MV, OP>::findBasis(
01037                 MV &X, Teuchos::RCP<MV> MX, 
01038                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
01039                 bool completeBasis, int howMany ) const {
01040     // For the inner product defined by the operator Op or the identity (Op == 0)
01041     //   -> Orthonormalize X 
01042     // Modify MX accordingly
01043     //
01044     // Note that when Op is 0, MX is not referenced
01045     //
01046     // Parameter variables
01047     //
01048     // X  : Vectors to be orthonormalized
01049     //
01050     // MX : Image of the multivector X under the operator Op
01051     //
01052     // Op  : Pointer to the operator for the inner product
01053     //
01054     //
01055 
01056 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01057     Teuchos::TimeMonitor normTimer( *timerNorm_ );
01058 #endif
01059 
01060     const ScalarType ONE  = SCT::one();
01061     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
01062 
01063     int xc = MVT::GetNumberVecs( X );
01064     ptrdiff_t xr = MVText::GetGlobalLength( X );
01065 
01066     if (howMany == -1) {
01067       howMany = xc;
01068     }
01069 
01070     /*******************************************************
01071      *  If _hasOp == false, we will not reference MX below *
01072      *******************************************************/
01073 
01074     // if Op==null, MX == X (via pointer)
01075     // Otherwise, either the user passed in MX or we will allocated and compute it
01076     if (this->_hasOp) {
01077       if (MX == Teuchos::null) {
01078         // we need to allocate space for MX
01079         MX = MVT::Clone(X,xc);
01080         OPT::Apply(*(this->_Op),X,*MX);
01081       }
01082     }
01083 
01084     /* if the user doesn't want to store the coefficienets, 
01085      * allocate some local memory for them 
01086      */
01087     if ( B == Teuchos::null ) {
01088       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
01089     }
01090 
01091     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
01092     ptrdiff_t mxr = (this->_hasOp) ? MVText::GetGlobalLength( *MX ) : xr;
01093 
01094     // check size of C, B
01095     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
01096                         "Belos::ICGSOrthoManager::findBasis(): X must be non-empty" );
01097     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
01098                         "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
01099     TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
01100                         "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
01101     TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument, 
01102                         "Belos::ICGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
01103     TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
01104                         "Belos::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
01105 
01106     /* xstart is which column we are starting the process with, based on howMany
01107      * columns before xstart are assumed to be Op-orthonormal already
01108      */
01109     int xstart = xc - howMany;
01110 
01111     for (int j = xstart; j < xc; j++) {
01112 
01113       // numX is 
01114       // * number of currently orthonormal columns of X
01115       // * the index of the current column of X
01116       int numX = j;
01117       bool addVec = false;
01118 
01119       // Get a view of the vector currently being worked on.
01120       std::vector<int> index(1);
01121       index[0] = numX;
01122       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
01123       Teuchos::RCP<MV> MXj;
01124       if ((this->_hasOp)) {
01125         // MXj is a view of the current vector in MX
01126         MXj = MVT::CloneViewNonConst( *MX, index );
01127       }
01128       else {
01129         // MXj is a pointer to Xj, and MUST NOT be modified
01130         MXj = Xj;
01131       }
01132 
01133       // Get a view of the previous vectors.
01134       std::vector<int> prev_idx( numX );
01135       Teuchos::RCP<const MV> prevX, prevMX;
01136 
01137       if (numX > 0) {
01138         for (int i=0; i<numX; i++) {
01139           prev_idx[i] = i;
01140         }
01141         prevX = MVT::CloneView( X, prev_idx );
01142         if (this->_hasOp) {
01143           prevMX = MVT::CloneView( *MX, prev_idx );
01144         }
01145       } 
01146 
01147       // Make storage for these Gram-Schmidt iterations.
01148       Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
01149       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01150       //
01151       // Save old MXj vector and compute Op-norm
01152       //
01153       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
01154       MVT::MvDot( *Xj, *MXj, oldDot );
01155       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
01156       TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
01157           "Belos::ICGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
01158 
01159       if (numX > 0) {
01160 
01161         Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
01162 
01163         for (int i=0; i<max_ortho_steps_; ++i) {
01164 
01165           // product <- prevX^T MXj
01166           {
01167 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01168             Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01169 #endif
01170             MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
01171           }
01172 
01173           // Xj <- Xj - prevX prevX^T MXj   
01174           //     = Xj - prevX product
01175           {
01176 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01177             Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01178 #endif
01179             MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
01180           }
01181 
01182           // Update MXj
01183           if (this->_hasOp) {
01184             // MXj <- Op*Xj_new
01185             //      = Op*(Xj_old - prevX prevX^T MXj)
01186             //      = MXj - prevMX product
01187 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01188             Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01189 #endif
01190             MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
01191           }
01192 
01193           // Set coefficients
01194           if ( i==0 )
01195             product = P2;
01196           else
01197             product += P2;
01198         }
01199 
01200       } // if (numX > 0)
01201 
01202       // Compute Op-norm with old MXj
01203       MVT::MvDot( *Xj, *oldMXj, newDot );
01204 
01205       // Check to see if the new vector is dependent.
01206       if (completeBasis) {
01207         //
01208         // We need a complete basis, so add random vectors if necessary
01209         //
01210         if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
01211 
01212           // Add a random vector and orthogonalize it against previous vectors in block.
01213           addVec = true;
01214 #ifdef ORTHO_DEBUG
01215           std::cout << "Belos::ICGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
01216 #endif
01217           //
01218           Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01219           Teuchos::RCP<MV> tempMXj;
01220           MVT::MvRandom( *tempXj );
01221           if (this->_hasOp) {
01222             tempMXj = MVT::Clone( X, 1 );
01223             OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01224           } 
01225           else {
01226             tempMXj = tempXj;
01227           }
01228           MVT::MvDot( *tempXj, *tempMXj, oldDot );
01229           //
01230           for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
01231             {
01232 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01233               Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01234 #endif
01235               MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,product);
01236             }
01237             {
01238 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01239               Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01240 #endif
01241               MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
01242             }
01243             if (this->_hasOp) {
01244 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01245               Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01246 #endif
01247               MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
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       // If we haven't left this method yet, then we can normalize the new vector Xj.
01275       // Normalize Xj.
01276       // Xj <- Xj / std::sqrt(newDot)
01277       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01278       {
01279         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01280         if (this->_hasOp) {
01281           // Update MXj.
01282           MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01283         }
01284       }
01285 
01286       // If we've added a random vector, enter a zero in the j'th diagonal element.
01287       if (addVec) {
01288         (*B)(j,j) = ZERO;
01289       }
01290       else {
01291         (*B)(j,j) = diag;
01292       }
01293 
01294       // Save the coefficients, if we are working on the original vector and not a randomly generated one
01295       if (!addVec) {
01296         for (int i=0; i<numX; i++) {
01297           (*B)(i,j) = product(i,0);
01298         }
01299       }
01300 
01301     } // for (j = 0; j < xc; ++j)
01302 
01303     return xc;
01304   }
01305 
01307   // Routine to compute the block orthogonalization
01308   template<class ScalarType, class MV, class OP>
01309   bool 
01310   ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX, 
01311                                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01312                                                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01313   {
01314     int nq = Q.size();
01315     int xc = MVT::GetNumberVecs( X );
01316     const ScalarType ONE  = SCT::one();
01317 
01318     std::vector<int> qcs( nq );
01319     for (int i=0; i<nq; i++) {
01320       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01321     }
01322 
01323     // Perform the Gram-Schmidt transformation for a block of vectors
01324 
01325     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01326     // Define the product Q^T * (Op*X)
01327     for (int i=0; i<nq; i++) {
01328       // Multiply Q' with MX
01329       {
01330 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01331         Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01332 #endif
01333         MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
01334       }
01335       // Multiply by Q and subtract the result in X
01336       {
01337 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01338         Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01339 #endif
01340         MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
01341       }
01342 
01343       // Update MX, with the least number of applications of Op as possible
01344       if (this->_hasOp) {
01345         if (xc <= qcs[i]) {
01346           OPT::Apply( *(this->_Op), X, *MX);
01347         }
01348         else {
01349           // this will possibly be used again below; don't delete it
01350           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01351           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01352           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01353         }
01354       }
01355     }
01356 
01357     // Do as many steps of classical Gram-Schmidt as required by max_ortho_steps_
01358     for (int j = 1; j < max_ortho_steps_; ++j) {
01359       
01360       for (int i=0; i<nq; i++) {
01361         Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
01362 
01363         // Apply another step of classical Gram-Schmidt
01364         {
01365 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01366           Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01367 #endif
01368           MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
01369         }
01370         *C[i] += C2;
01371         {
01372 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01373           Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01374 #endif
01375           MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
01376         }
01377 
01378         // Update MX, with the least number of applications of Op as possible
01379         if (this->_hasOp) {
01380           if (MQ[i].get()) {
01381             // MQ was allocated and computed above; use it
01382             MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01383           }
01384           else if (xc <= qcs[i]) {
01385             // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01386             OPT::Apply( *(this->_Op), X, *MX);
01387           }
01388         }
01389       } // for (int i=0; i<nq; i++)
01390     } // for (int j = 0; j < max_ortho_steps; ++j)
01391 
01392     return false;
01393   }
01394 
01396   // Routine to compute the block orthogonalization
01397   template<class ScalarType, class MV, class OP>
01398   bool 
01399   ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
01400                                                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01401                                                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01402   {
01403     int nq = Q.size();
01404     int xc = MVT::GetNumberVecs( X );
01405     bool dep_flg = false;
01406     const ScalarType ONE  = SCT::one();
01407 
01408     std::vector<int> qcs( nq );
01409     for (int i=0; i<nq; i++) {
01410       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01411     }
01412 
01413     // Perform the Gram-Schmidt transformation for a block of vectors
01414     
01415     // Compute the initial Op-norms
01416     std::vector<ScalarType> oldDot( xc );
01417     MVT::MvDot( X, *MX, oldDot );
01418 
01419     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01420     // Define the product Q^T * (Op*X)
01421     for (int i=0; i<nq; i++) {
01422       // Multiply Q' with MX
01423       {
01424 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01425         Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01426 #endif
01427         MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
01428       }
01429       // Multiply by Q and subtract the result in X
01430       {
01431 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01432         Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01433 #endif
01434         MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
01435       }
01436       // Update MX, with the least number of applications of Op as possible
01437       if (this->_hasOp) {
01438         if (xc <= qcs[i]) {
01439           OPT::Apply( *(this->_Op), X, *MX);
01440         }
01441         else {
01442           // this will possibly be used again below; don't delete it
01443           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01444           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01445           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01446         }
01447       }
01448     }
01449 
01450     // Do as many steps of classical Gram-Schmidt as required by max_ortho_steps_
01451     for (int j = 1; j < max_ortho_steps_; ++j) {
01452       
01453       for (int i=0; i<nq; i++) {
01454         Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
01455 
01456         // Apply another step of classical Gram-Schmidt
01457         {
01458 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01459           Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
01460 #endif
01461           MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
01462         }
01463         *C[i] += C2;
01464         {
01465 #ifdef BELOS_TEUCHOS_TIME_MONITOR
01466           Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
01467 #endif
01468           MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
01469         }
01470 
01471         // Update MX, with the least number of applications of Op as possible
01472         if (this->_hasOp) {
01473           if (MQ[i].get()) {
01474             // MQ was allocated and computed above; use it
01475             MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01476           }
01477           else if (xc <= qcs[i]) {
01478             // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01479             OPT::Apply( *(this->_Op), X, *MX);
01480           }
01481         }
01482       } // for (int i=0; i<nq; i++)
01483     } // for (int j = 0; j < max_ortho_steps; ++j)
01484 
01485     // Compute new Op-norms
01486     std::vector<ScalarType> newDot(xc);
01487     MVT::MvDot( X, *MX, newDot );
01488 
01489     // Check to make sure the new block of vectors are not dependent on previous vectors
01490     for (int i=0; i<xc; i++){
01491       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01492         dep_flg = true;
01493         break;
01494       }
01495     } // end for (i=0;...)
01496 
01497     return dep_flg;
01498   }
01499   
01500   template<class ScalarType, class MV, class OP>
01501   int
01502   ICGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01503                                                        Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01504                                                        Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01505                                                        Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
01506   {
01507     Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
01508 
01509     const ScalarType ONE  = SCT::one();
01510     const ScalarType ZERO  = SCT::zero();
01511     
01512     int nq = Q.size();
01513     int xc = MVT::GetNumberVecs( X );
01514     std::vector<int> indX( 1 );
01515     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01516 
01517     std::vector<int> qcs( nq );
01518     for (int i=0; i<nq; i++) {
01519       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01520     }
01521 
01522     // Create pointers for the previous vectors of X that have already been orthonormalized.
01523     Teuchos::RCP<const MV> lastQ;
01524     Teuchos::RCP<MV> Xj, MXj;
01525     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01526 
01527     // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
01528     for (int j=0; j<xc; j++) {
01529       
01530       bool dep_flg = false;
01531       
01532       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01533       if (j > 0) {
01534         std::vector<int> index( j );
01535         for (int ind=0; ind<j; ind++) {
01536           index[ind] = ind;
01537         }
01538         lastQ = MVT::CloneView( X, index );
01539 
01540         // Add these views to the Q and C arrays.
01541         Q.push_back( lastQ );
01542         C.push_back( B );
01543         qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01544       }
01545 
01546       // Get a view of the current vector in X to orthogonalize.
01547       indX[0] = j;
01548       Xj = MVT::CloneViewNonConst( X, indX );
01549       if (this->_hasOp) {
01550         MXj = MVT::CloneViewNonConst( *MX, indX );
01551       }
01552       else {
01553         MXj = Xj;
01554       }
01555 
01556       // Compute the initial Op-norms
01557       MVT::MvDot( *Xj, *MXj, oldDot );
01558 
01559       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
01560       // Define the product Q^T * (Op*X)
01561       for (int i=0; i<Q.size(); i++) {
01562 
01563         // Get a view of the current serial dense matrix
01564         Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01565 
01566         // Multiply Q' with MX
01567         MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,tempC);
01568         // Multiply by Q and subtract the result in Xj
01569         MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
01570 
01571         // Update MXj, with the least number of applications of Op as possible
01572         if (this->_hasOp) {
01573           if (xc <= qcs[i]) {
01574             OPT::Apply( *(this->_Op), *Xj, *MXj);
01575           }
01576           else {
01577             // this will possibly be used again below; don't delete it
01578             MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01579             OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01580             MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01581           }
01582         }
01583       }
01584 
01585       // Do any additional steps of classical Gram-Schmidt orthogonalization 
01586       for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
01587 
01588         for (int i=0; i<Q.size(); i++) {
01589           Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01590           Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01591 
01592           // Apply another step of classical Gram-Schmidt
01593           MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,C2);
01594           tempC += C2;
01595           MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
01596 
01597           // Update MXj, with the least number of applications of Op as possible
01598           if (this->_hasOp) {
01599             if (MQ[i].get()) {
01600               // MQ was allocated and computed above; use it
01601               MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01602             }
01603             else if (xc <= qcs[i]) {
01604               // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01605               OPT::Apply( *(this->_Op), *Xj, *MXj);
01606             }
01607           }
01608         } // for (int i=0; i<Q.size(); i++)
01609 
01610       } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
01611 
01612       // Check for linear dependence.
01613       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01614         dep_flg = true;
01615       }
01616 
01617       // Normalize the new vector if it's not dependent
01618       if (!dep_flg) {
01619         ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01620 
01621         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01622         if (this->_hasOp) {
01623           // Update MXj.
01624           MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01625         }
01626 
01627         // Enter value on diagonal of B.
01628         (*B)(j,j) = diag;
01629       }
01630       else {
01631         // Create a random vector and orthogonalize it against all previous columns of Q.
01632         Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01633         Teuchos::RCP<MV> tempMXj;
01634         MVT::MvRandom( *tempXj );
01635         if (this->_hasOp) {
01636           tempMXj = MVT::Clone( X, 1 );
01637           OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01638         } 
01639         else {
01640           tempMXj = tempXj;
01641         }
01642         MVT::MvDot( *tempXj, *tempMXj, oldDot );
01643         //
01644         for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
01645 
01646           for (int i=0; i<Q.size(); i++) {
01647             Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01648 
01649             // Apply another step of classical Gram-Schmidt
01650             MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*tempXj,tempMXj,product);
01651             MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
01652 
01653             // Update MXj, with the least number of applications of Op as possible
01654             if (this->_hasOp) {
01655               if (MQ[i].get()) {
01656                 // MQ was allocated and computed above; use it
01657                 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01658               }
01659               else if (xc <= qcs[i]) {
01660                 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01661                 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01662               }
01663             }
01664           } // for (int i=0; i<nq; i++)
01665 
01666         }
01667 
01668         // Compute the Op-norms after the correction step.
01669         MVT::MvDot( *tempXj, *tempMXj, newDot );
01670 
01671         // Copy vector into current column of Xj
01672         if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01673           ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01674 
01675           // Enter value on diagonal of B.
01676           (*B)(j,j) = ZERO;
01677 
01678           // Copy vector into current column of _basisvecs
01679           MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01680           if (this->_hasOp) {
01681             MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01682           }
01683         }
01684         else {
01685           return j;
01686         } 
01687       } // if (!dep_flg)
01688 
01689       // Remove the vectors from array
01690       if (j > 0) {
01691         Q.resize( nq );
01692         C.resize( nq );
01693         qcs.resize( nq );
01694       }
01695 
01696     } // for (int j=0; j<xc; j++)
01697 
01698     return xc;
01699   }
01700 
01701 } // namespace Belos
01702 
01703 #endif // BELOS_ICGS_ORTHOMANAGER_HPP
01704 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines