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