Belos Version of the Day
BelosSimpleOrthoManager.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 
00045 #ifndef __Belos_SimpleOrthoManager_hpp
00046 #define __Belos_SimpleOrthoManager_hpp
00047 
00048 #include <BelosConfigDefs.hpp>
00049 #include <BelosMultiVecTraits.hpp>
00050 #include <BelosOrthoManager.hpp>
00051 #include <BelosOutputManager.hpp>
00052 #include <Teuchos_ParameterList.hpp>
00053 #include <Teuchos_ParameterListAcceptorDefaultBase.hpp>
00054 #include <Teuchos_StandardCatchMacros.hpp>
00055 #include <Teuchos_TimeMonitor.hpp>
00056 
00057 namespace Belos {
00058 
00067   template<class Scalar, class MV>
00068   class SimpleOrthoManager : 
00069     public OrthoManager<Scalar, MV>, 
00070     public Teuchos::ParameterListAcceptorDefaultBase
00071   {
00072   public:
00073     typedef Scalar scalar_type;
00074     typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitude_type;
00075     typedef Teuchos::SerialDenseMatrix<int, Scalar> mat_type;
00076     typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> > mat_ptr;
00077     
00078   private:
00079     typedef MultiVecTraits<Scalar, MV> MVT;
00080     typedef Teuchos::ScalarTraits<Scalar> STS;
00081     typedef Teuchos::ScalarTraits<magnitude_type> STM;
00082 
00084     std::string label_;
00086     Teuchos::RCP<OutputManager<Scalar> > outMan_;
00088     bool reorthogonalize_;
00090     bool useMgs_;
00092     mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
00093 
00094 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00095 
00096     Teuchos::RCP<Teuchos::Time> timerOrtho_;
00098     Teuchos::RCP<Teuchos::Time> timerProject_;
00100     Teuchos::RCP<Teuchos::Time> timerNormalize_;
00101 
00110     static Teuchos::RCP<Teuchos::Time>
00111     makeTimer (const std::string& prefix, 
00112          const std::string& timerName)
00113     {
00114       const std::string timerLabel = 
00115   prefix.empty() ? timerName : (prefix + ": " + timerName);
00116       return Teuchos::TimeMonitor::getNewCounter (timerLabel);
00117     }
00118 #endif // BELOS_TEUCHOS_TIME_MONITOR
00119     
00120   public:
00121 
00128     Teuchos::RCP<const Teuchos::ParameterList> 
00129     getValidParameters () const
00130     {
00131       using Teuchos::ParameterList;
00132       using Teuchos::parameterList;
00133       using Teuchos::RCP;
00134 
00135       const std::string defaultNormalizationMethod ("MGS");
00136       const bool defaultReorthogonalization = false;
00137 
00138       if (defaultParams_.is_null()) {
00139   RCP<ParameterList> params = parameterList ("Simple");
00140   params->set ("Normalization", defaultNormalizationMethod,
00141          "Which normalization method to use. Valid values are \"MGS\""
00142          " (for Modified Gram-Schmidt) and \"CGS\" (for Classical "
00143          "Gram-Schmidt).");
00144   params->set ("Reorthogonalization", defaultReorthogonalization,
00145          "Whether to perform one (unconditional) reorthogonalization "
00146          "pass.");
00147   defaultParams_ = params;
00148       }
00149       return defaultParams_;
00150     }
00151 
00161     static TEUCHOS_DEPRECATED Teuchos::RCP<const Teuchos::ParameterList> 
00162     getDefaultParameters ()
00163     {
00164       using Teuchos::ParameterList;
00165       using Teuchos::parameterList;
00166       using Teuchos::RCP;
00167 
00168       const std::string defaultNormalizationMethod ("MGS");
00169       const bool defaultReorthogonalization = false;
00170 
00171       RCP<ParameterList> params = parameterList ("Simple");
00172       params->set ("Normalization", defaultNormalizationMethod,
00173        "Which normalization method to use. Valid values are \"MGS\""
00174        " (for Modified Gram-Schmidt) and \"CGS\" (for Classical "
00175        "Gram-Schmidt).");
00176       params->set ("Reorthogonalization", defaultReorthogonalization,
00177        "Whether to perform one (unconditional) reorthogonalization "
00178        "pass.");
00179       return params;
00180     }
00181 
00187     Teuchos::RCP<const Teuchos::ParameterList> 
00188     getFastParameters ()
00189     {
00190       using Teuchos::ParameterList;
00191       using Teuchos::parameterList;
00192       using Teuchos::RCP;
00193       using Teuchos::rcp;
00194 
00195       const std::string fastNormalizationMethod ("CGS");
00196       const bool fastReorthogonalization = false;
00197 
00198       // Start with a clone of the default parameters.
00199       RCP<ParameterList> fastParams = parameterList (*getValidParameters());
00200       fastParams->set ("Normalization", fastNormalizationMethod);
00201       fastParams->set ("Reorthogonalization", fastReorthogonalization);
00202 
00203       return fastParams;
00204     }
00205 
00206     void 
00207     setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
00208     {
00209       using Teuchos::ParameterList;
00210       using Teuchos::parameterList;
00211       using Teuchos::RCP;
00212       using Teuchos::Exceptions::InvalidParameter;
00213 
00214       RCP<const ParameterList> defaultParams = getValidParameters();
00215       RCP<ParameterList> params;
00216       if (plist.is_null ()) {
00217   params = parameterList (*defaultParams);
00218       } else {
00219   params = plist;
00220   params->validateParametersAndSetDefaults (*defaultParams);
00221       }
00222       const std::string normalizeImpl = params->get<std::string>("Normalization");
00223       const bool reorthogonalize = params->get<bool>("Reorthogonalization");
00224 
00225       if (normalizeImpl == "MGS" || 
00226     normalizeImpl == "Mgs" || 
00227     normalizeImpl == "mgs") {
00228   useMgs_ = true;
00229   params->set ("Normalization", std::string ("MGS")); // Standardize.
00230       } else {
00231   useMgs_ = false;
00232   params->set ("Normalization", std::string ("CGS")); // Standardize.
00233       }
00234       reorthogonalize_ = reorthogonalize;
00235 
00236       setMyParamList (params);
00237     }
00238 
00249     SimpleOrthoManager (const Teuchos::RCP<OutputManager<Scalar> >& outMan,
00250       const std::string& label,
00251       const Teuchos::RCP<Teuchos::ParameterList>& params) :
00252       label_ (label),
00253       outMan_ (outMan)
00254     {
00255 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00256       timerOrtho_ = makeTimer (label, "All orthogonalization");
00257       timerProject_ = makeTimer (label, "Projection");
00258       timerNormalize_ = makeTimer (label, "Normalization");
00259 #endif // BELOS_TEUCHOS_TIME_MONITOR
00260 
00261       setParameterList (params);
00262       if (! outMan_.is_null ()) {
00263   using std::endl;
00264   std::ostream& dbg = outMan_->stream(Debug);
00265   dbg << "Belos::SimpleOrthoManager constructor:" << endl
00266       << "-- Normalization method: " 
00267       << (useMgs_ ? "MGS" : "CGS") << endl
00268       << "-- Reorthogonalize (unconditionally)? " 
00269       << (reorthogonalize_ ? "Yes" : "No") << endl;
00270       }
00271     }
00272 
00276     SimpleOrthoManager (const std::string& label) :
00277       label_ (label)
00278     {
00279 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00280       timerOrtho_ = makeTimer (label, "All orthogonalization");
00281       timerProject_ = makeTimer (label, "Projection");
00282       timerNormalize_ = makeTimer (label, "Normalization");
00283 #endif // BELOS_TEUCHOS_TIME_MONITOR
00284 
00285       setParameterList (Teuchos::null);
00286     }
00287     
00289     virtual ~SimpleOrthoManager() {}
00290       
00291     void innerProd (const MV &X, const MV &Y, mat_type& Z) const {
00292       MVT::MvTransMv (STS::one (), X, Y, Z);
00293     }
00294 
00295     void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
00296       const int numCols = MVT::GetNumberVecs (X);
00297       // std::vector<T>::size_type is unsigned; int is signed.  Mixed
00298       // unsigned/signed comparisons trigger compiler warnings.
00299       if (normVec.size () < static_cast<size_t> (numCols)) {
00300   normVec.resize (numCols); // Resize normvec if necessary.
00301       }
00302       MVT::MvNorm (X, normVec);
00303     }
00304 
00305     void 
00306     project (MV &X, 
00307        Teuchos::Array<mat_ptr> C,
00308        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00309     {
00310 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00311       Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
00312       Teuchos::TimeMonitor timerMonitorProject(*timerProject_);
00313 #endif // BELOS_TEUCHOS_TIME_MONITOR
00314 
00315       allocateProjectionCoefficients (C, Q, X, true);
00316       rawProject (X, Q, C);
00317       if (reorthogonalize_) { // Unconditional reorthogonalization
00318   Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C2;
00319   allocateProjectionCoefficients (C2, Q, X, false);
00320   for (int k = 0; k < Q.size(); ++k)
00321     *C[k] += *C2[k];
00322       }
00323     }
00324 
00325     int 
00326     normalize (MV &X, mat_ptr B) const
00327     {
00328 #ifdef BELOS_TEUCHOS_TIME_MONITOR
00329       Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
00330       Teuchos::TimeMonitor timerMonitorProject(*timerNormalize_);
00331 #endif // BELOS_TEUCHOS_TIME_MONITOR
00332 
00333       if (useMgs_) {
00334   return normalizeMgs (X, B);
00335       } else {
00336   return normalizeCgs (X, B);
00337       }
00338     }
00339 
00340     int 
00341     projectAndNormalize (MV &X, 
00342        Teuchos::Array<mat_ptr> C,
00343        mat_ptr B,
00344        Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
00345     {
00346       // Don't need time monitors here: project() and normalize() have
00347       // their own.
00348       project (X, C, Q);
00349       return normalize (X, B);
00350     }
00351 
00352     magnitude_type 
00353     orthonormError(const MV &X) const 
00354     {
00355       const Scalar ONE = STS::one();
00356       const int ncols = MVT::GetNumberVecs(X);
00357       mat_type XTX (ncols, ncols);
00358       innerProd (X, X, XTX);
00359       for (int k = 0; k < ncols; ++k) {
00360   XTX(k,k) -= ONE;
00361       }
00362       return XTX.normFrobenius();
00363     }
00364 
00365     magnitude_type
00366     orthogError(const MV &X1, const MV &X2) const
00367     {
00368       const int ncols_X1 = MVT::GetNumberVecs (X1);
00369       const int ncols_X2 = MVT::GetNumberVecs (X2);
00370       mat_type X1_T_X2 (ncols_X1, ncols_X2);
00371       innerProd (X1, X2, X1_T_X2);
00372       return X1_T_X2.normFrobenius();
00373     }
00374 
00375     void setLabel (const std::string& label) { label_ = label; }
00376     const std::string& getLabel() const { return label_; }
00377 
00378   private:
00379 
00380     int
00381     normalizeMgs (MV &X, 
00382       Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > B) const
00383     {
00384       using Teuchos::Range1D;
00385       using Teuchos::RCP;
00386       using Teuchos::rcp;
00387       using Teuchos::View;
00388 
00389       const int numCols = MVT::GetNumberVecs (X);
00390       if (numCols == 0) {
00391   return 0;
00392       }
00393 
00394       if (B.is_null ()) {
00395   B = rcp (new mat_type (numCols, numCols));
00396       } else if (B->numRows () != numCols || B->numCols () != numCols) {
00397   B->shape (numCols, numCols);
00398       }
00399 
00400       // Modified Gram-Schmidt orthogonalization
00401       std::vector<magnitude_type> normVec (1);
00402       for (int j = 0; j < numCols; ++j) {
00403   RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(j, j));
00404   MV& X_j = *X_cur;
00405   for (int i = 0; i < j; ++i) {
00406     RCP<const MV> X_prv = MVT::CloneView (X, Range1D(i, i));
00407     const MV& X_i = *X_prv;
00408     mat_type B_ij (View, *B, 1, 1, i, j);
00409     innerProd (X_i, X_j, B_ij);
00410     MVT::MvTimesMatAddMv (-STS::one(), X_i, B_ij, STS::one(), X_j);
00411     if (reorthogonalize_) { // Unconditional reorthogonalization
00412       // innerProd() overwrites B(i,j), so save the
00413       // first-pass projection coefficient and update
00414       // B(i,j) afterwards.
00415       const Scalar B_ij_first = (*B)(i, j);
00416       innerProd (X_i, X_j, B_ij);
00417       MVT::MvTimesMatAddMv (-STS::one(), X_i, B_ij, STS::one(), X_j);
00418       (*B)(i, j) += B_ij_first;
00419     }
00420   }
00421   // Normalize column j of X
00422   norm (X_j, normVec);
00423   const magnitude_type theNorm = normVec[0];
00424   (*B)(j, j) = theNorm;
00425   if (normVec[0] != STM::zero()) {
00426     MVT::MvScale (X_j, STS::one() / theNorm);
00427   } else {
00428     return j; // break out early
00429   }
00430       }
00431       return numCols; // full rank, as far as we know
00432     }
00433 
00434 
00435     int 
00436     normalizeCgs (MV &X, mat_ptr B) const
00437     {
00438       using Teuchos::Range1D;
00439       using Teuchos::RCP;
00440       using Teuchos::rcp;
00441       using Teuchos::View;
00442 
00443       const int numCols = MVT::GetNumberVecs (X);
00444       if (numCols == 0) {
00445   return 0;
00446       }
00447 
00448       if (B.is_null ()) {
00449   B = rcp (new mat_type (numCols, numCols));
00450       } else if (B->numRows() != numCols || B->numCols() != numCols) {
00451   B->shape (numCols, numCols);
00452       }
00453       mat_type& B_ref = *B;
00454 
00455       // Classical Gram-Schmidt orthogonalization 
00456       std::vector<magnitude_type> normVec (1);
00457 
00458       // Space for reorthogonalization
00459       mat_type B2 (numCols, numCols);
00460 
00461       // Do the first column first.
00462       {
00463   RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(0, 0));
00464   // Normalize column 0 of X
00465   norm (*X_cur, normVec);
00466   const magnitude_type theNorm = normVec[0];
00467   B_ref(0,0) = theNorm;
00468   if (theNorm != STM::zero ()) {
00469     const Scalar invNorm = STS::one () / theNorm;
00470     MVT::MvScale (*X_cur, invNorm);
00471   }
00472   else {
00473     return 0; // break out early
00474   }
00475       }
00476 
00477       // Orthogonalize the remaining columns of X
00478       for (int j = 1; j < numCols; ++j) {
00479   RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(j, j));
00480   RCP<const MV> X_prv = MVT::CloneView (X, Range1D(0, j-1));        
00481   mat_type B_prvcur (View, B_ref, j, 1, 0, j);
00482 
00483   // Project X_cur against X_prv (first pass)
00484   innerProd (*X_prv, *X_cur, B_prvcur);
00485   MVT::MvTimesMatAddMv (-STS::one(), *X_prv, B_prvcur, STS::one(), *X_cur);
00486   // Unconditional reorthogonalization: 
00487   // project X_cur against X_prv (second pass)
00488   if (reorthogonalize_) {
00489     mat_type B2_prvcur (View, B2, j, 1, 0, j);
00490     innerProd (*X_prv, *X_cur, B2_prvcur);
00491     MVT::MvTimesMatAddMv (-STS::one(), *X_prv, B2_prvcur, STS::one(), *X_cur);
00492     B_prvcur += B2_prvcur;
00493   }
00494   // Normalize column j of X
00495   norm (*X_cur, normVec);
00496   const magnitude_type theNorm = normVec[0];
00497   B_ref(j,j) = theNorm;
00498   if (theNorm != STM::zero ()) {
00499     const Scalar invNorm = STS::one () / theNorm;
00500     MVT::MvScale (*X_cur, invNorm);
00501   }
00502   else {
00503     return j; // break out early
00504   }
00505       }
00506       return numCols; // full rank, as far as we know
00507     }
00508 
00509 
00510     void
00511     allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C, 
00512             Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
00513             const MV& X,
00514             const bool attemptToRecycle = true) const
00515     {
00516       using Teuchos::rcp;
00517 
00518       const int num_Q_blocks = Q.size();
00519       const int ncols_X = MVT::GetNumberVecs (X);
00520       C.resize (num_Q_blocks);
00521       // # of block(s) that had to be (re)allocated (either allocated
00522       // freshly, or resized).
00523       int numAllocated = 0;
00524       if (attemptToRecycle) {
00525   for (int i = 0; i < num_Q_blocks; ++i) {
00526     const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
00527     // Create a new C[i] if necessary, otherwise resize if
00528     // necessary, otherwise fill with zeros.
00529     if (C[i].is_null ()) {
00530       C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
00531       numAllocated++;
00532     }
00533     else {
00534       mat_type& Ci = *C[i];
00535       if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X) {
00536         Ci.shape (ncols_Qi, ncols_X);
00537         numAllocated++;
00538       }
00539       else {
00540         Ci.putScalar (STS::zero());
00541       }
00542     }
00543   }
00544       }
00545       else { // Just allocate; don't try to check if we can recycle
00546   for (int i = 0; i < num_Q_blocks; ++i) {
00547     const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
00548     C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
00549     numAllocated++;
00550   }
00551       }
00552       if (! outMan_.is_null()) {
00553   using std::endl;
00554   std::ostream& dbg = outMan_->stream(Debug);
00555   dbg << "SimpleOrthoManager::allocateProjectionCoefficients: " 
00556       << "Allocated " << numAllocated << " blocks out of " 
00557       << num_Q_blocks << endl;
00558       }
00559     }
00560 
00561     void
00562     rawProject (MV& X, 
00563     Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
00564     Teuchos::ArrayView<mat_ptr> C) const
00565     { 
00566       // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
00567       const int num_Q_blocks = Q.size();
00568       for (int i = 0; i < num_Q_blocks; ++i) {
00569   mat_type& Ci = *C[i];
00570   const MV& Qi = *Q[i];
00571   innerProd (Qi, X, Ci);
00572   MVT::MvTimesMatAddMv (-STS::one(), Qi, Ci, STS::one(), X);
00573       }
00574     }
00575 
00576   };
00577 } // namespace Belos
00578 
00579 #endif // __Belos_SimpleOrthoManager_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines