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