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