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