Belos Version of the Day
BelosDGKSOrthoManager.hpp
Go to the documentation of this file.
00001 //@HEADER
00002 // ************************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                  Copyright 2004 Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 
00047 #ifndef BELOS_DGKS_ORTHOMANAGER_HPP
00048 #define BELOS_DGKS_ORTHOMANAGER_HPP
00049 
00057 // #define ORTHO_DEBUG
00058 
00059 #include "BelosConfigDefs.hpp"
00060 #include "BelosMultiVecTraits.hpp"
00061 #include "BelosOperatorTraits.hpp"
00062 #include "BelosMatOrthoManager.hpp"
00063 
00064 #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   getDefaultDgksParameters()
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 DGKS orthogonalization.
00085     // Documentation will be embedded in the parameter list.
00086     const int defaultMaxNumOrthogPasses = 2;
00087     const magnitude_type one = STM::one();
00088     const magnitude_type eps = STM::eps();
00089     const magnitude_type squareRootTwo = 
00090       STM::squareroot (as<magnitude_type> (2));
00091     const magnitude_type defaultBlkTol = 
00092       as<magnitude_type> (10) * STM::squareroot (eps);
00093     const magnitude_type defaultDepTol = one / squareRootTwo;
00094     const magnitude_type defaultSingTol = as<magnitude_type> (10) * eps;
00095 
00096     Teuchos::RCP<Teuchos::ParameterList> params = 
00097       Teuchos::parameterList ("DGKS");
00098     params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00099      "Maximum number of orthogonalization passes "
00100      "(includes the first).  Default is 2, since "
00101      "\"twice is enough\" for Krylov methods.");
00102     params->set ("blkTol", defaultBlkTol, 
00103      "Block reorthogonalization threshhold.");
00104     params->set ("depTol", defaultDepTol, 
00105      "(Non-block) reorthogonalization threshold.");
00106     params->set ("singTol", defaultSingTol, 
00107      "Singular block detection threshold.");
00108     return params;
00109   }
00110 
00116   template<class ScalarType>
00117   Teuchos::RCP<const Teuchos::ParameterList> TEUCHOS_DEPRECATED
00118   getFastDgksParameters()
00119   {
00120     using Teuchos::ParameterList;
00121     using Teuchos::RCP;
00122     using Teuchos::rcp;
00123     using Teuchos::ScalarTraits;
00124     typedef typename ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00125     typedef ScalarTraits<magnitude_type> STM;
00126 
00127     RCP<const ParameterList> defaultParams = getDefaultDgksParameters<ScalarType> ();
00128     // Start with a clone of the default parameters.
00129     RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
00130 
00131     const int maxBlkOrtho = 1;
00132     const magnitude_type blkTol = STM::zero();
00133     const magnitude_type depTol = STM::zero();
00134     const magnitude_type singTol = STM::zero();
00135 
00136     params->set ("maxNumOrthogPasses", maxBlkOrtho);
00137     params->set ("blkTol", blkTol);
00138     params->set ("depTol", depTol);
00139     params->set ("singTol", singTol);
00140 
00141     return params;
00142   }
00143 
00153   template<class ScalarType>
00154   void TEUCHOS_DEPRECATED
00155   readDgksParameters (const Teuchos::RCP<const Teuchos::ParameterList>& params,
00156           int& maxNumOrthogPasses,
00157           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& blkTol,
00158           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& depTol,
00159           typename Teuchos::ScalarTraits<ScalarType>::magnitudeType& singTol)
00160   {
00161     using Teuchos::ParameterList;
00162     using Teuchos::RCP;
00163     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType magnitude_type;
00164     const magnitude_type zero = Teuchos::ScalarTraits<magnitude_type>::zero();
00165     RCP<const ParameterList> defaultParams = getDefaultDgksParameters<ScalarType>();
00166 
00167     // Using temporary variables and fetching all values before
00168     // setting the output arguments ensures the strong exception
00169     // guarantee for this function: if an exception is thrown, no
00170     // externally visible side effects (in this case, setting the
00171     // output arguments) have taken place.
00172     int _maxNumOrthogPasses;
00173     magnitude_type _blkTol, _depTol, _singTol;
00174     if (params.is_null())
00175       {
00176   _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00177   _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00178   _depTol = defaultParams->get<magnitude_type> ("depTol");  
00179   _singTol = defaultParams->get<magnitude_type> ("singTol");
00180       }
00181     else
00182       {
00183   try {
00184     _maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00185     if (_maxNumOrthogPasses < 1)
00186       _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00187   } catch (Teuchos::Exceptions::InvalidParameter&) {
00188     _maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
00189   }
00190 
00191   try {
00192     _blkTol = params->get<magnitude_type> ("blkTol");
00193     if (_blkTol < zero)
00194       _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00195   } catch (Teuchos::Exceptions::InvalidParameter&) {
00196     _blkTol = defaultParams->get<magnitude_type> ("blkTol");
00197   }
00198 
00199   try {
00200     _depTol = params->get<magnitude_type> ("depTol");
00201     if (_depTol < zero)
00202       _depTol = defaultParams->get<magnitude_type> ("depTol");
00203   } catch (Teuchos::Exceptions::InvalidParameter&) {
00204     _depTol = defaultParams->get<magnitude_type> ("depTol");
00205   }
00206 
00207   try {
00208     _singTol = params->get<magnitude_type> ("singTol");
00209     if (_singTol < zero)
00210       _singTol = defaultParams->get<magnitude_type> ("singTol");
00211   } catch (Teuchos::Exceptions::InvalidParameter&) {
00212     _singTol = defaultParams->get<magnitude_type> ("singTol");
00213   }
00214       }
00215     maxNumOrthogPasses = _maxNumOrthogPasses;
00216     blkTol = _blkTol;
00217     depTol = _depTol;
00218     singTol = _singTol;
00219   }
00220 
00221   template<class ScalarType, class MV, class OP>
00222   class DGKSOrthoManager : 
00223     public MatOrthoManager<ScalarType,MV,OP>,
00224     public Teuchos::ParameterListAcceptorDefaultBase
00225   {
00226   private:
00227     typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00228     typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
00229     typedef Teuchos::ScalarTraits<ScalarType>  SCT;
00230     typedef MultiVecTraits<ScalarType,MV>      MVT;
00231     typedef OperatorTraits<ScalarType,MV,OP>   OPT;
00232 
00233   public:
00235 
00236 
00238     DGKSOrthoManager( const std::string& label = "Belos",
00239                       Teuchos::RCP<const OP> Op = Teuchos::null,
00240           const int max_blk_ortho = 2,
00241           const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
00242           const MagnitudeType dep_tol = MGT::one()/MGT::squareroot( 2*MGT::one() ),
00243           const MagnitudeType sing_tol = 10*MGT::eps() )
00244       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00245   max_blk_ortho_( max_blk_ortho ),
00246   blk_tol_( blk_tol ),
00247   dep_tol_( dep_tol ),
00248   sing_tol_( sing_tol ),
00249   label_( label )
00250     {
00251 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00252       std::string orthoLabel = label_ + ": Orthogonalization";
00253       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer( orthoLabel );
00254 #endif
00255     }    
00256 
00258     DGKSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
00259           const std::string& label = "Belos",
00260                       Teuchos::RCP<const OP> Op = Teuchos::null)
00261       : MatOrthoManager<ScalarType,MV,OP>(Op), 
00262   max_blk_ortho_ (2),
00263   blk_tol_ (Teuchos::as<MagnitudeType>(10) * MGT::squareroot(MGT::eps())),
00264   dep_tol_ (MGT::one() / MGT::squareroot (Teuchos::as<MagnitudeType>(2))),
00265   sing_tol_ (Teuchos::as<MagnitudeType>(10) * MGT::eps()),
00266   label_( label )
00267     {
00268       setParameterList (plist);
00269 
00270 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00271       std::string orthoLabel = label_ + ": Orthogonalization";
00272       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer( orthoLabel );
00273 #endif
00274     }    
00275 
00277     ~DGKSOrthoManager() {}
00279 
00281 
00282 
00283     void 
00284     setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00285     {
00286       using Teuchos::ParameterList;
00287       using Teuchos::parameterList;
00288       using Teuchos::RCP;
00289 
00290       RCP<const ParameterList> defaultParams = getValidParameters();
00291       RCP<ParameterList> params;
00292       if (plist.is_null()) {
00293   // No need to validate default parameters.
00294   params = parameterList (*defaultParams);
00295       } else {
00296   params = plist;
00297   params->validateParametersAndSetDefaults (*defaultParams);
00298       }
00299   
00300       // Using temporary variables and fetching all values before
00301       // setting the output arguments ensures the strong exception
00302       // guarantee for this function: if an exception is thrown, no
00303       // externally visible side effects (in this case, setting the
00304       // output arguments) have taken place.
00305       const int maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
00306       const MagnitudeType blkTol = params->get<MagnitudeType> ("blkTol");
00307       const MagnitudeType depTol = params->get<MagnitudeType> ("depTol");
00308       const MagnitudeType singTol = params->get<MagnitudeType> ("singTol");
00309 
00310       max_blk_ortho_ = maxNumOrthogPasses;
00311       blk_tol_ = blkTol;
00312       dep_tol_ = depTol;
00313       sing_tol_ = singTol;
00314 
00315       setMyParamList (params);
00316     }
00317 
00318     Teuchos::RCP<const Teuchos::ParameterList>
00319     getValidParameters () const
00320     {
00321       using Teuchos::as;
00322       using Teuchos::ParameterList;
00323       using Teuchos::parameterList;
00324       using Teuchos::RCP;
00325 
00326       if (defaultParams_.is_null()) {
00327   RCP<ParameterList> params = parameterList ("DGKS");
00328   const MagnitudeType eps = MGT::eps ();
00329 
00330   // Default parameter values for DGKS orthogonalization.
00331   // Documentation will be embedded in the parameter list.
00332   const int defaultMaxNumOrthogPasses = 2;
00333   const MagnitudeType defaultBlkTol = 
00334     as<MagnitudeType> (10) * MGT::squareroot (eps);
00335   const MagnitudeType defaultDepTol =
00336     MGT::one() / MGT::squareroot (as<MagnitudeType> (2));
00337 
00338   const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
00339 
00340   params->set ("maxNumOrthogPasses", defaultMaxNumOrthogPasses,
00341          "Maximum number of orthogonalization passes (includes the "
00342          "first).  Default is 2, since \"twice is enough\" for Krylov "
00343          "methods.");
00344   params->set ("blkTol", defaultBlkTol, "Block reorthogonalization "
00345          "threshhold.");
00346   params->set ("depTol", defaultDepTol, 
00347          "(Non-block) reorthogonalization threshold.");
00348   params->set ("singTol", defaultSingTol, "Singular block detection "
00349          "threshold.");
00350   defaultParams_ = params;
00351       }
00352       return defaultParams_;
00353     }
00354 
00356 
00357     Teuchos::RCP<const Teuchos::ParameterList>
00358     getFastParameters() const
00359     {
00360       using Teuchos::ParameterList;
00361       using Teuchos::RCP;
00362       using Teuchos::rcp;
00363 
00364       RCP<const ParameterList> defaultParams = getValidParameters ();
00365       // Start with a clone of the default parameters.
00366       RCP<ParameterList> params = rcp (new ParameterList (*defaultParams));
00367 
00368       const int maxBlkOrtho = 1;
00369       const MagnitudeType blkTol = MGT::zero();
00370       const MagnitudeType depTol = MGT::zero();
00371       const MagnitudeType singTol = MGT::zero();
00372 
00373       params->set ("maxNumOrthogPasses", maxBlkOrtho);
00374       params->set ("blkTol", blkTol);
00375       params->set ("depTol", depTol);
00376       params->set ("singTol", singTol);
00377 
00378       return params;
00379     }
00380 
00382 
00383 
00385     void setBlkTol( const MagnitudeType blk_tol ) { 
00386       // Update the parameter list as well.
00387       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00388       if (! params.is_null()) {
00389   // If it's null, then we haven't called setParameterList()
00390   // yet.  It's entirely possible to construct the parameter
00391   // list on demand, so we don't try to create the parameter
00392   // list here.
00393   params->set ("blkTol", blk_tol);
00394       }
00395       blk_tol_ = blk_tol; 
00396     }
00397 
00399     void setDepTol( const MagnitudeType dep_tol ) { 
00400       // Update the parameter list as well.
00401       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00402       if (! params.is_null()) {
00403   params->set ("depTol", dep_tol);
00404       }
00405       dep_tol_ = dep_tol; 
00406     }
00407 
00409     void setSingTol( const MagnitudeType sing_tol ) { 
00410       // Update the parameter list as well.
00411       Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
00412       if (! params.is_null()) {
00413   params->set ("singTol", sing_tol);
00414       }
00415       sing_tol_ = sing_tol; 
00416     }
00417 
00419     MagnitudeType getBlkTol() const { return blk_tol_; } 
00420 
00422     MagnitudeType getDepTol() const { return dep_tol_; } 
00423 
00425     MagnitudeType getSingTol() const { return sing_tol_; } 
00426 
00428 
00429 
00431 
00432 
00460     void project ( MV &X, Teuchos::RCP<MV> MX, 
00461                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00462                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00463 
00464 
00467     void project ( MV &X, 
00468                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00469                    Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00470       project(X,Teuchos::null,C,Q);
00471     }
00472 
00473 
00474  
00499     int normalize ( MV &X, Teuchos::RCP<MV> MX, 
00500                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
00501 
00502 
00505     int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00506       return normalize(X,Teuchos::null,B);
00507     }
00508 
00509 
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 
00589     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00590     orthonormError(const MV &X) const {
00591       return orthonormError(X,Teuchos::null);
00592     }
00593 
00600     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00601     orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
00602 
00606     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00607     orthogError(const MV &X1, const MV &X2) const {
00608       return orthogError(X1,Teuchos::null,X2);
00609     }
00610 
00615     typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00616     orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
00617 
00619 
00621 
00622 
00625     void setLabel(const std::string& label);
00626 
00629     const std::string& getLabel() const { return label_; }
00630 
00632 
00633   private:
00634     
00636     int max_blk_ortho_;
00638     MagnitudeType blk_tol_;
00640     MagnitudeType dep_tol_;
00642     MagnitudeType sing_tol_;
00643 
00645     std::string label_;
00646 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00647     Teuchos::RCP<Teuchos::Time> timerOrtho_;
00648 #endif // BELOS_TEUCHOS_TIME_MONITOR
00649 
00651     mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
00652 
00654     int findBasis(MV &X, Teuchos::RCP<MV> MX, 
00655       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C, 
00656       bool completeBasis, int howMany = -1 ) const;
00657     
00659     bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
00660         Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00661         Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
00662 
00676     int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
00677            Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00678            Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00679            Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;    
00680   };
00681 
00683   // Set the label for this orthogonalization manager and create new timers if it's changed
00684   template<class ScalarType, class MV, class OP>
00685   void DGKSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
00686   {
00687     if (label != label_) {
00688       label_ = label;
00689       std::string orthoLabel = label_ + ": Orthogonalization";
00690 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00691       timerOrtho_ = Teuchos::TimeMonitor::getNewTimer(orthoLabel);
00692 #endif
00693     }
00694   }
00695 
00697   // Compute the distance from orthonormality
00698   template<class ScalarType, class MV, class OP>
00699   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00700   DGKSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
00701     const ScalarType ONE = SCT::one();
00702     int rank = MVT::GetNumberVecs(X);
00703     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00704     MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
00705     for (int i=0; i<rank; i++) {
00706       xTx(i,i) -= ONE;
00707     }
00708     return xTx.normFrobenius();
00709   }
00710 
00712   // Compute the distance from orthogonality
00713   template<class ScalarType, class MV, class OP>
00714   typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 
00715   DGKSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
00716     int r1 = MVT::GetNumberVecs(X1);
00717     int r2  = MVT::GetNumberVecs(X2);
00718     Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
00719     MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
00720     return xTx.normFrobenius();
00721   }
00722 
00724   // Find an Op-orthonormal basis for span(X) - span(W)
00725   template<class ScalarType, class MV, class OP>
00726   int DGKSOrthoManager<ScalarType, MV, OP>::projectAndNormalize(
00727                                     MV &X, Teuchos::RCP<MV> MX, 
00728                                     Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00729                                     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
00730                                     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q ) const 
00731   {
00732     using Teuchos::Array;
00733     using Teuchos::null;
00734     using Teuchos::is_null;
00735     using Teuchos::RCP;
00736     using Teuchos::rcp;
00737     using Teuchos::SerialDenseMatrix;
00738     typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
00739     typedef typename Array< RCP< const MV > >::size_type size_type;
00740     
00741 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00742     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00743 #endif
00744 
00745     ScalarType    ONE  = SCT::one();
00746     ScalarType    ZERO  = SCT::zero();
00747 
00748     int nq = Q.size();
00749     int xc = MVT::GetNumberVecs( X );
00750     int xr = MVT::GetVecLength( X );
00751     int rank = xc;
00752 
00753     // If the user doesn't want to store the normalization
00754     // coefficients, allocate some local memory for them.  This will
00755     // go away at the end of this method.
00756     if (is_null (B)) {
00757       B = rcp (new serial_dense_matrix_type (xc, xc));
00758     }
00759     // Likewise, if the user doesn't want to store the projection
00760     // coefficients, allocate some local memory for them.  Also make
00761     // sure that all the entries of C are the right size.  We're going
00762     // to overwrite them anyway, so we don't have to worry about the
00763     // contents (other than to resize them if they are the wrong
00764     // size).
00765     if (C.size() < nq)
00766       C.resize (nq);
00767     for (size_type k = 0; k < nq; ++k)
00768       {
00769   const int numRows = MVT::GetNumberVecs (*Q[k]);
00770   const int numCols = xc; // Number of vectors in X
00771   
00772   if (is_null (C[k]))
00773     C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
00774   else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
00775     {
00776       int err = C[k]->reshape (numRows, numCols);
00777       TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error, 
00778              "DGKS orthogonalization: failed to reshape "
00779              "C[" << k << "] (the array of block "
00780              "coefficients resulting from projecting X "
00781              "against Q[1:" << nq << "]).");
00782     }
00783       }
00784 
00785     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00786     if (this->_hasOp) {
00787       if (MX == Teuchos::null) {
00788         // we need to allocate space for MX
00789         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00790         OPT::Apply(*(this->_Op),X,*MX);
00791       }
00792     }
00793     else {
00794       // Op == I  -->  MX = X (ignore it if the user passed it in)
00795       MX = Teuchos::rcp( &X, false );
00796     }
00797 
00798     int mxc = MVT::GetNumberVecs( *MX );
00799     int mxr = MVT::GetVecLength( *MX );
00800 
00801     // short-circuit
00802     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
00803 
00804     int numbas = 0;
00805     for (int i=0; i<nq; i++) {
00806       numbas += MVT::GetNumberVecs( *Q[i] );
00807     }
00808 
00809     // check size of B
00810     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
00811                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
00812     // check size of X and MX
00813     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00814                         "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
00815     // check size of X w.r.t. MX 
00816     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 
00817                         "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
00818     // check feasibility
00819     //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 
00820     //                    "Belos::DGKSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
00821 
00822     // Some flags for checking dependency returns from the internal orthogonalization methods
00823     bool dep_flg = false;
00824 
00825     // Make a temporary copy of X and MX, just in case a block dependency is detected.
00826     Teuchos::RCP<MV> tmpX, tmpMX;
00827     tmpX = MVT::CloneCopy(X);
00828     if (this->_hasOp) {
00829       tmpMX = MVT::CloneCopy(*MX);
00830     }
00831 
00832     // Use the cheaper block orthogonalization.
00833     dep_flg = blkOrtho( X, MX, C, Q );
00834 
00835     // If a dependency has been detected in this block, then perform
00836     // the more expensive single-vector orthogonalization.
00837     if (dep_flg) {
00838       rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00839 
00840       // Copy tmpX back into X.
00841       MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00842       if (this->_hasOp) {
00843   MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00844       }
00845     } 
00846     else {
00847       // There is no dependency, so orthonormalize new block X
00848       rank = findBasis( X, MX, B, false );
00849       if (rank < xc) {
00850   // A dependency was found during orthonormalization of X,
00851   // rerun orthogonalization using more expensive single-vector orthogonalization.
00852   rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
00853 
00854   // Copy tmpX back into X.
00855   MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
00856   if (this->_hasOp) {
00857     MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
00858   }
00859       }    
00860     }
00861 
00862     // this should not raise an std::exception; but our post-conditions oblige us to check
00863     TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 
00864                         "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
00865 
00866     // Return the rank of X.
00867     return rank;
00868   }
00869 
00870 
00871 
00873   // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
00874   template<class ScalarType, class MV, class OP>
00875   int DGKSOrthoManager<ScalarType, MV, OP>::normalize(
00876                                 MV &X, Teuchos::RCP<MV> MX, 
00877                                 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
00878 
00879 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00880     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00881 #endif
00882 
00883     // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
00884     return findBasis(X, MX, B, true);
00885   }
00886 
00887 
00888 
00890   template<class ScalarType, class MV, class OP>
00891   void DGKSOrthoManager<ScalarType, MV, OP>::project(
00892                           MV &X, Teuchos::RCP<MV> MX, 
00893                           Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
00894                           Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
00895     // For the inner product defined by the operator Op or the identity (Op == 0)
00896     //   -> Orthogonalize X against each Q[i]
00897     // Modify MX accordingly
00898     //
00899     // Note that when Op is 0, MX is not referenced
00900     //
00901     // Parameter variables
00902     //
00903     // X  : Vectors to be transformed
00904     //
00905     // MX : Image of the block vector X by the mass matrix
00906     //
00907     // Q  : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
00908     //
00909 
00910 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00911     Teuchos::TimeMonitor orthotimer(*timerOrtho_);
00912 #endif
00913 
00914     int xc = MVT::GetNumberVecs( X );
00915     int xr = MVT::GetVecLength( X );
00916     int nq = Q.size();
00917     std::vector<int> qcs(nq);
00918     // short-circuit
00919     if (nq == 0 || xc == 0 || xr == 0) {
00920       return;
00921     }
00922     int qr = MVT::GetVecLength ( *Q[0] );
00923     // if we don't have enough C, expand it with null references
00924     // if we have too many, resize to throw away the latter ones
00925     // if we have exactly as many as we have Q, this call has no effect
00926     C.resize(nq);
00927 
00928 
00929     /******   DO NO MODIFY *MX IF _hasOp == false   ******/
00930     if (this->_hasOp) {
00931       if (MX == Teuchos::null) {
00932         // we need to allocate space for MX
00933         MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00934         OPT::Apply(*(this->_Op),X,*MX);
00935       }
00936     }
00937     else {
00938       // Op == I  -->  MX = X (ignore it if the user passed it in)
00939       MX = Teuchos::rcp( &X, false );
00940     }
00941     int mxc = MVT::GetNumberVecs( *MX );
00942     int mxr = MVT::GetVecLength( *MX );
00943 
00944     // check size of X and Q w.r.t. common sense
00945     TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 
00946                         "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
00947     // check size of X w.r.t. MX and Q
00948     TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 
00949                         "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
00950 
00951     // tally up size of all Q and check/allocate C
00952     int baslen = 0;
00953     for (int i=0; i<nq; i++) {
00954       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 
00955                           "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
00956       qcs[i] = MVT::GetNumberVecs( *Q[i] );
00957       TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 
00958                           "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
00959       baslen += qcs[i];
00960 
00961       // check size of C[i]
00962       if ( C[i] == Teuchos::null ) {
00963         C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00964       }
00965       else {
00966         TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 
00967                            "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
00968       }
00969     }
00970 
00971     // Use the cheaper block orthogonalization, don't check for rank deficiency.
00972     blkOrtho( X, MX, C, Q );
00973 
00974   }  
00975 
00977   // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 
00978   // the rank is numvectors(X)
00979   template<class ScalarType, class MV, class OP>
00980   int DGKSOrthoManager<ScalarType, MV, OP>::findBasis(
00981                   MV &X, Teuchos::RCP<MV> MX, 
00982                   Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00983                   bool completeBasis, int howMany ) const {
00984     // For the inner product defined by the operator Op or the identity (Op == 0)
00985     //   -> Orthonormalize X 
00986     // Modify MX accordingly
00987     //
00988     // Note that when Op is 0, MX is not referenced
00989     //
00990     // Parameter variables
00991     //
00992     // X  : Vectors to be orthonormalized
00993     //
00994     // MX : Image of the multivector X under the operator Op
00995     //
00996     // Op  : Pointer to the operator for the inner product
00997     //
00998     //
00999 
01000     const ScalarType ONE  = SCT::one();
01001     const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
01002 
01003     int xc = MVT::GetNumberVecs( X );
01004     int xr = MVT::GetVecLength( X );
01005 
01006     if (howMany == -1) {
01007       howMany = xc;
01008     }
01009 
01010     /*******************************************************
01011      *  If _hasOp == false, we will not reference MX below *
01012      *******************************************************/
01013 
01014     // if Op==null, MX == X (via pointer)
01015     // Otherwise, either the user passed in MX or we will allocated and compute it
01016     if (this->_hasOp) {
01017       if (MX == Teuchos::null) {
01018         // we need to allocate space for MX
01019         MX = MVT::Clone(X,xc);
01020         OPT::Apply(*(this->_Op),X,*MX);
01021       }
01022     }
01023 
01024     /* if the user doesn't want to store the coefficienets, 
01025      * allocate some local memory for them 
01026      */
01027     if ( B == Teuchos::null ) {
01028       B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
01029     }
01030 
01031     int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
01032     int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX )  : xr;
01033 
01034     // check size of C, B
01035     TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 
01036                         "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
01037     TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 
01038                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
01039     TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 
01040                         "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
01041     TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 
01042                         "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
01043     TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument, 
01044                         "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
01045 
01046     /* xstart is which column we are starting the process with, based on howMany
01047      * columns before xstart are assumed to be Op-orthonormal already
01048      */
01049     int xstart = xc - howMany;
01050 
01051     for (int j = xstart; j < xc; j++) {
01052 
01053       // numX is 
01054       // * number of currently orthonormal columns of X
01055       // * the index of the current column of X
01056       int numX = j;
01057       bool addVec = false;
01058 
01059       // Get a view of the vector currently being worked on.
01060       std::vector<int> index(1);
01061       index[0] = numX;
01062       Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
01063       Teuchos::RCP<MV> MXj;
01064       if ((this->_hasOp)) {
01065         // MXj is a view of the current vector in MX
01066         MXj = MVT::CloneViewNonConst( *MX, index );
01067       }
01068       else {
01069         // MXj is a pointer to Xj, and MUST NOT be modified
01070         MXj = Xj;
01071       }
01072 
01073       // Get a view of the previous vectors.
01074       std::vector<int> prev_idx( numX );
01075       Teuchos::RCP<const MV> prevX, prevMX;
01076 
01077       if (numX > 0) {
01078         for (int i=0; i<numX; i++) {
01079           prev_idx[i] = i;
01080         }
01081         prevX = MVT::CloneView( X, prev_idx );
01082         if (this->_hasOp) {
01083           prevMX = MVT::CloneView( *MX, prev_idx );
01084         }
01085       } 
01086 
01087       // Make storage for these Gram-Schmidt iterations.
01088       Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
01089       std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01090       //
01091       // Save old MXj vector and compute Op-norm
01092       //
01093       Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 
01094       MVT::MvDot( *Xj, *MXj, oldDot );
01095       // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
01096       TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError, 
01097         "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
01098 
01099       if (numX > 0) {
01100   // Apply the first step of Gram-Schmidt
01101   
01102   // product <- prevX^T MXj
01103   MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,product);
01104   
01105   // Xj <- Xj - prevX prevX^T MXj   
01106   //     = Xj - prevX product
01107   MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
01108   
01109   // Update MXj
01110   if (this->_hasOp) {
01111     // MXj <- Op*Xj_new
01112     //      = Op*(Xj_old - prevX prevX^T MXj)
01113     //      = MXj - prevMX product
01114     MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
01115   }
01116   
01117   // Compute new Op-norm
01118   MVT::MvDot( *Xj, *MXj, newDot );
01119   
01120   // Check if a correction is needed.
01121   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(dep_tol_*oldDot[0]) ) {
01122     // Apply the second step of Gram-Schmidt
01123     // This is the same as above
01124     Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
01125     
01126     MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
01127     product += P2;
01128     MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
01129     if ((this->_hasOp)) {
01130       MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
01131     }
01132   } // if (newDot[0] < dep_tol_*oldDot[0])
01133   
01134       } // if (numX > 0)
01135 
01136         // Compute Op-norm with old MXj
01137       MVT::MvDot( *Xj, *oldMXj, newDot );
01138 
01139       // Check to see if the new vector is dependent.
01140       if (completeBasis) {
01141   //
01142   // We need a complete basis, so add random vectors if necessary
01143   //
01144   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
01145     
01146     // Add a random vector and orthogonalize it against previous vectors in block.
01147     addVec = true;
01148 #ifdef ORTHO_DEBUG
01149     std::cout << "Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
01150 #endif
01151     //
01152     Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01153     Teuchos::RCP<MV> tempMXj;
01154     MVT::MvRandom( *tempXj );
01155     if (this->_hasOp) {
01156       tempMXj = MVT::Clone( X, 1 );
01157       OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01158     } 
01159     else {
01160       tempMXj = tempXj;
01161     }
01162     MVT::MvDot( *tempXj, *tempMXj, oldDot );
01163     //
01164     for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
01165       MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,product);
01166       MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
01167       if (this->_hasOp) {
01168         MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
01169       }
01170     }
01171     // Compute new Op-norm
01172     MVT::MvDot( *tempXj, *tempMXj, newDot );
01173     //
01174     if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01175       // Copy vector into current column of _basisvecs
01176       MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
01177       if (this->_hasOp) {
01178         MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
01179       }
01180     }
01181     else {
01182       return numX;
01183     } 
01184   }
01185       }
01186       else {
01187   //
01188   // We only need to detect dependencies.
01189   //
01190   if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
01191     return numX;
01192   }
01193       }
01194       
01195       // If we haven't left this method yet, then we can normalize the new vector Xj.
01196       // Normalize Xj.
01197       // Xj <- Xj / std::sqrt(newDot)
01198       ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01199 
01200       if (SCT::magnitude(diag) > ZERO) {      
01201         MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01202         if (this->_hasOp) {
01203     // Update MXj.
01204     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01205         }
01206       }
01207 
01208       // If we've added a random vector, enter a zero in the j'th diagonal element.
01209       if (addVec) {
01210   (*B)(j,j) = ZERO;
01211       }
01212       else {
01213   (*B)(j,j) = diag;
01214       }
01215 
01216       // Save the coefficients, if we are working on the original vector and not a randomly generated one
01217       if (!addVec) {
01218   for (int i=0; i<numX; i++) {
01219     (*B)(i,j) = product(i,0);
01220   }
01221       }
01222 
01223     } // for (j = 0; j < xc; ++j)
01224 
01225     return xc;
01226   }
01227 
01229   // Routine to compute the block orthogonalization
01230   template<class ScalarType, class MV, class OP>
01231   bool 
01232   DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX, 
01233                Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01234                Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
01235   {
01236     int nq = Q.size();
01237     int xc = MVT::GetNumberVecs( X );
01238     bool dep_flg = false;
01239     const ScalarType ONE  = SCT::one();
01240 
01241     std::vector<int> qcs( nq );
01242     for (int i=0; i<nq; i++) {
01243       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01244     }
01245 
01246     // Perform the Gram-Schmidt transformation for a block of vectors
01247     
01248     // Compute the initial Op-norms
01249     std::vector<ScalarType> oldDot( xc );
01250     MVT::MvDot( X, *MX, oldDot );
01251 
01252     Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
01253     // Define the product Q^T * (Op*X)
01254     for (int i=0; i<nq; i++) {
01255       // Multiply Q' with MX
01256       MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
01257       // Multiply by Q and subtract the result in X
01258       MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
01259 
01260       // Update MX, with the least number of applications of Op as possible
01261       if (this->_hasOp) {
01262         if (xc <= qcs[i]) {
01263           OPT::Apply( *(this->_Op), X, *MX);
01264         }
01265         else {
01266           // this will possibly be used again below; don't delete it
01267           MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01268           OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01269           MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
01270         }
01271       }
01272     }
01273 
01274     // Do as many steps of classical Gram-Schmidt as required by max_blk_ortho_
01275     for (int j = 1; j < max_blk_ortho_; ++j) {
01276       
01277       for (int i=0; i<nq; i++) {
01278   Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
01279         
01280   // Apply another step of classical Gram-Schmidt
01281   MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
01282   *C[i] += C2;
01283   MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
01284         
01285   // Update MX, with the least number of applications of Op as possible
01286   if (this->_hasOp) {
01287     if (MQ[i].get()) {
01288       // MQ was allocated and computed above; use it
01289       MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
01290     }
01291     else if (xc <= qcs[i]) {
01292       // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01293       OPT::Apply( *(this->_Op), X, *MX);
01294     }
01295   }
01296       } // for (int i=0; i<nq; i++)
01297     } // for (int j = 0; j < max_blk_ortho; ++j)
01298   
01299     // Compute new Op-norms
01300     std::vector<ScalarType> newDot(xc);
01301     MVT::MvDot( X, *MX, newDot );
01302  
01303     // Check to make sure the new block of vectors are not dependent on previous vectors
01304     for (int i=0; i<xc; i++){
01305       if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
01306   dep_flg = true;
01307   break;
01308       }
01309     } // end for (i=0;...)
01310 
01311     return dep_flg;
01312   }
01313   
01314 
01315   template<class ScalarType, class MV, class OP>
01316   int
01317   DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX, 
01318                    Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 
01319                    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 
01320                    Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
01321   {
01322     Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
01323 
01324     const ScalarType ONE  = SCT::one();
01325     const ScalarType ZERO  = SCT::zero();
01326     
01327     int nq = Q.size();
01328     int xc = MVT::GetNumberVecs( X );
01329     std::vector<int> indX( 1 );
01330     std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
01331 
01332     std::vector<int> qcs( nq );
01333     for (int i=0; i<nq; i++) {
01334       qcs[i] = MVT::GetNumberVecs( *Q[i] );
01335     }
01336 
01337     // Create pointers for the previous vectors of X that have already been orthonormalized.
01338     Teuchos::RCP<const MV> lastQ;
01339     Teuchos::RCP<MV> Xj, MXj;
01340     Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
01341 
01342     // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
01343     for (int j=0; j<xc; j++) {
01344       
01345       bool dep_flg = false;
01346       
01347       // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
01348       if (j > 0) {
01349   std::vector<int> index( j );
01350   for (int ind=0; ind<j; ind++) {
01351     index[ind] = ind;
01352   }
01353   lastQ = MVT::CloneView( X, index );
01354 
01355   // Add these views to the Q and C arrays.
01356   Q.push_back( lastQ );
01357   C.push_back( B );
01358   qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
01359       }
01360       
01361       // Get a view of the current vector in X to orthogonalize.
01362       indX[0] = j;
01363       Xj = MVT::CloneViewNonConst( X, indX );
01364       if (this->_hasOp) {
01365   MXj = MVT::CloneViewNonConst( *MX, indX );
01366       }
01367       else {
01368   MXj = Xj;
01369       }
01370       
01371       // Compute the initial Op-norms
01372       MVT::MvDot( *Xj, *MXj, oldDot );
01373       
01374       Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
01375       // Define the product Q^T * (Op*X)
01376       for (int i=0; i<Q.size(); i++) {
01377 
01378   // Get a view of the current serial dense matrix
01379   Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01380 
01381   // Multiply Q' with MX
01382   MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,tempC);
01383   // Multiply by Q and subtract the result in Xj
01384   MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
01385   
01386   // Update MXj, with the least number of applications of Op as possible
01387   if (this->_hasOp) {
01388     if (xc <= qcs[i]) {
01389       OPT::Apply( *(this->_Op), *Xj, *MXj);
01390     }
01391     else {
01392       // this will possibly be used again below; don't delete it
01393       MQ[i] = MVT::Clone( *Q[i], qcs[i] );
01394       OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
01395       MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
01396     }
01397   }
01398       }
01399       
01400       // Compute the Op-norms
01401       MVT::MvDot( *Xj, *MXj, newDot );
01402       
01403       // Do one step of classical Gram-Schmidt orthogonalization 
01404       // with a second correction step if needed.
01405       
01406       if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
01407   
01408   for (int i=0; i<Q.size(); i++) {
01409     Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
01410     Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
01411     
01412     // Apply another step of classical Gram-Schmidt
01413     MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,C2);
01414     tempC += C2;
01415     MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
01416     
01417     // Update MXj, with the least number of applications of Op as possible
01418     if (this->_hasOp) {
01419       if (MQ[i].get()) {
01420         // MQ was allocated and computed above; use it
01421         MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
01422       }
01423       else if (xc <= qcs[i]) {
01424         // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01425         OPT::Apply( *(this->_Op), *Xj, *MXj);
01426       }
01427     }
01428   } // for (int i=0; i<Q.size(); i++)
01429   
01430   // Compute the Op-norms after the correction step.
01431   MVT::MvDot( *Xj, *MXj, newDot );
01432   
01433       } // if ()
01434       
01435       // Check for linear dependence.
01436       if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
01437   dep_flg = true;
01438       }
01439       
01440       // Normalize the new vector if it's not dependent
01441       if (!dep_flg) {
01442   ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01443   
01444   MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
01445   if (this->_hasOp) {
01446     // Update MXj.
01447     MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
01448   }
01449   
01450   // Enter value on diagonal of B.
01451   (*B)(j,j) = diag;
01452       }
01453       else {
01454   // Create a random vector and orthogonalize it against all previous columns of Q.
01455   Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
01456   Teuchos::RCP<MV> tempMXj;
01457   MVT::MvRandom( *tempXj );
01458   if (this->_hasOp) {
01459     tempMXj = MVT::Clone( X, 1 );
01460     OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
01461   } 
01462   else {
01463     tempMXj = tempXj;
01464   }
01465   MVT::MvDot( *tempXj, *tempMXj, oldDot );
01466   //
01467   for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
01468     
01469     for (int i=0; i<Q.size(); i++) {
01470       Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
01471       
01472       // Apply another step of classical Gram-Schmidt
01473       MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*tempXj,tempMXj,product);
01474       MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
01475       
01476       // Update MXj, with the least number of applications of Op as possible
01477       if (this->_hasOp) {
01478         if (MQ[i].get()) {
01479     // MQ was allocated and computed above; use it
01480     MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
01481         }
01482         else if (xc <= qcs[i]) {
01483     // MQ was not allocated and computed above; it was cheaper to use X before and it still is
01484     OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
01485         }
01486       }
01487     } // for (int i=0; i<nq; i++)
01488     
01489   }
01490   
01491   // Compute the Op-norms after the correction step.
01492   MVT::MvDot( *tempXj, *tempMXj, newDot );
01493   
01494   // Copy vector into current column of Xj
01495   if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
01496     ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
01497     
01498     // Enter value on diagonal of B.
01499     (*B)(j,j) = ZERO;
01500 
01501     // Copy vector into current column of _basisvecs
01502     MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
01503     if (this->_hasOp) {
01504       MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
01505     }
01506   }
01507   else {
01508     return j;
01509   } 
01510       } // if (!dep_flg)
01511 
01512       // Remove the vectors from array
01513       if (j > 0) {
01514   Q.resize( nq );
01515   C.resize( nq );
01516   qcs.resize( nq );
01517       }
01518 
01519     } // for (int j=0; j<xc; j++)
01520 
01521     return xc;
01522   }
01523 
01524 } // namespace Belos
01525 
01526 #endif // BELOS_DGKS_ORTHOMANAGER_HPP
01527 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines