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