45#ifndef __Belos_SimpleOrthoManager_hpp
46#define __Belos_SimpleOrthoManager_hpp
52#include <Teuchos_ParameterList.hpp>
53#include <Teuchos_StandardCatchMacros.hpp>
54#include <Teuchos_TimeMonitor.hpp>
66 template<
class Scalar,
class MV>
72 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType
magnitude_type;
73 typedef Teuchos::SerialDenseMatrix<int, Scalar>
mat_type;
74 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> >
mat_ptr;
78 typedef Teuchos::ScalarTraits<Scalar> STS;
79 typedef Teuchos::ScalarTraits<magnitude_type> STM;
84 Teuchos::RCP<OutputManager<Scalar> > outMan_;
86 bool reorthogonalize_;
90 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
92#ifdef BELOS_TEUCHOS_TIME_MONITOR
94 Teuchos::RCP<Teuchos::Time> timerOrtho_;
96 Teuchos::RCP<Teuchos::Time> timerProject_;
98 Teuchos::RCP<Teuchos::Time> timerNormalize_;
108 static Teuchos::RCP<Teuchos::Time>
109 makeTimer (
const std::string& prefix,
110 const std::string& timerName)
112 const std::string timerLabel =
113 prefix.empty() ? timerName : (prefix +
": " + timerName);
114 return Teuchos::TimeMonitor::getNewCounter (timerLabel);
126 Teuchos::RCP<const Teuchos::ParameterList>
129 using Teuchos::ParameterList;
130 using Teuchos::parameterList;
133 const std::string defaultNormalizationMethod (
"MGS");
134 const bool defaultReorthogonalization =
false;
136 if (defaultParams_.is_null()) {
137 RCP<ParameterList> params = parameterList (
"Simple");
138 params->set (
"Normalization", defaultNormalizationMethod,
139 "Which normalization method to use. Valid values are \"MGS\""
140 " (for Modified Gram-Schmidt) and \"CGS\" (for Classical "
142 params->set (
"Reorthogonalization", defaultReorthogonalization,
143 "Whether to perform one (unconditional) reorthogonalization "
145 defaultParams_ = params;
147 return defaultParams_;
155 Teuchos::RCP<const Teuchos::ParameterList>
158 using Teuchos::ParameterList;
159 using Teuchos::parameterList;
163 const std::string fastNormalizationMethod (
"CGS");
164 const bool fastReorthogonalization =
false;
168 fastParams->set (
"Normalization", fastNormalizationMethod);
169 fastParams->set (
"Reorthogonalization", fastReorthogonalization);
177 using Teuchos::ParameterList;
178 using Teuchos::parameterList;
180 using Teuchos::Exceptions::InvalidParameter;
183 RCP<ParameterList> params;
184 if (plist.is_null ()) {
185 params = parameterList (*defaultParams);
188 params->validateParametersAndSetDefaults (*defaultParams);
190 const std::string normalizeImpl = params->get<std::string>(
"Normalization");
191 const bool reorthogonalize = params->get<
bool>(
"Reorthogonalization");
193 if (normalizeImpl ==
"MGS" ||
194 normalizeImpl ==
"Mgs" ||
195 normalizeImpl ==
"mgs") {
197 params->set (
"Normalization", std::string (
"MGS"));
200 params->set (
"Normalization", std::string (
"CGS"));
202 reorthogonalize_ = reorthogonalize;
204 this->setMyParamList (params);
218 const std::string& label,
219 const Teuchos::RCP<Teuchos::ParameterList>& params) :
223#ifdef BELOS_TEUCHOS_TIME_MONITOR
224 timerOrtho_ = makeTimer (label,
"All orthogonalization");
225 timerProject_ = makeTimer (label,
"Projection");
226 timerNormalize_ = makeTimer (label,
"Normalization");
230 if (! outMan_.is_null ()) {
232 std::ostream& dbg = outMan_->stream(
Debug);
233 dbg <<
"Belos::SimpleOrthoManager constructor:" << endl
234 <<
"-- Normalization method: "
235 << (useMgs_ ?
"MGS" :
"CGS") << endl
236 <<
"-- Reorthogonalize (unconditionally)? "
237 << (reorthogonalize_ ?
"Yes" :
"No") << endl;
247#ifdef BELOS_TEUCHOS_TIME_MONITOR
248 timerOrtho_ = makeTimer (label,
"All orthogonalization");
249 timerProject_ = makeTimer (label,
"Projection");
250 timerNormalize_ = makeTimer (label,
"Normalization");
263 void norm (
const MV& X, std::vector<magnitude_type>& normVec)
const {
267 if (normVec.size () <
static_cast<size_t> (numCols)) {
268 normVec.resize (numCols);
275 Teuchos::Array<mat_ptr> C,
276 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
278#ifdef BELOS_TEUCHOS_TIME_MONITOR
279 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
280 Teuchos::TimeMonitor timerMonitorProject(*timerProject_);
283 allocateProjectionCoefficients (C, Q, X,
true);
284 rawProject (X, Q, C);
285 if (reorthogonalize_) {
286 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C2;
287 allocateProjectionCoefficients (C2, Q, X,
false);
288 for (
int k = 0; k < Q.size(); ++k)
296#ifdef BELOS_TEUCHOS_TIME_MONITOR
297 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
298 Teuchos::TimeMonitor timerMonitorProject(*timerNormalize_);
302 return normalizeMgs (X, B);
304 return normalizeCgs (X, B);
311 Teuchos::Array<mat_ptr> C,
313 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
326 const Scalar ONE = STS::one();
330 for (
int k = 0; k < ncols; ++k) {
333 return XTX.normFrobenius();
341 mat_type X1_T_X2 (ncols_X1, ncols_X2);
343 return X1_T_X2.normFrobenius();
346 void setLabel (
const std::string& label) { label_ = label; }
347 const std::string&
getLabel()
const {
return label_; }
353 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > B)
const
355 using Teuchos::Range1D;
366 B = rcp (
new mat_type (numCols, numCols));
367 }
else if (B->numRows () != numCols || B->numCols () != numCols) {
368 B->shape (numCols, numCols);
372 std::vector<magnitude_type> normVec (1);
373 for (
int j = 0; j < numCols; ++j) {
376 for (
int i = 0; i < j; ++i) {
378 const MV& X_i = *X_prv;
379 mat_type B_ij (View, *B, 1, 1, i, j);
382 if (reorthogonalize_) {
386 const Scalar B_ij_first = (*B)(i, j);
389 (*B)(i, j) += B_ij_first;
395 (*B)(j, j) = theNorm;
396 if (normVec[0] != STM::zero()) {
407 normalizeCgs (MV &X,
mat_ptr B)
const
409 using Teuchos::Range1D;
420 B = rcp (
new mat_type (numCols, numCols));
421 }
else if (B->numRows() != numCols || B->numCols() != numCols) {
422 B->shape (numCols, numCols);
427 std::vector<magnitude_type> normVec (1);
436 norm (*X_cur, normVec);
438 B_ref(0,0) = theNorm;
439 if (theNorm != STM::zero ()) {
440 const Scalar invNorm = STS::one () / theNorm;
449 for (
int j = 1; j < numCols; ++j) {
452 mat_type B_prvcur (View, B_ref, j, 1, 0, j);
459 if (reorthogonalize_) {
460 mat_type B2_prvcur (View, B2, j, 1, 0, j);
463 B_prvcur += B2_prvcur;
466 norm (*X_cur, normVec);
468 B_ref(j,j) = theNorm;
469 if (theNorm != STM::zero ()) {
470 const Scalar invNorm = STS::one () / theNorm;
482 allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C,
483 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
485 const bool attemptToRecycle =
true)
const
489 const int num_Q_blocks = Q.size();
491 C.resize (num_Q_blocks);
494 int numAllocated = 0;
495 if (attemptToRecycle) {
496 for (
int i = 0; i < num_Q_blocks; ++i) {
500 if (C[i].is_null ()) {
501 C[i] = rcp (
new mat_type (ncols_Qi, ncols_X));
506 if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X) {
507 Ci.shape (ncols_Qi, ncols_X);
511 Ci.putScalar (STS::zero());
517 for (
int i = 0; i < num_Q_blocks; ++i) {
519 C[i] = rcp (
new mat_type (ncols_Qi, ncols_X));
523 if (! outMan_.is_null()) {
525 std::ostream& dbg = outMan_->stream(
Debug);
526 dbg <<
"SimpleOrthoManager::allocateProjectionCoefficients: "
527 <<
"Allocated " << numAllocated <<
" blocks out of "
528 << num_Q_blocks << endl;
534 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
535 Teuchos::ArrayView<mat_ptr> C)
const
538 const int num_Q_blocks = Q.size();
539 for (
int i = 0; i < num_Q_blocks; ++i) {
541 const MV& Qi = *Q[i];
Belos header file which uses auto-configuration information to include necessary C++ headers.
Declaration of basic traits for the multivector type.
Templated virtual class for providing orthogonalization/orthonormalization methods.
Class which manages the output and verbosity of the Belos solvers.
Traits class which defines basic operations on multivectors.
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, Teuchos::SerialDenseMatrix< int, ScalarType > &B)
Compute a dense matrix B through the matrix-matrix multiply .
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
static void MvNorm(const MV &mv, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, NormType type=TwoNorm)
Compute the norm of each individual vector of mv. Upon return, normvec[i] holds the value of ,...
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Simple OrthoManager implementation for benchmarks.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a default list of parameters.
magnitude_type orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors.
int normalize(MV &X, mat_ptr B) const
This method takes a multivector X and attempts to compute an orthonormal basis for ,...
void norm(const MV &X, std::vector< magnitude_type > &normVec) const
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > mat_ptr
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
virtual int projectAndNormalizeImpl(MV &X, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
SimpleOrthoManager(const std::string &label="Belos")
Constructor.
magnitude_type orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Project X against the (orthogonal) entries of Q.
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Provides the inner product defining the orthogonality concepts.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get a "fast" list of parameters.
SimpleOrthoManager(const Teuchos::RCP< OutputManager< Scalar > > &outMan, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Constructor.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
virtual ~SimpleOrthoManager()
Virtual destructor for memory safety of derived classes.