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