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