Belos  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
BelosIMGSOrthoManager.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 
10 
15 #ifndef BELOS_IMGS_ORTHOMANAGER_HPP
16 #define BELOS_IMGS_ORTHOMANAGER_HPP
17 
25 // #define ORTHO_DEBUG
26 
27 #include "BelosConfigDefs.hpp"
28 #include "BelosMultiVecTraits.hpp"
29 #include "BelosOperatorTraits.hpp"
30 #include "BelosMatOrthoManager.hpp"
31 #include "Teuchos_SerialDenseMatrix.hpp"
32 #include "Teuchos_SerialDenseVector.hpp"
33 
34 #include "Teuchos_as.hpp"
35 #ifdef BELOS_TEUCHOS_TIME_MONITOR
36 #include "Teuchos_TimeMonitor.hpp"
37 #endif // BELOS_TEUCHOS_TIME_MONITOR
38 
39 namespace Belos {
40 
42  template<class ScalarType, class MV, class OP>
43  Teuchos::RCP<Teuchos::ParameterList> getIMGSDefaultParameters ();
44 
46  template<class ScalarType, class MV, class OP>
47  Teuchos::RCP<Teuchos::ParameterList> getIMGSFastParameters();
48 
49  template<class ScalarType, class MV, class OP>
50  class IMGSOrthoManager :
51  public MatOrthoManager<ScalarType,MV,OP>
52  {
53  private:
54  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
55  typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
56  typedef Teuchos::ScalarTraits<ScalarType> SCT;
57  typedef MultiVecTraits<ScalarType,MV> MVT;
58  typedef OperatorTraits<ScalarType,MV,OP> OPT;
59 
60  public:
62 
63 
65  IMGSOrthoManager( const std::string& label = "Belos",
66  Teuchos::RCP<const OP> Op = Teuchos::null,
67  const int max_ortho_steps = max_ortho_steps_default_,
68  const MagnitudeType blk_tol = blk_tol_default_,
69  const MagnitudeType sing_tol = sing_tol_default_ )
70  : MatOrthoManager<ScalarType,MV,OP>(Op),
71  max_ortho_steps_( max_ortho_steps ),
72  blk_tol_( blk_tol ),
73  sing_tol_( sing_tol ),
74  label_( label )
75  {
76 #ifdef BELOS_TEUCHOS_TIME_MONITOR
77  std::stringstream ss;
78  ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
79 
80  std::string orthoLabel = ss.str() + ": Orthogonalization";
81  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
82 
83  std::string updateLabel = ss.str() + ": Ortho (Update)";
84  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
85 
86  std::string normLabel = ss.str() + ": Ortho (Norm)";
87  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
88 
89  std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
90  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
91 #endif
92  }
93 
95  IMGSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
96  const std::string& label = "Belos",
97  Teuchos::RCP<const OP> Op = Teuchos::null) :
98  MatOrthoManager<ScalarType,MV,OP>(Op),
99  max_ortho_steps_ (max_ortho_steps_default_),
100  blk_tol_ (blk_tol_default_),
101  sing_tol_ (sing_tol_default_),
102  label_ (label)
103  {
104  setParameterList (plist);
105 
106 #ifdef BELOS_TEUCHOS_TIME_MONITOR
107  std::stringstream ss;
108  ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
109 
110  std::string orthoLabel = ss.str() + ": Orthogonalization";
111  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
112 
113  std::string updateLabel = ss.str() + ": Ortho (Update)";
114  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
115 
116  std::string normLabel = ss.str() + ": Ortho (Norm)";
117  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
118 
119  std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
120  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
121 #endif
122  }
123 
125  ~IMGSOrthoManager() {}
127 
129 
130  void
131  setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
132  {
133  using Teuchos::Exceptions::InvalidParameterName;
134  using Teuchos::ParameterList;
135  using Teuchos::parameterList;
136  using Teuchos::RCP;
137 
138  RCP<const ParameterList> defaultParams = getValidParameters();
139  RCP<ParameterList> params;
140  if (plist.is_null()) {
141  params = parameterList (*defaultParams);
142  } else {
143  params = plist;
144  // Some users might want to specify "blkTol" as "depTol". Due
145  // to this case, we don't invoke
146  // validateParametersAndSetDefaults on params. Instead, we go
147  // through the parameter list one parameter at a time and look
148  // for alternatives.
149  }
150 
151  // Using temporary variables and fetching all values before
152  // setting the output arguments ensures the strong exception
153  // guarantee for this function: if an exception is thrown, no
154  // externally visible side effects (in this case, setting the
155  // output arguments) have taken place.
156  int maxNumOrthogPasses;
157  MagnitudeType blkTol;
158  MagnitudeType singTol;
159 
160  try {
161  maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
162  } catch (InvalidParameterName&) {
163  maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
164  params->set ("maxNumOrthogPasses", maxNumOrthogPasses);
165  }
166 
167  // Handling of the "blkTol" parameter is a special case. This
168  // is because some users may prefer to call this parameter
169  // "depTol" for consistency with DGKS. However, our default
170  // parameter list calls this "blkTol", and we don't want the
171  // default list's value to override the user's value. Thus, we
172  // first check the user's parameter list for both names, and
173  // only then access the default parameter list.
174  try {
175  blkTol = params->get<MagnitudeType> ("blkTol");
176  } catch (InvalidParameterName&) {
177  try {
178  blkTol = params->get<MagnitudeType> ("depTol");
179  // "depTol" is the wrong name, so remove it and replace with
180  // "blkTol". We'll set "blkTol" below.
181  params->remove ("depTol");
182  } catch (InvalidParameterName&) {
183  blkTol = defaultParams->get<MagnitudeType> ("blkTol");
184  }
185  params->set ("blkTol", blkTol);
186  }
187 
188  try {
189  singTol = params->get<MagnitudeType> ("singTol");
190  } catch (InvalidParameterName&) {
191  singTol = defaultParams->get<MagnitudeType> ("singTol");
192  params->set ("singTol", singTol);
193  }
194 
195  max_ortho_steps_ = maxNumOrthogPasses;
196  blk_tol_ = blkTol;
197  sing_tol_ = singTol;
198 
199  this->setMyParamList (params);
200  }
201 
202  Teuchos::RCP<const Teuchos::ParameterList>
203  getValidParameters () const
204  {
205  if (defaultParams_.is_null()) {
206  defaultParams_ = Belos::getIMGSDefaultParameters<ScalarType, MV, OP>();
207  }
208 
209  return defaultParams_;
210  }
211 
213 
215 
216 
218  void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
219 
221  void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
222 
224  MagnitudeType getBlkTol() const { return blk_tol_; }
225 
227  MagnitudeType getSingTol() const { return sing_tol_; }
228 
230 
231 
233 
234 
262  void project ( MV &X, Teuchos::RCP<MV> MX,
263  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
264  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
265 
266 
269  void project ( MV &X,
270  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
271  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
272  project(X,Teuchos::null,C,Q);
273  }
274 
275 
276 
301  int normalize ( MV &X, Teuchos::RCP<MV> MX,
302  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
303 
304 
307  int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
308  return normalize(X,Teuchos::null,B);
309  }
310 
311  protected:
353  virtual int
354  projectAndNormalizeWithMxImpl (MV &X,
355  Teuchos::RCP<MV> MX,
356  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
357  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
358  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
359 
360  public:
362 
364 
368  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
369  orthonormError(const MV &X) const {
370  return orthonormError(X,Teuchos::null);
371  }
372 
377  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
378  orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
379 
383  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
384  orthogError(const MV &X1, const MV &X2) const {
385  return orthogError(X1,Teuchos::null,X2);
386  }
387 
392  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
393  orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
394 
396 
398 
399 
402  void setLabel(const std::string& label);
403 
406  const std::string& getLabel() const { return label_; }
407 
409 
411 
412 
414  static const int max_ortho_steps_default_;
416  static const MagnitudeType blk_tol_default_;
418  static const MagnitudeType sing_tol_default_;
419 
421  static const int max_ortho_steps_fast_;
423  static const MagnitudeType blk_tol_fast_;
425  static const MagnitudeType sing_tol_fast_;
426 
428 
429  private:
430 
432  int max_ortho_steps_;
434  MagnitudeType blk_tol_;
436  MagnitudeType sing_tol_;
437 
439  std::string label_;
440 #ifdef BELOS_TEUCHOS_TIME_MONITOR
441  Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
442 #endif // BELOS_TEUCHOS_TIME_MONITOR
443 
445  mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
446 
448  int findBasis(MV &X, Teuchos::RCP<MV> MX,
449  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
450  bool completeBasis, int howMany = -1 ) const;
451 
453  bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
454  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
455  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
456 
458  bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
459  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
460  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
461 
475  int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
476  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
477  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
478  Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;
479  };
480 
481  // Set static variables.
482  template<class ScalarType, class MV, class OP>
483  const int IMGSOrthoManager<ScalarType,MV,OP>::max_ortho_steps_default_ = 1;
484 
485  template<class ScalarType, class MV, class OP>
486  const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
487  IMGSOrthoManager<ScalarType,MV,OP>::blk_tol_default_
488  = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
489  Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps() );
490 
491  template<class ScalarType, class MV, class OP>
492  const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
493  IMGSOrthoManager<ScalarType,MV,OP>::sing_tol_default_
494  = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
495 
496  template<class ScalarType, class MV, class OP>
497  const int IMGSOrthoManager<ScalarType,MV,OP>::max_ortho_steps_fast_ = 1;
498 
499  template<class ScalarType, class MV, class OP>
500  const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
501  IMGSOrthoManager<ScalarType,MV,OP>::blk_tol_fast_
502  = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
503 
504  template<class ScalarType, class MV, class OP>
505  const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
506  IMGSOrthoManager<ScalarType,MV,OP>::sing_tol_fast_
507  = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
508 
510  // Set the label for this orthogonalization manager and create new timers if it's changed
511  template<class ScalarType, class MV, class OP>
512  void IMGSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
513  {
514  if (label != label_) {
515  label_ = label;
516 #ifdef BELOS_TEUCHOS_TIME_MONITOR
517  std::stringstream ss;
518  ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
519 
520  std::string orthoLabel = ss.str() + ": Orthogonalization";
521  timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
522 
523  std::string updateLabel = ss.str() + ": Ortho (Update)";
524  timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
525 
526  std::string normLabel = ss.str() + ": Ortho (Norm)";
527  timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
528 
529  std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
530  timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
531 #endif
532  }
533  }
534 
536  // Compute the distance from orthonormality
537  template<class ScalarType, class MV, class OP>
538  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
539  IMGSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
540  const ScalarType ONE = SCT::one();
541  int rank = MVT::GetNumberVecs(X);
542  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
543  MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
544  for (int i=0; i<rank; i++) {
545  xTx(i,i) -= ONE;
546  }
547  return xTx.normFrobenius();
548  }
549 
551  // Compute the distance from orthogonality
552  template<class ScalarType, class MV, class OP>
553  typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
554  IMGSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
555  int r1 = MVT::GetNumberVecs(X1);
556  int r2 = MVT::GetNumberVecs(X2);
557  Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
558  MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
559  return xTx.normFrobenius();
560  }
561 
563  // Find an Op-orthonormal basis for span(X) - span(W)
564  template<class ScalarType, class MV, class OP>
565  int
566  IMGSOrthoManager<ScalarType, MV, OP>::
567  projectAndNormalizeWithMxImpl(MV &X,
568  Teuchos::RCP<MV> MX,
569  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
570  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
571  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
572  {
573  using Teuchos::Array;
574  using Teuchos::null;
575  using Teuchos::is_null;
576  using Teuchos::RCP;
577  using Teuchos::rcp;
578  using Teuchos::SerialDenseMatrix;
579  typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
580  typedef typename Array< RCP< const MV > >::size_type size_type;
581 
582 #ifdef BELOS_TEUCHOS_TIME_MONITOR
583  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
584 #endif
585 
586  ScalarType ONE = SCT::one();
587  const MagnitudeType ZERO = MGT::zero();
588 
589  int nq = Q.size();
590  int xc = MVT::GetNumberVecs( X );
591  ptrdiff_t xr = MVT::GetGlobalLength( X );
592  int rank = xc;
593 
594  // If the user doesn't want to store the normalization
595  // coefficients, allocate some local memory for them. This will
596  // go away at the end of this method.
597  if (is_null (B)) {
598  B = rcp (new serial_dense_matrix_type (xc, xc));
599  }
600  // Likewise, if the user doesn't want to store the projection
601  // coefficients, allocate some local memory for them. Also make
602  // sure that all the entries of C are the right size. We're going
603  // to overwrite them anyway, so we don't have to worry about the
604  // contents (other than to resize them if they are the wrong
605  // size).
606  if (C.size() < nq)
607  C.resize (nq);
608  for (size_type k = 0; k < nq; ++k)
609  {
610  const int numRows = MVT::GetNumberVecs (*Q[k]);
611  const int numCols = xc; // Number of vectors in X
612 
613  if (is_null (C[k]))
614  C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
615  else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
616  {
617  int err = C[k]->reshape (numRows, numCols);
618  TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
619  "IMGS orthogonalization: failed to reshape "
620  "C[" << k << "] (the array of block "
621  "coefficients resulting from projecting X "
622  "against Q[1:" << nq << "]).");
623  }
624  }
625 
626  /****** DO NOT MODIFY *MX IF _hasOp == false ******/
627  if (this->_hasOp) {
628  if (MX == Teuchos::null) {
629  // we need to allocate space for MX
630  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
631  OPT::Apply(*(this->_Op),X,*MX);
632  }
633  }
634  else {
635  // Op == I --> MX = X (ignore it if the user passed it in)
636  MX = Teuchos::rcp( &X, false );
637  }
638 
639  int mxc = MVT::GetNumberVecs( *MX );
640  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
641 
642  // short-circuit
643  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
644 
645  int numbas = 0;
646  for (int i=0; i<nq; i++) {
647  numbas += MVT::GetNumberVecs( *Q[i] );
648  }
649 
650  // check size of B
651  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
652  "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
653  // check size of X and MX
654  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
655  "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
656  // check size of X w.r.t. MX
657  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
658  "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
659  // check feasibility
660  //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
661  // "Belos::IMGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
662 
663  // Some flags for checking dependency returns from the internal orthogonalization methods
664  bool dep_flg = false;
665 
666  // Make a temporary copy of X and MX, just in case a block dependency is detected.
667  Teuchos::RCP<MV> tmpX, tmpMX;
668  tmpX = MVT::CloneCopy(X);
669  if (this->_hasOp) {
670  tmpMX = MVT::CloneCopy(*MX);
671  }
672 
673  if (xc == 1) {
674 
675  // Use the cheaper block orthogonalization.
676  // NOTE: Don't check for dependencies because the update has one vector.
677  dep_flg = blkOrtho1( X, MX, C, Q );
678 
679  // Normalize the new block X
680  if ( B == Teuchos::null ) {
681  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
682  }
683  std::vector<ScalarType> diag(xc);
684  {
685 #ifdef BELOS_TEUCHOS_TIME_MONITOR
686  Teuchos::TimeMonitor normTimer( *timerNorm_ );
687 #endif
688  MVT::MvDot( X, *MX, diag );
689  }
690  (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
691 
692  if (SCT::magnitude((*B)(0,0)) > ZERO) {
693  rank = 1;
694  MVT::MvScale( X, ONE/(*B)(0,0) );
695  if (this->_hasOp) {
696  // Update MXj.
697  MVT::MvScale( *MX, ONE/(*B)(0,0) );
698  }
699  }
700  }
701  else {
702 
703  // Use the cheaper block orthogonalization.
704  dep_flg = blkOrtho( X, MX, C, Q );
705 
706  // If a dependency has been detected in this block, then perform
707  // the more expensive nonblock (single vector at a time)
708  // orthogonalization.
709  if (dep_flg) {
710  rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
711 
712  // Copy tmpX back into X.
713  MVT::Assign( *tmpX, X );
714  if (this->_hasOp) {
715  MVT::Assign( *tmpMX, *MX );
716  }
717  }
718  else {
719  // There is no dependency, so orthonormalize new block X
720  rank = findBasis( X, MX, B, false );
721  if (rank < xc) {
722  // A dependency was found during orthonormalization of X,
723  // rerun orthogonalization using more expensive single-
724  // vector orthogonalization.
725  rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
726 
727  // Copy tmpX back into X.
728  MVT::Assign( *tmpX, X );
729  if (this->_hasOp) {
730  MVT::Assign( *tmpMX, *MX );
731  }
732  }
733  }
734  } // if (xc == 1) {
735 
736  // this should not raise an std::exception; but our post-conditions oblige us to check
737  TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
738  "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
739 
740  // Return the rank of X.
741  return rank;
742  }
743 
744 
745 
747  // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
748  template<class ScalarType, class MV, class OP>
749  int IMGSOrthoManager<ScalarType, MV, OP>::normalize(
750  MV &X, Teuchos::RCP<MV> MX,
751  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
752 
753 #ifdef BELOS_TEUCHOS_TIME_MONITOR
754  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
755 #endif
756 
757  // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
758  return findBasis(X, MX, B, true);
759  }
760 
761 
762 
764  template<class ScalarType, class MV, class OP>
765  void IMGSOrthoManager<ScalarType, MV, OP>::project(
766  MV &X, Teuchos::RCP<MV> MX,
767  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
768  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
769  // For the inner product defined by the operator Op or the identity (Op == 0)
770  // -> Orthogonalize X against each Q[i]
771  // Modify MX accordingly
772  //
773  // Note that when Op is 0, MX is not referenced
774  //
775  // Parameter variables
776  //
777  // X : Vectors to be transformed
778  //
779  // MX : Image of the block of vectors X by the mass matrix
780  //
781  // Q : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
782  //
783 
784 #ifdef BELOS_TEUCHOS_TIME_MONITOR
785  Teuchos::TimeMonitor orthotimer(*timerOrtho_);
786 #endif
787 
788  int xc = MVT::GetNumberVecs( X );
789  ptrdiff_t xr = MVT::GetGlobalLength( X );
790  int nq = Q.size();
791  std::vector<int> qcs(nq);
792  // short-circuit
793  if (nq == 0 || xc == 0 || xr == 0) {
794  return;
795  }
796  ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
797  // if we don't have enough C, expand it with null references
798  // if we have too many, resize to throw away the latter ones
799  // if we have exactly as many as we have Q, this call has no effect
800  C.resize(nq);
801 
802 
803  /****** DO NOT MODIFY *MX IF _hasOp == false ******/
804  if (this->_hasOp) {
805  if (MX == Teuchos::null) {
806  // we need to allocate space for MX
807  MX = MVT::Clone(X,MVT::GetNumberVecs(X));
808  OPT::Apply(*(this->_Op),X,*MX);
809  }
810  }
811  else {
812  // Op == I --> MX = X (ignore it if the user passed it in)
813  MX = Teuchos::rcp( &X, false );
814  }
815  int mxc = MVT::GetNumberVecs( *MX );
816  ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
817 
818  // check size of X and Q w.r.t. common sense
819  TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
820  "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
821  // check size of X w.r.t. MX and Q
822  TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
823  "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
824 
825  // tally up size of all Q and check/allocate C
826  for (int i=0; i<nq; i++) {
827  TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
828  "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
829  qcs[i] = MVT::GetNumberVecs( *Q[i] );
830  TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
831  "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
832 
833  // check size of C[i]
834  if ( C[i] == Teuchos::null ) {
835  C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
836  }
837  else {
838  TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
839  "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
840  }
841  }
842 
843  // Use the cheaper block orthogonalization, don't check for rank deficiency.
844  blkOrtho( X, MX, C, Q );
845 
846  }
847 
849  // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
850  // the rank is numvectors(X)
851  template<class ScalarType, class MV, class OP>
852  int IMGSOrthoManager<ScalarType, MV, OP>::findBasis(
853  MV &X, Teuchos::RCP<MV> MX,
854  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
855  bool completeBasis, int howMany ) const {
856  // For the inner product defined by the operator Op or the identity (Op == 0)
857  // -> Orthonormalize X
858  // Modify MX accordingly
859  //
860  // Note that when Op is 0, MX is not referenced
861  //
862  // Parameter variables
863  //
864  // X : Vectors to be orthonormalized
865  //
866  // MX : Image of the multivector X under the operator Op
867  //
868  // Op : Pointer to the operator for the inner product
869  //
870  //
871 
872  const ScalarType ONE = SCT::one();
873  const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
874 
875  int xc = MVT::GetNumberVecs( X );
876  ptrdiff_t xr = MVT::GetGlobalLength( X );
877 
878  if (howMany == -1) {
879  howMany = xc;
880  }
881 
882  /*******************************************************
883  * If _hasOp == false, we will not reference MX below *
884  *******************************************************/
885 
886  // if Op==null, MX == X (via pointer)
887  // Otherwise, either the user passed in MX or we will allocated and compute it
888  if (this->_hasOp) {
889  if (MX == Teuchos::null) {
890  // we need to allocate space for MX
891  MX = MVT::Clone(X,xc);
892  OPT::Apply(*(this->_Op),X,*MX);
893  }
894  }
895 
896  /* if the user doesn't want to store the coefficienets,
897  * allocate some local memory for them
898  */
899  if ( B == Teuchos::null ) {
900  B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
901  }
902 
903  int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
904  ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
905 
906  // check size of C, B
907  TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
908  "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
909  TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
910  "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
911  TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
912  "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
913  TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
914  "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
915  TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
916  "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
917 
918  /* xstart is which column we are starting the process with, based on howMany
919  * columns before xstart are assumed to be Op-orthonormal already
920  */
921  int xstart = xc - howMany;
922 
923  for (int j = xstart; j < xc; j++) {
924 
925  // numX is
926  // * number of currently orthonormal columns of X
927  // * the index of the current column of X
928  int numX = j;
929  bool addVec = false;
930 
931  // Get a view of the vector currently being worked on.
932  std::vector<int> index(1);
933  index[0] = numX;
934  Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
935  Teuchos::RCP<MV> MXj;
936  if ((this->_hasOp)) {
937  // MXj is a view of the current vector in MX
938  MXj = MVT::CloneViewNonConst( *MX, index );
939  }
940  else {
941  // MXj is a pointer to Xj, and MUST NOT be modified
942  MXj = Xj;
943  }
944 
945  Teuchos::RCP<MV> oldMXj;
946  if (numX > 0) {
947  oldMXj = MVT::CloneCopy( *MXj );
948  }
949 
950  // Make storage for these Gram-Schmidt iterations.
951  Teuchos::SerialDenseVector<int,ScalarType> product(numX);
952  Teuchos::SerialDenseVector<int,ScalarType> P2(1);
953  Teuchos::RCP<const MV> prevX, prevMX;
954 
955  std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
956  //
957  // Save old MXj vector and compute Op-norm
958  //
959  {
960 #ifdef BELOS_TEUCHOS_TIME_MONITOR
961  Teuchos::TimeMonitor normTimer( *timerNorm_ );
962 #endif
963  MVT::MvDot( *Xj, *MXj, oldDot );
964  }
965  // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
966  TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
967  "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
968 
969  // Perform MGS one vector at a time
970  for (int ii=0; ii<numX; ii++) {
971 
972  index[0] = ii;
973  prevX = MVT::CloneView( X, index );
974  if (this->_hasOp) {
975  prevMX = MVT::CloneView( *MX, index );
976  }
977 
978  for (int i=0; i<max_ortho_steps_; ++i) {
979 
980  // product <- prevX^T MXj
981  {
982 #ifdef BELOS_TEUCHOS_TIME_MONITOR
983  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
984 #endif
985  MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
986  }
987 
988  // Xj <- Xj - prevX prevX^T MXj
989  // = Xj - prevX product
990  {
991 #ifdef BELOS_TEUCHOS_TIME_MONITOR
992  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
993 #endif
994  MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
995  }
996 
997  // Update MXj
998  if (this->_hasOp) {
999  // MXj <- Op*Xj_new
1000  // = Op*(Xj_old - prevX prevX^T MXj)
1001  // = MXj - prevMX product
1002 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1003  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1004 #endif
1005  MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1006  }
1007 
1008  // Set coefficients
1009  if ( i==0 )
1010  product[ii] = P2[0];
1011  else
1012  product[ii] += P2[0];
1013 
1014  } // for (int i=0; i<max_ortho_steps_; ++i)
1015 
1016  } // for (int ii=0; ii<numX; ++ii)
1017 
1018  // Compute Op-norm with old MXj
1019  if (numX > 0) {
1020 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1021  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1022 #endif
1023  MVT::MvDot( *Xj, *oldMXj, newDot );
1024  }
1025  else {
1026  newDot[0] = oldDot[0];
1027  }
1028 
1029  // Check to see if the new vector is dependent.
1030  if (completeBasis) {
1031  //
1032  // We need a complete basis, so add random vectors if necessary
1033  //
1034  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1035 
1036  // Add a random vector and orthogonalize it against previous vectors in block.
1037  addVec = true;
1038 #ifdef ORTHO_DEBUG
1039  std::cout << "Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1040 #endif
1041  //
1042  Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1043  Teuchos::RCP<MV> tempMXj;
1044  MVT::MvRandom( *tempXj );
1045  if (this->_hasOp) {
1046  tempMXj = MVT::Clone( X, 1 );
1047  OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1048  }
1049  else {
1050  tempMXj = tempXj;
1051  }
1052  {
1053 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1054  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1055 #endif
1056  MVT::MvDot( *tempXj, *tempMXj, oldDot );
1057  }
1058  //
1059  // Perform MGS one vector at a time
1060  for (int ii=0; ii<numX; ii++) {
1061 
1062  index[0] = ii;
1063  prevX = MVT::CloneView( X, index );
1064  if (this->_hasOp) {
1065  prevMX = MVT::CloneView( *MX, index );
1066  }
1067 
1068  for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1069  {
1070 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1071  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1072 #endif
1073  MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,P2);
1074  }
1075  {
1076 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1077  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1078 #endif
1079  MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
1080  }
1081  if (this->_hasOp) {
1082 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1083  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1084 #endif
1085  MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
1086  }
1087 
1088  // Set coefficients
1089  if ( num_orth==0 )
1090  product[ii] = P2[0];
1091  else
1092  product[ii] += P2[0];
1093  }
1094  }
1095 
1096  // Compute new Op-norm
1097  {
1098 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1099  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1100 #endif
1101  MVT::MvDot( *tempXj, *tempMXj, newDot );
1102  }
1103  //
1104  if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1105  // Copy vector into current column of _basisvecs
1106  MVT::Assign( *tempXj, *Xj );
1107  if (this->_hasOp) {
1108  MVT::Assign( *tempMXj, *MXj );
1109  }
1110  }
1111  else {
1112  return numX;
1113  }
1114  }
1115  }
1116  else {
1117  //
1118  // We only need to detect dependencies.
1119  //
1120  if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1121  return numX;
1122  }
1123  }
1124 
1125 
1126  // If we haven't left this method yet, then we can normalize the new vector Xj.
1127  // Normalize Xj.
1128  // Xj <- Xj / std::sqrt(newDot)
1129  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1130  if (SCT::magnitude(diag) > ZERO) {
1131  MVT::MvScale( *Xj, ONE/diag );
1132  if (this->_hasOp) {
1133  // Update MXj.
1134  MVT::MvScale( *MXj, ONE/diag );
1135  }
1136  }
1137 
1138  // If we've added a random vector, enter a zero in the j'th diagonal element.
1139  if (addVec) {
1140  (*B)(j,j) = ZERO;
1141  }
1142  else {
1143  (*B)(j,j) = diag;
1144  }
1145 
1146  // Save the coefficients, if we are working on the original vector and not a randomly generated one
1147  if (!addVec) {
1148  for (int i=0; i<numX; i++) {
1149  (*B)(i,j) = product(i);
1150  }
1151  }
1152 
1153  } // for (j = 0; j < xc; ++j)
1154 
1155  return xc;
1156  }
1157 
1159  // Routine to compute the block orthogonalization
1160  template<class ScalarType, class MV, class OP>
1161  bool
1162  IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1163  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1164  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1165  {
1166  int nq = Q.size();
1167  int xc = MVT::GetNumberVecs( X );
1168  const ScalarType ONE = SCT::one();
1169 
1170  std::vector<int> qcs( nq );
1171  for (int i=0; i<nq; i++) {
1172  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1173  }
1174 
1175  // Perform the Gram-Schmidt transformation for a block of vectors
1176  std::vector<int> index(1);
1177  Teuchos::RCP<const MV> tempQ;
1178 
1179  Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1180  // Define the product Q^T * (Op*X)
1181  for (int i=0; i<nq; i++) {
1182 
1183  // Perform MGS one vector at a time
1184  for (int ii=0; ii<qcs[i]; ii++) {
1185 
1186  index[0] = ii;
1187  tempQ = MVT::CloneView( *Q[i], index );
1188  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1189 
1190  // Multiply Q' with MX
1191  {
1192 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1193  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1194 #endif
1195  MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,X,MX,tempC);
1196  }
1197  // Multiply by Q and subtract the result in X
1198  {
1199 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1200  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1201 #endif
1202  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1203  }
1204  }
1205 
1206  // Update MX, with the least number of applications of Op as possible
1207  if (this->_hasOp) {
1208  if (xc <= qcs[i]) {
1209  OPT::Apply( *(this->_Op), X, *MX);
1210  }
1211  else {
1212  // this will possibly be used again below; don't delete it
1213  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1214  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1215  MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1216  }
1217  }
1218  }
1219 
1220  // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1221  for (int j = 1; j < max_ortho_steps_; ++j) {
1222 
1223  for (int i=0; i<nq; i++) {
1224 
1225  Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1226 
1227  // Perform MGS one vector at a time
1228  for (int ii=0; ii<qcs[i]; ii++) {
1229 
1230  index[0] = ii;
1231  tempQ = MVT::CloneView( *Q[i], index );
1232  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1233  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1234 
1235  // Apply another step of modified Gram-Schmidt
1236  {
1237 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1238  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1239 #endif
1240  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1241  }
1242  tempC += tempC2;
1243  {
1244 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1245  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1246 #endif
1247  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1248  }
1249 
1250  }
1251 
1252  // Update MX, with the least number of applications of Op as possible
1253  if (this->_hasOp) {
1254  if (MQ[i].get()) {
1255  // MQ was allocated and computed above; use it
1256  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1257  }
1258  else if (xc <= qcs[i]) {
1259  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1260  OPT::Apply( *(this->_Op), X, *MX);
1261  }
1262  }
1263  } // for (int i=0; i<nq; i++)
1264  } // for (int j = 0; j < max_ortho_steps; ++j)
1265 
1266  return false;
1267  }
1268 
1270  // Routine to compute the block orthogonalization
1271  template<class ScalarType, class MV, class OP>
1272  bool
1273  IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1274  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1275  Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1276  {
1277  int nq = Q.size();
1278  int xc = MVT::GetNumberVecs( X );
1279  bool dep_flg = false;
1280  const ScalarType ONE = SCT::one();
1281 
1282  std::vector<int> qcs( nq );
1283  for (int i=0; i<nq; i++) {
1284  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1285  }
1286 
1287  // Perform the Gram-Schmidt transformation for a block of vectors
1288 
1289  // Compute the initial Op-norms
1290  std::vector<ScalarType> oldDot( xc );
1291  {
1292 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1293  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1294 #endif
1295  MVT::MvDot( X, *MX, oldDot );
1296  }
1297 
1298  std::vector<int> index(1);
1299  Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1300  Teuchos::RCP<const MV> tempQ;
1301 
1302  // Define the product Q^T * (Op*X)
1303  for (int i=0; i<nq; i++) {
1304 
1305  // Perform MGS one vector at a time
1306  for (int ii=0; ii<qcs[i]; ii++) {
1307 
1308  index[0] = ii;
1309  tempQ = MVT::CloneView( *Q[i], index );
1310  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1311 
1312  // Multiply Q' with MX
1313  {
1314 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1315  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1316 #endif
1317  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC);
1318  }
1319  // Multiply by Q and subtract the result in X
1320  {
1321 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1322  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1323 #endif
1324  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1325  }
1326  }
1327 
1328  // Update MX, with the least number of applications of Op as possible
1329  if (this->_hasOp) {
1330  if (xc <= qcs[i]) {
1331  OPT::Apply( *(this->_Op), X, *MX);
1332  }
1333  else {
1334  // this will possibly be used again below; don't delete it
1335  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1336  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1337  MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1338  }
1339  }
1340  }
1341 
1342  // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1343  for (int j = 1; j < max_ortho_steps_; ++j) {
1344 
1345  for (int i=0; i<nq; i++) {
1346  Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1347 
1348  // Perform MGS one vector at a time
1349  for (int ii=0; ii<qcs[i]; ii++) {
1350 
1351  index[0] = ii;
1352  tempQ = MVT::CloneView( *Q[i], index );
1353  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1354  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1355 
1356  // Apply another step of modified Gram-Schmidt
1357  {
1358 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1359  Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1360 #endif
1361  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1362  }
1363  tempC += tempC2;
1364  {
1365 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1366  Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1367 #endif
1368  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1369  }
1370  }
1371 
1372  // Update MX, with the least number of applications of Op as possible
1373  if (this->_hasOp) {
1374  if (MQ[i].get()) {
1375  // MQ was allocated and computed above; use it
1376  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1377  }
1378  else if (xc <= qcs[i]) {
1379  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1380  OPT::Apply( *(this->_Op), X, *MX);
1381  }
1382  }
1383  } // for (int i=0; i<nq; i++)
1384  } // for (int j = 0; j < max_ortho_steps; ++j)
1385 
1386  // Compute new Op-norms
1387  std::vector<ScalarType> newDot(xc);
1388  {
1389 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1390  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1391 #endif
1392  MVT::MvDot( X, *MX, newDot );
1393  }
1394 
1395  // Check to make sure the new block of vectors are not dependent on previous vectors
1396  for (int i=0; i<xc; i++){
1397  if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1398  dep_flg = true;
1399  break;
1400  }
1401  } // end for (i=0;...)
1402 
1403  return dep_flg;
1404  }
1405 
1406  template<class ScalarType, class MV, class OP>
1407  int
1408  IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1409  Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1410  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1411  Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
1412  {
1413  Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1414 
1415  const ScalarType ONE = SCT::one();
1416  const ScalarType ZERO = SCT::zero();
1417 
1418  int nq = Q.size();
1419  int xc = MVT::GetNumberVecs( X );
1420  std::vector<int> indX( 1 );
1421  std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1422 
1423  std::vector<int> qcs( nq );
1424  for (int i=0; i<nq; i++) {
1425  qcs[i] = MVT::GetNumberVecs( *Q[i] );
1426  }
1427 
1428  // Create pointers for the previous vectors of X that have already been orthonormalized.
1429  Teuchos::RCP<const MV> lastQ;
1430  Teuchos::RCP<MV> Xj, MXj;
1431  Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1432 
1433  // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
1434  for (int j=0; j<xc; j++) {
1435 
1436  bool dep_flg = false;
1437 
1438  // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
1439  if (j > 0) {
1440  std::vector<int> index( j );
1441  for (int ind=0; ind<j; ind++) {
1442  index[ind] = ind;
1443  }
1444  lastQ = MVT::CloneView( X, index );
1445 
1446  // Add these views to the Q and C arrays.
1447  Q.push_back( lastQ );
1448  C.push_back( B );
1449  qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1450  }
1451 
1452  // Get a view of the current vector in X to orthogonalize.
1453  indX[0] = j;
1454  Xj = MVT::CloneViewNonConst( X, indX );
1455  if (this->_hasOp) {
1456  MXj = MVT::CloneViewNonConst( *MX, indX );
1457  }
1458  else {
1459  MXj = Xj;
1460  }
1461 
1462  // Compute the initial Op-norms
1463  {
1464 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1465  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1466 #endif
1467  MVT::MvDot( *Xj, *MXj, oldDot );
1468  }
1469 
1470  Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1471  Teuchos::RCP<const MV> tempQ;
1472 
1473  // Define the product Q^T * (Op*X)
1474  for (int i=0; i<Q.size(); i++) {
1475 
1476  // Perform MGS one vector at a time
1477  for (int ii=0; ii<qcs[i]; ii++) {
1478 
1479  indX[0] = ii;
1480  tempQ = MVT::CloneView( *Q[i], indX );
1481  // Get a view of the current serial dense matrix
1482  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1483 
1484  // Multiply Q' with MX
1485  MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,*Xj,MXj,tempC);
1486 
1487  // Multiply by Q and subtract the result in Xj
1488  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1489  }
1490 
1491  // Update MXj, with the least number of applications of Op as possible
1492  if (this->_hasOp) {
1493  if (xc <= qcs[i]) {
1494  OPT::Apply( *(this->_Op), *Xj, *MXj);
1495  }
1496  else {
1497  // this will possibly be used again below; don't delete it
1498  MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1499  OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1500  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1501  MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1502  }
1503  }
1504  }
1505 
1506  // Do any additional steps of modified Gram-Schmidt orthogonalization
1507  for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1508 
1509  for (int i=0; i<Q.size(); i++) {
1510  Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1511 
1512  // Perform MGS one vector at a time
1513  for (int ii=0; ii<qcs[i]; ii++) {
1514 
1515  indX[0] = ii;
1516  tempQ = MVT::CloneView( *Q[i], indX );
1517  // Get a view of the current serial dense matrix
1518  Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1519 
1520  // Apply another step of modified Gram-Schmidt
1521  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *Xj, MXj, tempC2);
1522  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1523  }
1524 
1525  // Add the coefficients into C[i]
1526  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1527  tempC += C2;
1528 
1529  // Update MXj, with the least number of applications of Op as possible
1530  if (this->_hasOp) {
1531  if (MQ[i].get()) {
1532  // MQ was allocated and computed above; use it
1533  MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1534  }
1535  else if (xc <= qcs[i]) {
1536  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1537  OPT::Apply( *(this->_Op), *Xj, *MXj);
1538  }
1539  }
1540  } // for (int i=0; i<Q.size(); i++)
1541 
1542  } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
1543 
1544  // Compute the Op-norms after the correction step.
1545  {
1546 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1547  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1548 #endif
1549  MVT::MvDot( *Xj, *MXj, newDot );
1550  }
1551 
1552  // Check for linear dependence.
1553  if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1554  dep_flg = true;
1555  }
1556 
1557  // Normalize the new vector if it's not dependent
1558  if (!dep_flg) {
1559  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1560 
1561  MVT::MvScale( *Xj, ONE/diag );
1562  if (this->_hasOp) {
1563  // Update MXj.
1564  MVT::MvScale( *MXj, ONE/diag );
1565  }
1566 
1567  // Enter value on diagonal of B.
1568  (*B)(j,j) = diag;
1569  }
1570  else {
1571  // Create a random vector and orthogonalize it against all previous columns of Q.
1572  Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1573  Teuchos::RCP<MV> tempMXj;
1574  MVT::MvRandom( *tempXj );
1575  if (this->_hasOp) {
1576  tempMXj = MVT::Clone( X, 1 );
1577  OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1578  }
1579  else {
1580  tempMXj = tempXj;
1581  }
1582  {
1583 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1584  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1585 #endif
1586  MVT::MvDot( *tempXj, *tempMXj, oldDot );
1587  }
1588  //
1589  for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1590 
1591  for (int i=0; i<Q.size(); i++) {
1592  Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1593 
1594  // Perform MGS one vector at a time
1595  for (int ii=0; ii<qcs[i]; ii++) {
1596 
1597  indX[0] = ii;
1598  tempQ = MVT::CloneView( *Q[i], indX );
1599  Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1600 
1601  // Apply another step of modified Gram-Schmidt
1602  MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *tempXj, tempMXj, tempC );
1603  MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1604 
1605  }
1606 
1607  // Update MXj, with the least number of applications of Op as possible
1608  if (this->_hasOp) {
1609  if (MQ[i].get()) {
1610  // MQ was allocated and computed above; use it
1611  MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1612  }
1613  else if (xc <= qcs[i]) {
1614  // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1615  OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1616  }
1617  }
1618  } // for (int i=0; i<nq; i++)
1619  } // for (int num_orth=0; num_orth<max_orth_steps_; num_orth++)
1620 
1621  // Compute the Op-norms after the correction step.
1622  {
1623 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1624  Teuchos::TimeMonitor normTimer( *timerNorm_ );
1625 #endif
1626  MVT::MvDot( *tempXj, *tempMXj, newDot );
1627  }
1628 
1629  // Copy vector into current column of Xj
1630  if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1631  ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1632 
1633  // Enter value on diagonal of B.
1634  (*B)(j,j) = ZERO;
1635 
1636  // Copy vector into current column of _basisvecs
1637  MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1638  if (this->_hasOp) {
1639  MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1640  }
1641  }
1642  else {
1643  return j;
1644  }
1645  } // if (!dep_flg)
1646 
1647  // Remove the vectors from array
1648  if (j > 0) {
1649  Q.resize( nq );
1650  C.resize( nq );
1651  qcs.resize( nq );
1652  }
1653 
1654  } // for (int j=0; j<xc; j++)
1655 
1656  return xc;
1657  }
1658 
1659  template<class ScalarType, class MV, class OP>
1660  Teuchos::RCP<Teuchos::ParameterList> getIMGSDefaultParameters ()
1661  {
1662  using Teuchos::ParameterList;
1663  using Teuchos::parameterList;
1664  using Teuchos::RCP;
1665 
1666  RCP<ParameterList> params = parameterList ("IMGS");
1667 
1668  // Default parameter values for IMGS orthogonalization.
1669  // Documentation will be embedded in the parameter list.
1670  params->set ("maxNumOrthogPasses", IMGSOrthoManager<ScalarType, MV, OP>::max_ortho_steps_default_,
1671  "Maximum number of orthogonalization passes (includes the "
1672  "first). Default is 2, since \"twice is enough\" for Krylov "
1673  "methods.");
1674  params->set ("blkTol", IMGSOrthoManager<ScalarType, MV, OP>::blk_tol_default_,
1675  "Block reorthogonalization threshold.");
1676  params->set ("singTol", IMGSOrthoManager<ScalarType, MV, OP>::sing_tol_default_,
1677  "Singular block detection threshold.");
1678 
1679  return params;
1680  }
1681 
1682  template<class ScalarType, class MV, class OP>
1683  Teuchos::RCP<Teuchos::ParameterList> getIMGSFastParameters ()
1684  {
1685  using Teuchos::ParameterList;
1686  using Teuchos::RCP;
1687 
1688  RCP<ParameterList> params = getIMGSDefaultParameters<ScalarType, MV, OP>();
1689 
1690  params->set ("maxNumOrthogPasses",
1691  IMGSOrthoManager<ScalarType, MV, OP>::max_ortho_steps_fast_);
1692  params->set ("blkTol",
1693  IMGSOrthoManager<ScalarType, MV, OP>::blk_tol_fast_);
1694  params->set ("singTol",
1695  IMGSOrthoManager<ScalarType, MV, OP>::sing_tol_fast_);
1696 
1697  return params;
1698  }
1699 
1700 } // namespace Belos
1701 
1702 #endif // BELOS_IMGS_ORTHOMANAGER_HPP
1703 
Belos::IMGSOrthoManager::setParameterList
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
Definition: BelosIMGSOrthoManager.hpp:131
Teuchos_SerialDenseMatrix.hpp
is_null
bool is_null(const boost::shared_ptr< T > &p)
Belos::IMGSOrthoManager::orthogError
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...
Definition: BelosIMGSOrthoManager.hpp:384
Belos::IMGSOrthoManager::max_ortho_steps_fast_
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
Definition: BelosIMGSOrthoManager.hpp:421
Belos::IMGSOrthoManager::setBlkTol
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
Definition: BelosIMGSOrthoManager.hpp:218
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Belos::IMGSOrthoManager::project
void project(MV &X, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
This method calls project(X,Teuchos::null,C,Q); see documentation for that function.
Definition: BelosIMGSOrthoManager.hpp:269
Teuchos::TimeMonitor::getNewCounter
static RCP< Time > getNewCounter(const std::string &name)
Teuchos::is_null
bool is_null(const std::shared_ptr< T > &p)
Belos::IMGSOrthoManager::orthonormError
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector, measured as the Frobenius norm of ...
Definition: BelosIMGSOrthoManager.hpp:369
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
BelosMultiVecTraits.hpp
Declaration of basic traits for the multivector type.
BelosOperatorTraits.hpp
Class which defines basic traits for the operator type.
Teuchos::SerialDenseMatrix::normFrobenius
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Belos::IMGSOrthoManager::getLabel
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
Definition: BelosIMGSOrthoManager.hpp:406
Belos::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition: BelosMultiVecTraits.hpp:97
Belos::MatOrthoManager::innerProd
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Provides the inner product defining the orthogonality concepts, using the provided operator...
Definition: BelosMatOrthoManager.hpp:83
Belos::IMGSOrthoManager::IMGSOrthoManager
IMGSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_ortho_steps=max_ortho_steps_default_, const MagnitudeType blk_tol=blk_tol_default_, const MagnitudeType sing_tol=sing_tol_default_)
Constructor specifying re-orthogonalization tolerance.
Definition: BelosIMGSOrthoManager.hpp:65
Teuchos::rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Belos::IMGSOrthoManager::normalize
int normalize(MV &X, Teuchos::RCP< MV > MX, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method takes a multivector X and attempts to compute an orthonormal basis for ...
Definition: BelosIMGSOrthoManager.hpp:749
Belos::IMGSOrthoManager::getBlkTol
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
Definition: BelosIMGSOrthoManager.hpp:224
Teuchos::ParameterListAcceptorDefaultBase::setMyParamList
void setMyParamList(const RCP< ParameterList > &paramList)
Belos::IMGSOrthoManager::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition: BelosIMGSOrthoManager.hpp:203
Belos::IMGSOrthoManager::setLabel
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
Definition: BelosIMGSOrthoManager.hpp:512
Belos::IMGSOrthoManager
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
Definition: BelosIMGSOrthoManager.hpp:50
Belos::IMGSOrthoManager::project
void project(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a list of (mutually and internally) orthonormal bases Q, this method takes a multivector X and ...
Definition: BelosIMGSOrthoManager.hpp:765
Belos::IMGSOrthoManager::sing_tol_default_
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
Definition: BelosIMGSOrthoManager.hpp:418
Belos::getIMGSDefaultParameters
Teuchos::RCP< Teuchos::ParameterList > getIMGSDefaultParameters()
&quot;Default&quot; parameters for robustness and accuracy.
Definition: BelosIMGSOrthoManager.hpp:1660
Teuchos_SerialDenseVector.hpp
Belos::IMGSOrthoManager::max_ortho_steps_default_
static const int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
Definition: BelosIMGSOrthoManager.hpp:414
Belos::IMGSOrthoManager::IMGSOrthoManager
IMGSOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &plist, const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor that takes a list of parameters.
Definition: BelosIMGSOrthoManager.hpp:95
Belos::IMGSOrthoManager::getSingTol
MagnitudeType getSingTol() const
Return parameter for singular block detection.
Definition: BelosIMGSOrthoManager.hpp:227
Belos::OperatorTraits
Class which defines basic traits for the operator type.
Definition: BelosOperatorTraits.hpp:77
Belos::IMGSOrthoManager::normalize
int normalize(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method calls normalize(X,Teuchos::null,B); see documentation for that function.
Definition: BelosIMGSOrthoManager.hpp:307
Teuchos_TimeMonitor.hpp
Belos::IMGSOrthoManager::projectAndNormalizeWithMxImpl
virtual int projectAndNormalizeWithMxImpl(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for ...
Definition: BelosIMGSOrthoManager.hpp:567
Belos::MatOrthoManager
Belos&#39;s templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Definition: BelosIteration.hpp:38
Belos::IMGSOrthoManager::setSingTol
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
Definition: BelosIMGSOrthoManager.hpp:221
BelosConfigDefs.hpp
Belos header file which uses auto-configuration information to include necessary C++ headers...
Belos::getIMGSFastParameters
Teuchos::RCP< Teuchos::ParameterList > getIMGSFastParameters()
&quot;Fast&quot; but possibly unsafe or less accurate parameters.
Definition: BelosIMGSOrthoManager.hpp:1683
Teuchos_as.hpp
Belos::IMGSOrthoManager::blk_tol_default_
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
Definition: BelosIMGSOrthoManager.hpp:416
Belos::IMGSOrthoManager::sing_tol_fast_
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
Definition: BelosIMGSOrthoManager.hpp:425
Belos::IMGSOrthoManager::blk_tol_fast_
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
Definition: BelosIMGSOrthoManager.hpp:423
BelosMatOrthoManager.hpp
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::RCP::is_null
bool is_null() const
Generated for Belos by doxygen 1.8.5