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