45 #ifndef __Belos_SimpleOrthoManager_hpp 46 #define __Belos_SimpleOrthoManager_hpp 52 #include <Teuchos_ParameterList.hpp> 53 #include <Teuchos_ParameterListAcceptorDefaultBase.hpp> 54 #include <Teuchos_StandardCatchMacros.hpp> 55 #include <Teuchos_TimeMonitor.hpp> 67 template<
class Scalar,
class MV>
70 public Teuchos::ParameterListAcceptorDefaultBase
74 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType
magnitude_type;
75 typedef Teuchos::SerialDenseMatrix<int, Scalar>
mat_type;
76 typedef Teuchos::RCP<Teuchos::SerialDenseMatrix<int, Scalar> >
mat_ptr;
80 typedef Teuchos::ScalarTraits<Scalar>
STS;
81 typedef Teuchos::ScalarTraits<magnitude_type>
STM;
86 Teuchos::RCP<OutputManager<Scalar> >
outMan_;
94 #ifdef BELOS_TEUCHOS_TIME_MONITOR 95 Teuchos::RCP<Teuchos::Time> timerOrtho_;
98 Teuchos::RCP<Teuchos::Time> timerProject_;
100 Teuchos::RCP<Teuchos::Time> timerNormalize_;
110 static Teuchos::RCP<Teuchos::Time>
111 makeTimer (
const std::string& prefix,
112 const std::string& timerName)
114 const std::string timerLabel =
115 prefix.empty() ? timerName : (prefix +
": " + timerName);
116 return Teuchos::TimeMonitor::getNewCounter (timerLabel);
118 #endif // BELOS_TEUCHOS_TIME_MONITOR 128 Teuchos::RCP<const Teuchos::ParameterList>
131 using Teuchos::ParameterList;
132 using Teuchos::parameterList;
135 const std::string defaultNormalizationMethod (
"MGS");
136 const bool defaultReorthogonalization =
false;
139 RCP<ParameterList> params = parameterList (
"Simple");
140 params->set (
"Normalization", defaultNormalizationMethod,
141 "Which normalization method to use. Valid values are \"MGS\"" 142 " (for Modified Gram-Schmidt) and \"CGS\" (for Classical " 144 params->set (
"Reorthogonalization", defaultReorthogonalization,
145 "Whether to perform one (unconditional) reorthogonalization " 157 Teuchos::RCP<const Teuchos::ParameterList>
160 using Teuchos::ParameterList;
161 using Teuchos::parameterList;
165 const std::string fastNormalizationMethod (
"CGS");
166 const bool fastReorthogonalization =
false;
170 fastParams->set (
"Normalization", fastNormalizationMethod);
171 fastParams->set (
"Reorthogonalization", fastReorthogonalization);
179 using Teuchos::ParameterList;
180 using Teuchos::parameterList;
182 using Teuchos::Exceptions::InvalidParameter;
185 RCP<ParameterList> params;
186 if (plist.is_null ()) {
187 params = parameterList (*defaultParams);
190 params->validateParametersAndSetDefaults (*defaultParams);
192 const std::string normalizeImpl = params->get<std::string>(
"Normalization");
193 const bool reorthogonalize = params->get<
bool>(
"Reorthogonalization");
195 if (normalizeImpl ==
"MGS" ||
196 normalizeImpl ==
"Mgs" ||
197 normalizeImpl ==
"mgs") {
199 params->set (
"Normalization", std::string (
"MGS"));
202 params->set (
"Normalization", std::string (
"CGS"));
206 setMyParamList (params);
220 const std::string& label,
221 const Teuchos::RCP<Teuchos::ParameterList>& params) :
225 #ifdef BELOS_TEUCHOS_TIME_MONITOR 226 timerOrtho_ = makeTimer (label,
"All orthogonalization");
227 timerProject_ = makeTimer (label,
"Projection");
228 timerNormalize_ = makeTimer (label,
"Normalization");
229 #endif // BELOS_TEUCHOS_TIME_MONITOR 235 dbg <<
"Belos::SimpleOrthoManager constructor:" << endl
236 <<
"-- Normalization method: " 237 << (
useMgs_ ?
"MGS" :
"CGS") << endl
238 <<
"-- Reorthogonalize (unconditionally)? " 249 #ifdef BELOS_TEUCHOS_TIME_MONITOR 250 timerOrtho_ = makeTimer (label,
"All orthogonalization");
251 timerProject_ = makeTimer (label,
"Projection");
252 timerNormalize_ = makeTimer (label,
"Normalization");
253 #endif // BELOS_TEUCHOS_TIME_MONITOR 265 void norm (
const MV& X, std::vector<magnitude_type>& normVec)
const {
269 if (normVec.size () <
static_cast<size_t> (numCols)) {
270 normVec.resize (numCols);
277 Teuchos::Array<mat_ptr> C,
278 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 280 #ifdef BELOS_TEUCHOS_TIME_MONITOR 281 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
282 Teuchos::TimeMonitor timerMonitorProject(*timerProject_);
283 #endif // BELOS_TEUCHOS_TIME_MONITOR 288 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > > C2;
290 for (
int k = 0; k < Q.size(); ++k)
298 #ifdef BELOS_TEUCHOS_TIME_MONITOR 299 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
300 Teuchos::TimeMonitor timerMonitorProject(*timerNormalize_);
301 #endif // BELOS_TEUCHOS_TIME_MONITOR 313 Teuchos::Array<mat_ptr> C,
315 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 328 const Scalar ONE = STS::one();
332 for (
int k = 0; k < ncols; ++k) {
335 return XTX.normFrobenius();
343 mat_type X1_T_X2 (ncols_X1, ncols_X2);
345 return X1_T_X2.normFrobenius();
355 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > B)
const 357 using Teuchos::Range1D;
368 B = rcp (
new mat_type (numCols, numCols));
369 }
else if (B->numRows () != numCols || B->numCols () != numCols) {
370 B->shape (numCols, numCols);
374 std::vector<magnitude_type> normVec (1);
375 for (
int j = 0; j < numCols; ++j) {
378 for (
int i = 0; i < j; ++i) {
380 const MV& X_i = *X_prv;
381 mat_type B_ij (View, *B, 1, 1, i, j);
388 const Scalar B_ij_first = (*B)(i, j);
391 (*B)(i, j) += B_ij_first;
397 (*B)(j, j) = theNorm;
398 if (normVec[0] != STM::zero()) {
411 using Teuchos::Range1D;
422 B = rcp (
new mat_type (numCols, numCols));
423 }
else if (B->numRows() != numCols || B->numCols() != numCols) {
424 B->shape (numCols, numCols);
429 std::vector<magnitude_type> normVec (1);
438 norm (*X_cur, normVec);
440 B_ref(0,0) = theNorm;
441 if (theNorm != STM::zero ()) {
442 const Scalar invNorm = STS::one () / theNorm;
451 for (
int j = 1; j < numCols; ++j) {
454 mat_type B_prvcur (View, B_ref, j, 1, 0, j);
462 mat_type B2_prvcur (View, B2, j, 1, 0, j);
465 B_prvcur += B2_prvcur;
468 norm (*X_cur, normVec);
470 B_ref(j,j) = theNorm;
471 if (theNorm != STM::zero ()) {
472 const Scalar invNorm = STS::one () / theNorm;
485 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
487 const bool attemptToRecycle =
true)
const 491 const int num_Q_blocks = Q.size();
493 C.resize (num_Q_blocks);
496 int numAllocated = 0;
497 if (attemptToRecycle) {
498 for (
int i = 0; i < num_Q_blocks; ++i) {
502 if (C[i].is_null ()) {
503 C[i] = rcp (
new mat_type (ncols_Qi, ncols_X));
508 if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X) {
509 Ci.shape (ncols_Qi, ncols_X);
513 Ci.putScalar (STS::zero());
519 for (
int i = 0; i < num_Q_blocks; ++i) {
521 C[i] = rcp (
new mat_type (ncols_Qi, ncols_X));
528 dbg <<
"SimpleOrthoManager::allocateProjectionCoefficients: " 529 <<
"Allocated " << numAllocated <<
" blocks out of " 530 << num_Q_blocks << endl;
536 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
537 Teuchos::ArrayView<mat_ptr> C)
const 540 const int num_Q_blocks = Q.size();
541 for (
int i = 0; i < num_Q_blocks; ++i) {
543 const MV& Qi = *Q[i];
552 #endif // __Belos_SimpleOrthoManager_hpp Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
virtual int projectAndNormalizeImpl(MV &X, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a default list of parameters.
virtual ~SimpleOrthoManager()
Virtual destructor for memory safety of derived classes.
Class which manages the output and verbosity of the Belos solvers.
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > mat_ptr
int normalizeMgs(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > B) const
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).
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
Declaration of basic traits for the multivector type.
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.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
SimpleOrthoManager(const Teuchos::RCP< OutputManager< Scalar > > &outMan, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > ¶ms)
Constructor.
void rawProject(MV &X, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q, Teuchos::ArrayView< mat_ptr > C) const
int normalizeCgs(MV &X, mat_ptr B) const
Teuchos::ScalarTraits< Scalar > STS
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 int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Teuchos::RCP< OutputManager< Scalar > > outMan_
Output manager (used mainly for debugging).
Traits class which defines basic operations on multivectors.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
MultiVecTraits< Scalar, MV > MVT
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get a "fast" list of parameters.
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).
magnitude_type orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
Simple OrthoManager implementation for benchmarks.
Teuchos::ScalarTraits< magnitude_type > STM
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::SerialDenseMatrix< int, Scalar > mat_type
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
bool useMgs_
Whether to use MGS or CGS in the normalize() step.
magnitude_type orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors.
SimpleOrthoManager(const std::string &label="Belos")
Constructor.
void allocateProjectionCoefficients(Teuchos::Array< mat_ptr > &C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q, const MV &X, const bool attemptToRecycle=true) const
Templated virtual class for providing orthogonalization/orthonormalization methods.
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Provides the inner product defining the orthogonality concepts.
std::string label_
Label for Belos timer display.
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 void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
Belos header file which uses auto-configuration information to include necessary C++ headers...
bool reorthogonalize_
Whether or not to do (unconditional) reorthogonalization.