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