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 //
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>
49 #include <BelosMultiVecTraits.hpp>
50 #include <BelosOrthoManager.hpp>
51 #include <BelosOutputManager.hpp>
53 #include <Teuchos_ParameterListAcceptorDefaultBase.hpp>
54 #include <Teuchos_StandardCatchMacros.hpp>
55 #include <Teuchos_TimeMonitor.hpp>
56 
57 namespace Belos {
58 
67  template<class Scalar, class MV>
69  public OrthoManager<Scalar, MV>,
71  {
72  public:
73  typedef Scalar scalar_type;
77 
78  private:
82 
84  std::string label_;
88  bool reorthogonalize_;
90  bool useMgs_;
92  mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
93 
94 #ifdef BELOS_TEUCHOS_TIME_MONITOR
95  Teuchos::RCP<Teuchos::Time> timerOrtho_;
98  Teuchos::RCP<Teuchos::Time> timerProject_;
100  Teuchos::RCP<Teuchos::Time> timerNormalize_;
101 
111  makeTimer (const std::string& prefix,
112  const std::string& timerName)
113  {
114  const std::string timerLabel =
115  prefix.empty() ? timerName : (prefix + ": " + timerName);
116  return Teuchos::TimeMonitor::getNewCounter (timerLabel);
117  }
118 #endif // BELOS_TEUCHOS_TIME_MONITOR
119 
120  public:
121 
130  {
132  using Teuchos::parameterList;
133  using Teuchos::RCP;
134 
135  const std::string defaultNormalizationMethod ("MGS");
136  const bool defaultReorthogonalization = false;
137 
138  if (defaultParams_.is_null()) {
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 "
143  "Gram-Schmidt).");
144  params->set ("Reorthogonalization", defaultReorthogonalization,
145  "Whether to perform one (unconditional) reorthogonalization "
146  "pass.");
147  defaultParams_ = params;
148  }
149  return defaultParams_;
150  }
151 
159  {
161  using Teuchos::parameterList;
162  using Teuchos::RCP;
163  using Teuchos::rcp;
164 
165  const std::string fastNormalizationMethod ("CGS");
166  const bool fastReorthogonalization = false;
167 
168  // Start with a clone of the default parameters.
169  RCP<ParameterList> fastParams = parameterList (*getValidParameters());
170  fastParams->set ("Normalization", fastNormalizationMethod);
171  fastParams->set ("Reorthogonalization", fastReorthogonalization);
172 
173  return fastParams;
174  }
175 
176  void
178  {
180  using Teuchos::parameterList;
181  using Teuchos::RCP;
183 
184  RCP<const ParameterList> defaultParams = getValidParameters();
185  RCP<ParameterList> params;
186  if (plist.is_null ()) {
187  params = parameterList (*defaultParams);
188  } else {
189  params = plist;
190  params->validateParametersAndSetDefaults (*defaultParams);
191  }
192  const std::string normalizeImpl = params->get<std::string>("Normalization");
193  const bool reorthogonalize = params->get<bool>("Reorthogonalization");
194 
195  if (normalizeImpl == "MGS" ||
196  normalizeImpl == "Mgs" ||
197  normalizeImpl == "mgs") {
198  useMgs_ = true;
199  params->set ("Normalization", std::string ("MGS")); // Standardize.
200  } else {
201  useMgs_ = false;
202  params->set ("Normalization", std::string ("CGS")); // Standardize.
203  }
204  reorthogonalize_ = reorthogonalize;
205 
206  setMyParamList (params);
207  }
208 
220  const std::string& label,
221  const Teuchos::RCP<Teuchos::ParameterList>& params) :
222  label_ (label),
223  outMan_ (outMan)
224  {
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
230 
231  setParameterList (params);
232  if (! outMan_.is_null ()) {
233  using std::endl;
234  std::ostream& dbg = outMan_->stream(Debug);
235  dbg << "Belos::SimpleOrthoManager constructor:" << endl
236  << "-- Normalization method: "
237  << (useMgs_ ? "MGS" : "CGS") << endl
238  << "-- Reorthogonalize (unconditionally)? "
239  << (reorthogonalize_ ? "Yes" : "No") << endl;
240  }
241  }
242 
246  SimpleOrthoManager (const std::string& label = "Belos") :
247  label_ (label)
248  {
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
254 
255  setParameterList (Teuchos::null);
256  }
257 
259  virtual ~SimpleOrthoManager() {}
260 
261  void innerProd (const MV &X, const MV &Y, mat_type& Z) const {
262  MVT::MvTransMv (STS::one (), X, Y, Z);
263  }
264 
265  void norm (const MV& X, std::vector<magnitude_type>& normVec) const {
266  const int numCols = MVT::GetNumberVecs (X);
267  // std::vector<T>::size_type is unsigned; int is signed. Mixed
268  // unsigned/signed comparisons trigger compiler warnings.
269  if (normVec.size () < static_cast<size_t> (numCols)) {
270  normVec.resize (numCols); // Resize normvec if necessary.
271  }
272  MVT::MvNorm (X, normVec);
273  }
274 
275  void
276  project (MV &X,
279  {
280 #ifdef BELOS_TEUCHOS_TIME_MONITOR
281  Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
282  Teuchos::TimeMonitor timerMonitorProject(*timerProject_);
283 #endif // BELOS_TEUCHOS_TIME_MONITOR
284 
285  allocateProjectionCoefficients (C, Q, X, true);
286  rawProject (X, Q, C);
287  if (reorthogonalize_) { // Unconditional reorthogonalization
289  allocateProjectionCoefficients (C2, Q, X, false);
290  for (int k = 0; k < Q.size(); ++k)
291  *C[k] += *C2[k];
292  }
293  }
294 
295  int
296  normalize (MV &X, mat_ptr B) const
297  {
298 #ifdef BELOS_TEUCHOS_TIME_MONITOR
299  Teuchos::TimeMonitor timerMonitorOrtho(*timerOrtho_);
300  Teuchos::TimeMonitor timerMonitorProject(*timerNormalize_);
301 #endif // BELOS_TEUCHOS_TIME_MONITOR
302 
303  if (useMgs_) {
304  return normalizeMgs (X, B);
305  } else {
306  return normalizeCgs (X, B);
307  }
308  }
309 
310  protected:
311  virtual int
314  mat_ptr B,
316  {
317  // Don't need time monitors here: project() and normalize() have
318  // their own.
319  this->project (X, C, Q);
320  return this->normalize (X, B);
321  }
322 
323  public:
324 
326  orthonormError(const MV &X) const
327  {
328  const Scalar ONE = STS::one();
329  const int ncols = MVT::GetNumberVecs(X);
330  mat_type XTX (ncols, ncols);
331  innerProd (X, X, XTX);
332  for (int k = 0; k < ncols; ++k) {
333  XTX(k,k) -= ONE;
334  }
335  return XTX.normFrobenius();
336  }
337 
339  orthogError(const MV &X1, const MV &X2) const
340  {
341  const int ncols_X1 = MVT::GetNumberVecs (X1);
342  const int ncols_X2 = MVT::GetNumberVecs (X2);
343  mat_type X1_T_X2 (ncols_X1, ncols_X2);
344  innerProd (X1, X2, X1_T_X2);
345  return X1_T_X2.normFrobenius();
346  }
347 
348  void setLabel (const std::string& label) { label_ = label; }
349  const std::string& getLabel() const { return label_; }
350 
351  private:
352 
353  int
354  normalizeMgs (MV &X,
356  {
357  using Teuchos::Range1D;
358  using Teuchos::RCP;
359  using Teuchos::rcp;
360  using Teuchos::View;
361 
362  const int numCols = MVT::GetNumberVecs (X);
363  if (numCols == 0) {
364  return 0;
365  }
366 
367  if (B.is_null ()) {
368  B = rcp (new mat_type (numCols, numCols));
369  } else if (B->numRows () != numCols || B->numCols () != numCols) {
370  B->shape (numCols, numCols);
371  }
372 
373  // Modified Gram-Schmidt orthogonalization
374  std::vector<magnitude_type> normVec (1);
375  for (int j = 0; j < numCols; ++j) {
376  RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(j, j));
377  MV& X_j = *X_cur;
378  for (int i = 0; i < j; ++i) {
379  RCP<const MV> X_prv = MVT::CloneView (X, Range1D(i, i));
380  const MV& X_i = *X_prv;
381  mat_type B_ij (View, *B, 1, 1, i, j);
382  innerProd (X_i, X_j, B_ij);
383  MVT::MvTimesMatAddMv (-STS::one(), X_i, B_ij, STS::one(), X_j);
384  if (reorthogonalize_) { // Unconditional reorthogonalization
385  // innerProd() overwrites B(i,j), so save the
386  // first-pass projection coefficient and update
387  // B(i,j) afterwards.
388  const Scalar B_ij_first = (*B)(i, j);
389  innerProd (X_i, X_j, B_ij);
390  MVT::MvTimesMatAddMv (-STS::one(), X_i, B_ij, STS::one(), X_j);
391  (*B)(i, j) += B_ij_first;
392  }
393  }
394  // Normalize column j of X
395  norm (X_j, normVec);
396  const magnitude_type theNorm = normVec[0];
397  (*B)(j, j) = theNorm;
398  if (normVec[0] != STM::zero()) {
399  MVT::MvScale (X_j, STS::one() / theNorm);
400  } else {
401  return j; // break out early
402  }
403  }
404  return numCols; // full rank, as far as we know
405  }
406 
407 
408  int
409  normalizeCgs (MV &X, mat_ptr B) const
410  {
411  using Teuchos::Range1D;
412  using Teuchos::RCP;
413  using Teuchos::rcp;
414  using Teuchos::View;
415 
416  const int numCols = MVT::GetNumberVecs (X);
417  if (numCols == 0) {
418  return 0;
419  }
420 
421  if (B.is_null ()) {
422  B = rcp (new mat_type (numCols, numCols));
423  } else if (B->numRows() != numCols || B->numCols() != numCols) {
424  B->shape (numCols, numCols);
425  }
426  mat_type& B_ref = *B;
427 
428  // Classical Gram-Schmidt orthogonalization
429  std::vector<magnitude_type> normVec (1);
430 
431  // Space for reorthogonalization
432  mat_type B2 (numCols, numCols);
433 
434  // Do the first column first.
435  {
436  RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(0, 0));
437  // Normalize column 0 of X
438  norm (*X_cur, normVec);
439  const magnitude_type theNorm = normVec[0];
440  B_ref(0,0) = theNorm;
441  if (theNorm != STM::zero ()) {
442  const Scalar invNorm = STS::one () / theNorm;
443  MVT::MvScale (*X_cur, invNorm);
444  }
445  else {
446  return 0; // break out early
447  }
448  }
449 
450  // Orthogonalize the remaining columns of X
451  for (int j = 1; j < numCols; ++j) {
452  RCP<MV> X_cur = MVT::CloneViewNonConst (X, Range1D(j, j));
453  RCP<const MV> X_prv = MVT::CloneView (X, Range1D(0, j-1));
454  mat_type B_prvcur (View, B_ref, j, 1, 0, j);
455 
456  // Project X_cur against X_prv (first pass)
457  innerProd (*X_prv, *X_cur, B_prvcur);
458  MVT::MvTimesMatAddMv (-STS::one(), *X_prv, B_prvcur, STS::one(), *X_cur);
459  // Unconditional reorthogonalization:
460  // project X_cur against X_prv (second pass)
461  if (reorthogonalize_) {
462  mat_type B2_prvcur (View, B2, j, 1, 0, j);
463  innerProd (*X_prv, *X_cur, B2_prvcur);
464  MVT::MvTimesMatAddMv (-STS::one(), *X_prv, B2_prvcur, STS::one(), *X_cur);
465  B_prvcur += B2_prvcur;
466  }
467  // Normalize column j of X
468  norm (*X_cur, normVec);
469  const magnitude_type theNorm = normVec[0];
470  B_ref(j,j) = theNorm;
471  if (theNorm != STM::zero ()) {
472  const Scalar invNorm = STS::one () / theNorm;
473  MVT::MvScale (*X_cur, invNorm);
474  }
475  else {
476  return j; // break out early
477  }
478  }
479  return numCols; // full rank, as far as we know
480  }
481 
482 
483  void
484  allocateProjectionCoefficients (Teuchos::Array<mat_ptr>& C,
486  const MV& X,
487  const bool attemptToRecycle = true) const
488  {
489  using Teuchos::rcp;
490 
491  const int num_Q_blocks = Q.size();
492  const int ncols_X = MVT::GetNumberVecs (X);
493  C.resize (num_Q_blocks);
494  // # of block(s) that had to be (re)allocated (either allocated
495  // freshly, or resized).
496  int numAllocated = 0;
497  if (attemptToRecycle) {
498  for (int i = 0; i < num_Q_blocks; ++i) {
499  const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
500  // Create a new C[i] if necessary, otherwise resize if
501  // necessary, otherwise fill with zeros.
502  if (C[i].is_null ()) {
503  C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
504  numAllocated++;
505  }
506  else {
507  mat_type& Ci = *C[i];
508  if (Ci.numRows() != ncols_Qi || Ci.numCols() != ncols_X) {
509  Ci.shape (ncols_Qi, ncols_X);
510  numAllocated++;
511  }
512  else {
513  Ci.putScalar (STS::zero());
514  }
515  }
516  }
517  }
518  else { // Just allocate; don't try to check if we can recycle
519  for (int i = 0; i < num_Q_blocks; ++i) {
520  const int ncols_Qi = MVT::GetNumberVecs (*Q[i]);
521  C[i] = rcp (new mat_type (ncols_Qi, ncols_X));
522  numAllocated++;
523  }
524  }
525  if (! outMan_.is_null()) {
526  using std::endl;
527  std::ostream& dbg = outMan_->stream(Debug);
528  dbg << "SimpleOrthoManager::allocateProjectionCoefficients: "
529  << "Allocated " << numAllocated << " blocks out of "
530  << num_Q_blocks << endl;
531  }
532  }
533 
534  void
535  rawProject (MV& X,
538  {
539  // "Modified Gram-Schmidt" version of Block Gram-Schmidt.
540  const int num_Q_blocks = Q.size();
541  for (int i = 0; i < num_Q_blocks; ++i) {
542  mat_type& Ci = *C[i];
543  const MV& Qi = *Q[i];
544  innerProd (Qi, X, Ci);
545  MVT::MvTimesMatAddMv (-STS::one(), Qi, Ci, STS::one(), X);
546  }
547  }
548 
549  };
550 } // namespace Belos
551 
552 #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 Fri Jun 5 2020 10:20:46 for Belos by doxygen 1.8.5