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_;
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;
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 "
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"));
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");
233 dbg <<
"Belos::SimpleOrthoManager constructor:" << endl
234 <<
"-- Normalization method: "
235 << (
useMgs_ ?
"MGS" :
"CGS") << endl
236 <<
"-- Reorthogonalize (unconditionally)? "
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_);
286 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C2;
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_);
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();
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);
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()) {
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);
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;
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));
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.
int normalizeMgs(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > B) const
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
MultiVecTraits< Scalar, MV > MVT
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
bool reorthogonalize_
Whether or not to do (unconditional) reorthogonalization.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > mat_ptr
bool useMgs_
Whether to use MGS or CGS in the normalize() step.
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 rawProject(MV &X, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q, Teuchos::ArrayView< mat_ptr > C) const
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.
Teuchos::ScalarTraits< Scalar > STS
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)
Teuchos::RCP< OutputManager< Scalar > > outMan_
Output manager (used mainly for debugging).
std::string label_
Label for Belos timer display.
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.
void allocateProjectionCoefficients(Teuchos::Array< mat_ptr > &C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q, const MV &X, const bool attemptToRecycle=true) const
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.
int normalizeCgs(MV &X, mat_ptr B) const
Teuchos::ScalarTraits< magnitude_type > STM