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