Belos Version of the Day
Loading...
Searching...
No Matches
BelosSimpleOrthoManager.hpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Belos: Block Linear Solvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
45#ifndef __Belos_SimpleOrthoManager_hpp
46#define __Belos_SimpleOrthoManager_hpp
47
48#include <BelosConfigDefs.hpp>
50#include <BelosOrthoManager.hpp>
52#include <Teuchos_ParameterList.hpp>
53#include <Teuchos_StandardCatchMacros.hpp>
54#include <Teuchos_TimeMonitor.hpp>
55
56namespace Belos {
57
66 template<class Scalar, class MV>
68 public OrthoManager<Scalar, MV>
69 {
70 public:
71 typedef Scalar scalar_type;
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;
75
76 private:
78 typedef Teuchos::ScalarTraits<Scalar> STS;
79 typedef Teuchos::ScalarTraits<magnitude_type> STM;
80
82 std::string label_;
84 Teuchos::RCP<OutputManager<Scalar> > outMan_;
86 bool reorthogonalize_;
88 bool useMgs_;
90 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
91
92#ifdef BELOS_TEUCHOS_TIME_MONITOR
94 Teuchos::RCP<Teuchos::Time> timerOrtho_;
96 Teuchos::RCP<Teuchos::Time> timerProject_;
98 Teuchos::RCP<Teuchos::Time> timerNormalize_;
99
108 static Teuchos::RCP<Teuchos::Time>
109 makeTimer (const std::string& prefix,
110 const std::string& timerName)
111 {
112 const std::string timerLabel =
113 prefix.empty() ? timerName : (prefix + ": " + timerName);
114 return Teuchos::TimeMonitor::getNewCounter (timerLabel);
115 }
116#endif // BELOS_TEUCHOS_TIME_MONITOR
117
118 public:
119
126 Teuchos::RCP<const Teuchos::ParameterList>
128 {
129 using Teuchos::ParameterList;
130 using Teuchos::parameterList;
131 using Teuchos::RCP;
132
133 const std::string defaultNormalizationMethod ("MGS");
134 const bool defaultReorthogonalization = false;
135
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 "
141 "Gram-Schmidt).");
142 params->set ("Reorthogonalization", defaultReorthogonalization,
143 "Whether to perform one (unconditional) reorthogonalization "
144 "pass.");
145 defaultParams_ = params;
146 }
147 return defaultParams_;
148 }
149
155 Teuchos::RCP<const Teuchos::ParameterList>
157 {
158 using Teuchos::ParameterList;
159 using Teuchos::parameterList;
160 using Teuchos::RCP;
161 using Teuchos::rcp;
162
163 const std::string fastNormalizationMethod ("CGS");
164 const bool fastReorthogonalization = false;
165
166 // Start with a clone of the default parameters.
167 RCP<ParameterList> fastParams = parameterList (*getValidParameters());
168 fastParams->set ("Normalization", fastNormalizationMethod);
169 fastParams->set ("Reorthogonalization", fastReorthogonalization);
170
171 return fastParams;
172 }
173
174 void
175 setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
176 {
177 using Teuchos::ParameterList;
178 using Teuchos::parameterList;
179 using Teuchos::RCP;
180 using Teuchos::Exceptions::InvalidParameter;
181
182 RCP<const ParameterList> defaultParams = getValidParameters();
183 RCP<ParameterList> params;
184 if (plist.is_null ()) {
185 params = parameterList (*defaultParams);
186 } else {
187 params = plist;
188 params->validateParametersAndSetDefaults (*defaultParams);
189 }
190 const std::string normalizeImpl = params->get<std::string>("Normalization");
191 const bool reorthogonalize = params->get<bool>("Reorthogonalization");
192
193 if (normalizeImpl == "MGS" ||
194 normalizeImpl == "Mgs" ||
195 normalizeImpl == "mgs") {
196 useMgs_ = true;
197 params->set ("Normalization", std::string ("MGS")); // Standardize.
198 } else {
199 useMgs_ = false;
200 params->set ("Normalization", std::string ("CGS")); // Standardize.
201 }
202 reorthogonalize_ = reorthogonalize;
203
204 this->setMyParamList (params);
205 }
206
217 SimpleOrthoManager (const Teuchos::RCP<OutputManager<Scalar> >& outMan,
218 const std::string& label,
219 const Teuchos::RCP<Teuchos::ParameterList>& params) :
220 label_ (label),
221 outMan_ (outMan)
222 {
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
228
229 setParameterList (params);
230 if (! outMan_.is_null ()) {
231 using std::endl;
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;
238 }
239 }
240
244 SimpleOrthoManager (const std::string& label = "Belos") :
245 label_ (label)
246 {
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
252
253 setParameterList (Teuchos::null);
254 }
255
258
259 void innerProd (const MV &X, const MV &Y, mat_type& Z) const {
260 MVT::MvTransMv (STS::one (), X, Y, Z);
261 }
262
263 void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
264 const int numCols = MVT::GetNumberVecs (X);
265 // std::vector<T>::size_type is unsigned; int is signed. Mixed
266 // unsigned/signed comparisons trigger compiler warnings.
267 if (normVec.size () < static_cast<size_t> (numCols)) {
268 normVec.resize (numCols); // Resize normvec if necessary.
269 }
270 MVT::MvNorm (X, normVec);
271 }
272
273 void
274 project (MV &X,
275 Teuchos::Array<mat_ptr> C,
276 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
277 {
278#ifdef BELOS_TEUCHOS_TIME_MONITOR
279 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
280 Teuchos::TimeMonitor timerMonitorProject(*timerProject_);
281#endif // BELOS_TEUCHOS_TIME_MONITOR
282
283 allocateProjectionCoefficients (C, Q, X, true);
284 rawProject (X, Q, C);
285 if (reorthogonalize_) { // Unconditional reorthogonalization
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)
289 *C[k] += *C2[k];
290 }
291 }
292
293 int
294 normalize (MV &X, mat_ptr B) const
295 {
296#ifdef BELOS_TEUCHOS_TIME_MONITOR
297 Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
298 Teuchos::TimeMonitor timerMonitorProject(*timerNormalize_);
299#endif // BELOS_TEUCHOS_TIME_MONITOR
300
301 if (useMgs_) {
302 return normalizeMgs (X, B);
303 } else {
304 return normalizeCgs (X, B);
305 }
306 }
307
308 protected:
309 virtual int
311 Teuchos::Array<mat_ptr> C,
312 mat_ptr B,
313 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
314 {
315 // Don't need time monitors here: project() and normalize() have
316 // their own.
317 this->project (X, C, Q);
318 return this->normalize (X, B);
319 }
320
321 public:
322
324 orthonormError(const MV &X) const
325 {
326 const Scalar ONE = STS::one();
327 const int ncols = MVT::GetNumberVecs(X);
328 mat_type XTX (ncols, ncols);
329 innerProd (X, X, XTX);
330 for (int k = 0; k < ncols; ++k) {
331 XTX(k,k) -= ONE;
332 }
333 return XTX.normFrobenius();
334 }
335
337 orthogError(const MV &X1, const MV &X2) const
338 {
339 const int ncols_X1 = MVT::GetNumberVecs (X1);
340 const int ncols_X2 = MVT::GetNumberVecs (X2);
341 mat_type X1_T_X2 (ncols_X1, ncols_X2);
342 innerProd (X1, X2, X1_T_X2);
343 return X1_T_X2.normFrobenius();
344 }
345
346 void setLabel (const std::string& label) { label_ = label; }
347 const std::string& getLabel() const { return label_; }
348
349 private:
350
351 int
352 normalizeMgs (MV &X,
353 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,Scalar> > B) const
354 {
355 using Teuchos::Range1D;
356 using Teuchos::RCP;
357 using Teuchos::rcp;
358 using Teuchos::View;
359
360 const int numCols = MVT::GetNumberVecs (X);
361 if (numCols == 0) {
362 return 0;
363 }
364
365 if (B.is_null ()) {
366 B = rcp (new mat_type (numCols, numCols));
367 } else if (B->numRows () != numCols || B->numCols () != numCols) {
368 B->shape (numCols, numCols);
369 }
370
371 // Modified Gram-Schmidt orthogonalization
372 std::vector<magnitude_type> normVec (1);
373 for (int j = 0; j < numCols; ++j) {
374 RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(j, j));
375 MV& X_j = *X_cur;
376 for (int i = 0; i < j; ++i) {
377 RCP<const MV> X_prv = MVT::CloneView (X, Range1D(i, i));
378 const MV& X_i = *X_prv;
379 mat_type B_ij (View, *B, 1, 1, i, j);
380 innerProd (X_i, X_j, B_ij);
381 MVT::MvTimesMatAddMv (-STS::one(), X_i, B_ij, STS::one(), X_j);
382 if (reorthogonalize_) { // Unconditional reorthogonalization
383 // innerProd() overwrites B(i,j), so save the
384 // first-pass projection coefficient and update
385 // B(i,j) afterwards.
386 const Scalar B_ij_first = (*B)(i, j);
387 innerProd (X_i, X_j, B_ij);
388 MVT::MvTimesMatAddMv (-STS::one(), X_i, B_ij, STS::one(), X_j);
389 (*B)(i, j) += B_ij_first;
390 }
391 }
392 // Normalize column j of X
393 norm (X_j, normVec);
394 const magnitude_type theNorm = normVec[0];
395 (*B)(j, j) = theNorm;
396 if (normVec[0] != STM::zero()) {
397 MVT::MvScale (X_j, STS::one() / theNorm);
398 } else {
399 return j; // break out early
400 }
401 }
402 return numCols; // full rank, as far as we know
403 }
404
405
406 int
407 normalizeCgs (MV &X, mat_ptr B) const
408 {
409 using Teuchos::Range1D;
410 using Teuchos::RCP;
411 using Teuchos::rcp;
412 using Teuchos::View;
413
414 const int numCols = MVT::GetNumberVecs (X);
415 if (numCols == 0) {
416 return 0;
417 }
418
419 if (B.is_null ()) {
420 B = rcp (new mat_type (numCols, numCols));
421 } else if (B->numRows() != numCols || B->numCols() != numCols) {
422 B->shape (numCols, numCols);
423 }
424 mat_type& B_ref = *B;
425
426 // Classical Gram-Schmidt orthogonalization
427 std::vector<magnitude_type> normVec (1);
428
429 // Space for reorthogonalization
430 mat_type B2 (numCols, numCols);
431
432 // Do the first column first.
433 {
434 RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(0, 0));
435 // Normalize column 0 of X
436 norm (*X_cur, normVec);
437 const magnitude_type theNorm = normVec[0];
438 B_ref(0,0) = theNorm;
439 if (theNorm != STM::zero ()) {
440 const Scalar invNorm = STS::one () / theNorm;
441 MVT::MvScale (*X_cur, invNorm);
442 }
443 else {
444 return 0; // break out early
445 }
446 }
447
448 // Orthogonalize the remaining columns of X
449 for (int j = 1; j < numCols; ++j) {
450 RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(j, j));
451 RCP<const MV> X_prv = MVT::CloneView (X, Range1D(0, j-1));
452 mat_type B_prvcur (View, B_ref, j, 1, 0, j);
453
454 // Project X_cur against X_prv (first pass)
455 innerProd (*X_prv, *X_cur, B_prvcur);
456 MVT::MvTimesMatAddMv (-STS::one(), *X_prv, B_prvcur, STS::one(), *X_cur);
457 // Unconditional reorthogonalization:
458 // project X_cur against X_prv (second pass)
459 if (reorthogonalize_) {
460 mat_type B2_prvcur (View, B2, j, 1, 0, j);
461 innerProd (*X_prv, *X_cur, B2_prvcur);
462 MVT::MvTimesMatAddMv (-STS::one(), *X_prv, B2_prvcur, STS::one(), *X_cur);
463 B_prvcur += B2_prvcur;
464 }
465 // Normalize column j of X
466 norm (*X_cur, normVec);
467 const magnitude_type theNorm = normVec[0];
468 B_ref(j,j) = theNorm;
469 if (theNorm != STM::zero ()) {
470 const Scalar invNorm = STS::one () / theNorm;
471 MVT::MvScale (*X_cur, invNorm);
472 }
473 else {
474 return j; // break out early
475 }
476 }
477 return numCols; // full rank, as far as we know
478 }
479
480
481 void
482 allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C,
483 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
484 const MV& X,
485 const bool attemptToRecycle = true) const
486 {
487 using Teuchos::rcp;
488
489 const int num_Q_blocks = Q.size();
490 const int ncols_X = MVT::GetNumberVecs (X);
491 C.resize (num_Q_blocks);
492 // # of block(s) that had to be (re)allocated (either allocated
493 // freshly, or resized).
494 int numAllocated = 0;
495 if (attemptToRecycle) {
496 for (int i = 0; i < num_Q_blocks; ++i) {
497 const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
498 // Create a new C[i] if necessary, otherwise resize if
499 // necessary, otherwise fill with zeros.
500 if (C[i].is_null ()) {
501 C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
502 numAllocated++;
503 }
504 else {
505 mat_type& Ci = *C[i];
506 if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X) {
507 Ci.shape (ncols_Qi, ncols_X);
508 numAllocated++;
509 }
510 else {
511 Ci.putScalar (STS::zero());
512 }
513 }
514 }
515 }
516 else { // Just allocate; don't try to check if we can recycle
517 for (int i = 0; i < num_Q_blocks; ++i) {
518 const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
519 C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
520 numAllocated++;
521 }
522 }
523 if (! outMan_.is_null()) {
524 using std::endl;
525 std::ostream& dbg = outMan_->stream(Debug);
526 dbg << "SimpleOrthoManager::allocateProjectionCoefficients: "
527 << "Allocated " << numAllocated << " blocks out of "
528 << num_Q_blocks << endl;
529 }
530 }
531
532 void
533 rawProject (MV& X,
534 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q,
535 Teuchos::ArrayView<mat_ptr> C) const
536 {
537 // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
538 const int num_Q_blocks = Q.size();
539 for (int i = 0; i < num_Q_blocks; ++i) {
540 mat_type& Ci = *C[i];
541 const MV& Qi = *Q[i];
542 innerProd (Qi, X, Ci);
543 MVT::MvTimesMatAddMv (-STS::one(), Qi, Ci, STS::one(), X);
544 }
545 }
546
547 };
548} // namespace Belos
549
550#endif // __Belos_SimpleOrthoManager_hpp
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 > &params)
Constructor.
Teuchos::SerialDenseMatrix< int, Scalar > mat_type
virtual ~SimpleOrthoManager()
Virtual destructor for memory safety of derived classes.

Generated for Belos by doxygen 1.9.6