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];