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 93 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);
116 #endif // BELOS_TEUCHOS_TIME_MONITOR 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");
227 #endif // BELOS_TEUCHOS_TIME_MONITOR 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");
251 #endif // BELOS_TEUCHOS_TIME_MONITOR 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_);
281 #endif // BELOS_TEUCHOS_TIME_MONITOR 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_);
299 #endif // BELOS_TEUCHOS_TIME_MONITOR 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];
550 #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
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.
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.
Traits class which defines basic operations on multivectors.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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.
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 .
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.
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.
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...