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