Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosSimpleOrthoManager.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Belos: Block Linear Solvers Package
4 //
5 // Copyright 2004-2016 NTESS and the Belos contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
13 #ifndef __Belos_SimpleOrthoManager_hpp
14 #define __Belos_SimpleOrthoManager_hpp
15 
16 #include <BelosConfigDefs.hpp>
17 #include <BelosMultiVecTraits.hpp>
18 #include <BelosOrthoManager.hpp>
19 #include <BelosOutputManager.hpp>
21 #include <Teuchos_StandardCatchMacros.hpp>
22 #include <Teuchos_TimeMonitor.hpp>
23 
24 namespace Belos {
25 
34  template<class Scalar, class MV>
36  public OrthoManager<Scalar, MV>
37  {
38  public:
39  typedef Scalar scalar_type;
43 
44  private:
48 
50  std::string label_;
54  bool reorthogonalize_;
56  bool useMgs_;
58  mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
59 
60 #ifdef BELOS_TEUCHOS_TIME_MONITOR
61  Teuchos::RCP<Teuchos::Time> timerOrtho_;
64  Teuchos::RCP<Teuchos::Time> timerProject_;
66  Teuchos::RCP<Teuchos::Time> timerNormalize_;
67 
77  makeTimer (const std::string& prefix,
78  const std::string& timerName)
79  {
80  const std::string timerLabel =
81  prefix.empty() ? timerName : (prefix + ": " + timerName);
82  return Teuchos::TimeMonitor::getNewCounter (timerLabel);
83  }
84 #endif // BELOS_TEUCHOS_TIME_MONITOR
85 
86  public:
87 
96  {
98  using Teuchos::parameterList;
99  using Teuchos::RCP;
100 
101  const std::string defaultNormalizationMethod ("MGS");
102  const bool defaultReorthogonalization = false;
103 
104  if (defaultParams_.is_null()) {
105  RCP<ParameterList> params = parameterList ("Simple");
106  params->set ("Normalization", defaultNormalizationMethod,
107  "Which normalization method to use. Valid values are \"MGS\""
108  " (for Modified Gram-Schmidt) and \"CGS\" (for Classical "
109  "Gram-Schmidt).");
110  params->set ("Reorthogonalization", defaultReorthogonalization,
111  "Whether to perform one (unconditional) reorthogonalization "
112  "pass.");
113  defaultParams_ = params;
114  }
115  return defaultParams_;
116  }
117 
125  {
127  using Teuchos::parameterList;
128  using Teuchos::RCP;
129  using Teuchos::rcp;
130 
131  const std::string fastNormalizationMethod ("CGS");
132  const bool fastReorthogonalization = false;
133 
134  // Start with a clone of the default parameters.
135  RCP<ParameterList> fastParams = parameterList (*getValidParameters());
136  fastParams->set ("Normalization", fastNormalizationMethod);
137  fastParams->set ("Reorthogonalization", fastReorthogonalization);
138 
139  return fastParams;
140  }
141 
142  void
144  {
146  using Teuchos::parameterList;
147  using Teuchos::RCP;
149 
150  RCP<const ParameterList> defaultParams = getValidParameters();
151  RCP<ParameterList> params;
152  if (plist.is_null ()) {
153  params = parameterList (*defaultParams);
154  } else {
155  params = plist;
156  params->validateParametersAndSetDefaults (*defaultParams);
157  }
158  const std::string normalizeImpl = params->get<std::string>("Normalization");
159  const bool reorthogonalize = params->get<bool>("Reorthogonalization");
160 
161  if (normalizeImpl == "MGS" ||
162  normalizeImpl == "Mgs" ||
163  normalizeImpl == "mgs") {
164  useMgs_ = true;
165  params->set ("Normalization", std::string ("MGS")); // Standardize.
166  } else {
167  useMgs_ = false;
168  params->set ("Normalization", std::string ("CGS")); // Standardize.
169  }
170  reorthogonalize_ = reorthogonalize;
171 
172  this->setMyParamList (params);
173  }
174 
186  const std::string& label,
187  const Teuchos::RCP<Teuchos::ParameterList>& params) :
188  label_ (label),
189  outMan_ (outMan)
190  {
191 #ifdef BELOS_TEUCHOS_TIME_MONITOR
192  timerOrtho_ = makeTimer (label, "All orthogonalization");
193  timerProject_ = makeTimer (label, "Projection");
194  timerNormalize_ = makeTimer (label, "Normalization");
195 #endif // BELOS_TEUCHOS_TIME_MONITOR
196 
197  setParameterList (params);
198  if (! outMan_.is_null ()) {
199  using std::endl;
200  std::ostream& dbg = outMan_->stream(Debug);
201  dbg << "Belos::SimpleOrthoManager constructor:" << endl
202  << "-- Normalization method: "
203  << (useMgs_ ? "MGS" : "CGS") << endl
204  << "-- Reorthogonalize (unconditionally)? "
205  << (reorthogonalize_ ? "Yes" : "No") << endl;
206  }
207  }
208 
212  SimpleOrthoManager (const std::string& label = "Belos") :
213  label_ (label)
214  {
215 #ifdef BELOS_TEUCHOS_TIME_MONITOR
216  timerOrtho_ = makeTimer (label, "All orthogonalization");
217  timerProject_ = makeTimer (label, "Projection");
218  timerNormalize_ = makeTimer (label, "Normalization");
219 #endif // BELOS_TEUCHOS_TIME_MONITOR
220 
221  setParameterList (Teuchos::null);
222  }
223 
225  virtual ~SimpleOrthoManager() {}
226 
227  void innerProd (const MV &X, const MV &Y, mat_type& Z) const {
228  MVT::MvTransMv (STS::one (), X, Y, Z);
229  }
230 
231  void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
232  const int numCols = MVT::GetNumberVecs (X);
233  // std::vector<T>::size_type is unsigned; int is signed. Mixed
234  // unsigned/signed comparisons trigger compiler warnings.
235  if (normVec.size () < static_cast<size_t> (numCols)) {
236  normVec.resize (numCols); // Resize normvec if necessary.
237  }
238  MVT::MvNorm (X, normVec);
239  }
240 
241  void
242  project (MV &X,
245  {
246 #ifdef BELOS_TEUCHOS_TIME_MONITOR
247  Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
248  Teuchos::TimeMonitor timerMonitorProject(*timerProject_);
249 #endif // BELOS_TEUCHOS_TIME_MONITOR
250 
251  allocateProjectionCoefficients (C, Q, X, true);
252  rawProject (X, Q, C);
253  if (reorthogonalize_) { // Unconditional reorthogonalization
255  allocateProjectionCoefficients (C2, Q, X, false);
256  for (int k = 0; k < Q.size(); ++k)
257  *C[k] += *C2[k];
258  }
259  }
260 
261  int
262  normalize (MV &X, mat_ptr B) const
263  {
264 #ifdef BELOS_TEUCHOS_TIME_MONITOR
265  Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
266  Teuchos::TimeMonitor timerMonitorProject(*timerNormalize_);
267 #endif // BELOS_TEUCHOS_TIME_MONITOR
268 
269  if (useMgs_) {
270  return normalizeMgs (X, B);
271  } else {
272  return normalizeCgs (X, B);
273  }
274  }
275 
276  protected:
277  virtual int
280  mat_ptr B,
282  {
283  // Don't need time monitors here: project() and normalize() have
284  // their own.
285  this->project (X, C, Q);
286  return this->normalize (X, B);
287  }
288 
289  public:
290 
292  orthonormError(const MV &X) const
293  {
294  const Scalar ONE = STS::one();
295  const int ncols = MVT::GetNumberVecs(X);
296  mat_type XTX (ncols, ncols);
297  innerProd (X, X, XTX);
298  for (int k = 0; k < ncols; ++k) {
299  XTX(k,k) -= ONE;
300  }
301  return XTX.normFrobenius();
302  }
303 
305  orthogError(const MV &X1, const MV &X2) const
306  {
307  const int ncols_X1 = MVT::GetNumberVecs (X1);
308  const int ncols_X2 = MVT::GetNumberVecs (X2);
309  mat_type X1_T_X2 (ncols_X1, ncols_X2);
310  innerProd (X1, X2, X1_T_X2);
311  return X1_T_X2.normFrobenius();
312  }
313 
314  void setLabel (const std::string& label) { label_ = label; }
315  const std::string& getLabel() const { return label_; }
316 
317  private:
318 
319  int
320  normalizeMgs (MV &X,
322  {
323  using Teuchos::Range1D;
324  using Teuchos::RCP;
325  using Teuchos::rcp;
326  using Teuchos::View;
327 
328  const int numCols = MVT::GetNumberVecs (X);
329  if (numCols == 0) {
330  return 0;
331  }
332 
333  if (B.is_null ()) {
334  B = rcp (new mat_type (numCols, numCols));
335  } else if (B->numRows () != numCols || B->numCols () != numCols) {
336  B->shape (numCols, numCols);
337  }
338 
339  // Modified Gram-Schmidt orthogonalization
340  std::vector<magnitude_type> normVec (1);
341  for (int j = 0; j < numCols; ++j) {
342  RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(j, j));
343  MV& X_j = *X_cur;
344  for (int i = 0; i < j; ++i) {
345  RCP<const MV> X_prv = MVT::CloneView (X, Range1D(i, i));
346  const MV& X_i = *X_prv;
347  mat_type B_ij (View, *B, 1, 1, i, j);
348  innerProd (X_i, X_j, B_ij);
349  MVT::MvTimesMatAddMv (-STS::one(), X_i, B_ij, STS::one(), X_j);
350  if (reorthogonalize_) { // Unconditional reorthogonalization
351  // innerProd() overwrites B(i,j), so save the
352  // first-pass projection coefficient and update
353  // B(i,j) afterwards.
354  const Scalar B_ij_first = (*B)(i, j);
355  innerProd (X_i, X_j, B_ij);
356  MVT::MvTimesMatAddMv (-STS::one(), X_i, B_ij, STS::one(), X_j);
357  (*B)(i, j) += B_ij_first;
358  }
359  }
360  // Normalize column j of X
361  norm (X_j, normVec);
362  const magnitude_type theNorm = normVec[0];
363  (*B)(j, j) = theNorm;
364  if (normVec[0] != STM::zero()) {
365  MVT::MvScale (X_j, STS::one() / theNorm);
366  } else {
367  return j; // break out early
368  }
369  }
370  return numCols; // full rank, as far as we know
371  }
372 
373 
374  int
375  normalizeCgs (MV &X, mat_ptr B) const
376  {
377  using Teuchos::Range1D;
378  using Teuchos::RCP;
379  using Teuchos::rcp;
380  using Teuchos::View;
381 
382  const int numCols = MVT::GetNumberVecs (X);
383  if (numCols == 0) {
384  return 0;
385  }
386 
387  if (B.is_null ()) {
388  B = rcp (new mat_type (numCols, numCols));
389  } else if (B->numRows() != numCols || B->numCols() != numCols) {
390  B->shape (numCols, numCols);
391  }
392  mat_type& B_ref = *B;
393 
394  // Classical Gram-Schmidt orthogonalization
395  std::vector<magnitude_type> normVec (1);
396 
397  // Space for reorthogonalization
398  mat_type B2 (numCols, numCols);
399 
400  // Do the first column first.
401  {
402  RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(0, 0));
403  // Normalize column 0 of X
404  norm (*X_cur, normVec);
405  const magnitude_type theNorm = normVec[0];
406  B_ref(0,0) = theNorm;
407  if (theNorm != STM::zero ()) {
408  const Scalar invNorm = STS::one () / theNorm;
409  MVT::MvScale (*X_cur, invNorm);
410  }
411  else {
412  return 0; // break out early
413  }
414  }
415 
416  // Orthogonalize the remaining columns of X
417  for (int j = 1; j < numCols; ++j) {
418  RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(j, j));
419  RCP<const MV> X_prv = MVT::CloneView (X, Range1D(0, j-1));
420  mat_type B_prvcur (View, B_ref, j, 1, 0, j);
421 
422  // Project X_cur against X_prv (first pass)
423  innerProd (*X_prv, *X_cur, B_prvcur);
424  MVT::MvTimesMatAddMv (-STS::one(), *X_prv, B_prvcur, STS::one(), *X_cur);
425  // Unconditional reorthogonalization:
426  // project X_cur against X_prv (second pass)
427  if (reorthogonalize_) {
428  mat_type B2_prvcur (View, B2, j, 1, 0, j);
429  innerProd (*X_prv, *X_cur, B2_prvcur);
430  MVT::MvTimesMatAddMv (-STS::one(), *X_prv, B2_prvcur, STS::one(), *X_cur);
431  B_prvcur += B2_prvcur;
432  }
433  // Normalize column j of X
434  norm (*X_cur, normVec);
435  const magnitude_type theNorm = normVec[0];
436  B_ref(j,j) = theNorm;
437  if (theNorm != STM::zero ()) {
438  const Scalar invNorm = STS::one () / theNorm;
439  MVT::MvScale (*X_cur, invNorm);
440  }
441  else {
442  return j; // break out early
443  }
444  }
445  return numCols; // full rank, as far as we know
446  }
447 
448 
449  void
450  allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C,
452  const MV& X,
453  const bool attemptToRecycle = true) const
454  {
455  using Teuchos::rcp;
456 
457  const int num_Q_blocks = Q.size();
458  const int ncols_X = MVT::GetNumberVecs (X);
459  C.resize (num_Q_blocks);
460  // # of block(s) that had to be (re)allocated (either allocated
461  // freshly, or resized).
462  int numAllocated = 0;
463  if (attemptToRecycle) {
464  for (int i = 0; i < num_Q_blocks; ++i) {
465  const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
466  // Create a new C[i] if necessary, otherwise resize if
467  // necessary, otherwise fill with zeros.
468  if (C[i].is_null ()) {
469  C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
470  numAllocated++;
471  }
472  else {
473  mat_type& Ci = *C[i];
474  if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X) {
475  Ci.shape (ncols_Qi, ncols_X);
476  numAllocated++;
477  }
478  else {
479  Ci.putScalar (STS::zero());
480  }
481  }
482  }
483  }
484  else { // Just allocate; don't try to check if we can recycle
485  for (int i = 0; i < num_Q_blocks; ++i) {
486  const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
487  C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
488  numAllocated++;
489  }
490  }
491  if (! outMan_.is_null()) {
492  using std::endl;
493  std::ostream& dbg = outMan_->stream(Debug);
494  dbg << "SimpleOrthoManager::allocateProjectionCoefficients: "
495  << "Allocated " << numAllocated << " blocks out of "
496  << num_Q_blocks << endl;
497  }
498  }
499 
500  void
501  rawProject (MV& X,
504  {
505  // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
506  const int num_Q_blocks = Q.size();
507  for (int i = 0; i < num_Q_blocks; ++i) {
508  mat_type& Ci = *C[i];
509  const MV& Qi = *Q[i];
510  innerProd (Qi, X, Ci);
511  MVT::MvTimesMatAddMv (-STS::one(), Qi, Ci, STS::one(), X);
512  }
513  }
514 
515  };
516 } // namespace Belos
517 
518 #endif // __Belos_SimpleOrthoManager_hpp
void project(MV &X, Teuchos::Array< mat_ptr > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
void innerProd(const MV &X, const MV &Y, mat_type &Z) const
Provides the inner product defining the orthogonality concepts.
Belos&#39;s basic output manager for sending information of select verbosity levels to the appropriate ou...
magnitude_type orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
virtual ~SimpleOrthoManager()
Virtual destructor for memory safety of derived classes.
Class which manages the output and verbosity of the Belos solvers.
bool is_null(const boost::shared_ptr< T > &p)
int normalize(MV &X, mat_ptr B) const
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitude_type
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, Scalar > > mat_ptr
static RCP< Time > getNewCounter(const std::string &name)
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).
Declaration of basic traits for the multivector type.
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 > &params)
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 .
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Traits class which defines basic operations on multivectors.
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters()
Get a &quot;fast&quot; 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).
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Simple OrthoManager implementation for benchmarks.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a default list of parameters.
void norm(const MV &X, std::vector< magnitude_type > &normVec) const
void setMyParamList(const RCP< ParameterList > &paramList)
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
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 .
void resize(size_type new_size, const value_type &x=value_type())
SimpleOrthoManager(const std::string &label="Belos")
Constructor.
Templated virtual class for providing orthogonalization/orthonormalization methods.
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)
magnitude_type orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors.
Belos header file which uses auto-configuration information to include necessary C++ headers...
virtual int projectAndNormalizeImpl(MV &X, Teuchos::Array< mat_ptr > C, mat_ptr B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
bool is_null() const

Generated on Tue Jul 23 2024 09:24:22 for Belos by doxygen 1.8.5