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 "BelosDenseMatTraits.hpp"
30#include "BelosOperatorTraits.hpp"
31#include "BelosMatOrthoManager.hpp"
32#include "BelosTeuchosDenseAdapter.hpp"
33#include "BelosKokkosDenseAdapter.hpp"
34
35#include "Teuchos_as.hpp"
36#ifdef BELOS_TEUCHOS_TIME_MONITOR
37#include "Teuchos_TimeMonitor.hpp"
38#endif // BELOS_TEUCHOS_TIME_MONITOR
39
40namespace Belos {
41
43 template<class ScalarType, class MV, class OP, class DM = DefaultDenseMatrix<int,ScalarType>>
44 Teuchos::RCP<Teuchos::ParameterList> getIMGSDefaultParameters ();
45
47 template<class ScalarType, class MV, class OP, class DM = DefaultDenseMatrix<int,ScalarType>>
48 Teuchos::RCP<Teuchos::ParameterList> getIMGSFastParameters();
49
50 template<class ScalarType, class MV, class OP, class DM = DefaultDenseMatrix<int,ScalarType>>
51 class IMGSOrthoManager :
52 public MatOrthoManager<ScalarType,MV,OP,DM>
53 {
54 private:
55 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
56 typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
57 typedef Teuchos::ScalarTraits<ScalarType> SCT;
58 typedef MultiVecTraits<ScalarType,MV,DM> MVT;
59 typedef DenseMatTraits<ScalarType,DM> DMT;
60 typedef OperatorTraits<ScalarType,MV,OP> OPT;
61
62 public:
64
65
67 IMGSOrthoManager( const std::string& label = "Belos",
68 Teuchos::RCP<const OP> Op = Teuchos::null,
69 const int max_ortho_steps = max_ortho_steps_default_,
70 const MagnitudeType blk_tol = blk_tol_default_,
71 const MagnitudeType sing_tol = sing_tol_default_ )
72 : MatOrthoManager<ScalarType,MV,OP,DM>(Op),
73 max_ortho_steps_( max_ortho_steps ),
74 blk_tol_( blk_tol ),
75 sing_tol_( sing_tol ),
76 label_( label )
77 {
78#ifdef BELOS_TEUCHOS_TIME_MONITOR
79 std::stringstream ss;
80 ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
81
82 std::string orthoLabel = ss.str() + ": Orthogonalization";
83 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
84
85 std::string updateLabel = ss.str() + ": Ortho (Update)";
86 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
87
88 std::string normLabel = ss.str() + ": Ortho (Norm)";
89 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
90
91 std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
92 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
93#endif
94 }
95
97 IMGSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
98 const std::string& label = "Belos",
99 Teuchos::RCP<const OP> Op = Teuchos::null) :
100 MatOrthoManager<ScalarType,MV,OP,DM>(Op),
101 max_ortho_steps_ (max_ortho_steps_default_),
102 blk_tol_ (blk_tol_default_),
103 sing_tol_ (sing_tol_default_),
104 label_ (label)
105 {
106 setParameterList (plist);
107
108#ifdef BELOS_TEUCHOS_TIME_MONITOR
109 std::stringstream ss;
110 ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
111
112 std::string orthoLabel = ss.str() + ": Orthogonalization";
113 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
114
115 std::string updateLabel = ss.str() + ": Ortho (Update)";
116 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
117
118 std::string normLabel = ss.str() + ": Ortho (Norm)";
119 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
120
121 std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
122 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
123#endif
124 }
125
127 ~IMGSOrthoManager() {}
129
131
132 void
133 setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
134 {
135 using Teuchos::Exceptions::InvalidParameterName;
136 using Teuchos::ParameterList;
137 using Teuchos::parameterList;
138 using Teuchos::RCP;
139
140 RCP<const ParameterList> defaultParams = getValidParameters();
141 RCP<ParameterList> params;
142 if (plist.is_null()) {
143 params = parameterList (*defaultParams);
144 } else {
145 params = plist;
146 // Some users might want to specify "blkTol" as "depTol". Due
147 // to this case, we don't invoke
148 // validateParametersAndSetDefaults on params. Instead, we go
149 // through the parameter list one parameter at a time and look
150 // for alternatives.
151 }
152
153 // Using temporary variables and fetching all values before
154 // setting the output arguments ensures the strong exception
155 // guarantee for this function: if an exception is thrown, no
156 // externally visible side effects (in this case, setting the
157 // output arguments) have taken place.
158 int maxNumOrthogPasses;
159 MagnitudeType blkTol;
160 MagnitudeType singTol;
161
162 try {
163 maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
164 } catch (InvalidParameterName&) {
165 maxNumOrthogPasses = defaultParams->get<int> ("maxNumOrthogPasses");
166 params->set ("maxNumOrthogPasses", maxNumOrthogPasses);
167 }
168
169 // Handling of the "blkTol" parameter is a special case. This
170 // is because some users may prefer to call this parameter
171 // "depTol" for consistency with DGKS. However, our default
172 // parameter list calls this "blkTol", and we don't want the
173 // default list's value to override the user's value. Thus, we
174 // first check the user's parameter list for both names, and
175 // only then access the default parameter list.
176 try {
177 blkTol = params->get<MagnitudeType> ("blkTol");
178 } catch (InvalidParameterName&) {
179 try {
180 blkTol = params->get<MagnitudeType> ("depTol");
181 // "depTol" is the wrong name, so remove it and replace with
182 // "blkTol". We'll set "blkTol" below.
183 params->remove ("depTol");
184 } catch (InvalidParameterName&) {
185 blkTol = defaultParams->get<MagnitudeType> ("blkTol");
186 }
187 params->set ("blkTol", blkTol);
188 }
189
190 try {
191 singTol = params->get<MagnitudeType> ("singTol");
192 } catch (InvalidParameterName&) {
193 singTol = defaultParams->get<MagnitudeType> ("singTol");
194 params->set ("singTol", singTol);
195 }
196
197 max_ortho_steps_ = maxNumOrthogPasses;
198 blk_tol_ = blkTol;
199 sing_tol_ = singTol;
200
201 this->setMyParamList (params);
202 }
203
204 Teuchos::RCP<const Teuchos::ParameterList>
205 getValidParameters () const
206 {
207 if (defaultParams_.is_null()) {
208 defaultParams_ = Belos::getIMGSDefaultParameters<ScalarType, MV, OP>();
209 }
210
211 return defaultParams_;
212 }
213
215
217
218
220 void setBlkTol( const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
221
223 void setSingTol( const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
224
226 MagnitudeType getBlkTol() const { return blk_tol_; }
227
229 MagnitudeType getSingTol() const { return sing_tol_; }
230
232
233
235
236
264 void project ( MV &X, Teuchos::RCP<MV> MX,
265 Teuchos::Array<Teuchos::RCP<DM> > C,
266 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
267
268
271 void project ( MV &X,
272 Teuchos::Array<Teuchos::RCP<DM> > C,
273 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
274 project(X,Teuchos::null,C,Q);
275 }
276
277
278
303 int normalize ( MV &X, Teuchos::RCP<MV> MX,
304 Teuchos::RCP<DM> B) const;
305
306
309 int normalize ( MV &X, Teuchos::RCP<DM> B ) const {
310 return normalize(X,Teuchos::null,B);
311 }
312
313 protected:
355 virtual int
356 projectAndNormalizeWithMxImpl (MV &X,
357 Teuchos::RCP<MV> MX,
358 Teuchos::Array<Teuchos::RCP<DM>> C,
359 Teuchos::RCP<DM> B,
360 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
361
362 public:
364
366
370 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
371 orthonormError(const MV &X) const {
372 return orthonormError(X,Teuchos::null);
373 }
374
379 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
380 orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
381
385 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
386 orthogError(const MV &X1, const MV &X2) const {
387 return orthogError(X1,Teuchos::null,X2);
388 }
389
394 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
395 orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
396
398
400
401
404 void setLabel(const std::string& label);
405
408 const std::string& getLabel() const { return label_; }
409
411
413
414
416 static const int max_ortho_steps_default_;
418 static const MagnitudeType blk_tol_default_;
420 static const MagnitudeType sing_tol_default_;
421
423 static const int max_ortho_steps_fast_;
425 static const MagnitudeType blk_tol_fast_;
427 static const MagnitudeType sing_tol_fast_;
428
430
431 private:
432
434 int max_ortho_steps_;
436 MagnitudeType blk_tol_;
438 MagnitudeType sing_tol_;
439
441 std::string label_;
442#ifdef BELOS_TEUCHOS_TIME_MONITOR
443 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
444#endif // BELOS_TEUCHOS_TIME_MONITOR
445
447 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
448
450 int findBasis(MV &X, Teuchos::RCP<MV> MX,
451 Teuchos::RCP<DM> C,
452 bool completeBasis, int howMany = -1 ) const;
453
455 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
456 Teuchos::Array<Teuchos::RCP<DM> > C,
457 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
458
460 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
461 Teuchos::Array<Teuchos::RCP<DM> > C,
462 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
463
477 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
478 Teuchos::Array<Teuchos::RCP<DM> > C,
479 Teuchos::RCP<DM> B,
480 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;
481 };
482
483 // Set static variables.
484 template<class ScalarType, class MV, class OP, class DM>
485 const int IMGSOrthoManager<ScalarType,MV,OP,DM>::max_ortho_steps_default_ = 1;
486
487 template<class ScalarType, class MV, class OP, class DM>
488 const typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType
489 IMGSOrthoManager<ScalarType,MV,OP,DM>::blk_tol_default_
490 = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType>::squareroot(
491 Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType>::eps() );
492
493 template<class ScalarType, class MV, class OP, class DM>
494 const typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType
495 IMGSOrthoManager<ScalarType,MV,OP,DM>::sing_tol_default_
496 = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType>::eps();
497
498 template<class ScalarType, class MV, class OP, class DM>
499 const int IMGSOrthoManager<ScalarType,MV,OP,DM>::max_ortho_steps_fast_ = 1;
500
501 template<class ScalarType, class MV, class OP, class DM>
502 const typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType
503 IMGSOrthoManager<ScalarType,MV,OP,DM>::blk_tol_fast_
504 = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType>::zero();
505
506 template<class ScalarType, class MV, class OP, class DM>
507 const typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType
508 IMGSOrthoManager<ScalarType,MV,OP,DM>::sing_tol_fast_
509 = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP,DM>::MagnitudeType>::zero();
510
512 // Set the label for this orthogonalization manager and create new timers if it's changed
513 template<class ScalarType, class MV, class OP, class DM>
514 void IMGSOrthoManager<ScalarType,MV,OP,DM>::setLabel(const std::string& label)
515 {
516 if (label != label_) {
517 label_ = label;
518#ifdef BELOS_TEUCHOS_TIME_MONITOR
519 std::stringstream ss;
520 ss << label_ + ": IMGS[" << max_ortho_steps_ << "]";
521
522 std::string orthoLabel = ss.str() + ": Orthogonalization";
523 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
524
525 std::string updateLabel = ss.str() + ": Ortho (Update)";
526 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
527
528 std::string normLabel = ss.str() + ": Ortho (Norm)";
529 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
530
531 std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
532 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
533#endif
534 }
535 }
536
538 // Compute the distance from orthonormality
539 template<class ScalarType, class MV, class OP, class DM>
540 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
541 IMGSOrthoManager<ScalarType,MV,OP,DM>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
542 const ScalarType ONE = SCT::one();
543 int rank = MVT::GetNumberVecs(X);
544 Teuchos::RCP<DM> xTx = DMT::Create(rank,rank);
545 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd(X,X,MX,*xTx);
546 DMT::SyncDeviceToHost(*xTx);
547 for (int i=0; i<rank; i++) {
548 DMT::Value(*xTx,i,i) -= ONE;
549 }
550 DMT::SyncHostToDevice( *xTx );
551 return DMT::NormFrobenius(*xTx);
552 }
553
555 // Compute the distance from orthogonality
556 template<class ScalarType, class MV, class OP, class DM>
557 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
558 IMGSOrthoManager<ScalarType,MV,OP,DM>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
559 int r1 = MVT::GetNumberVecs(X1);
560 int r2 = MVT::GetNumberVecs(X2);
561 Teuchos::RCP<DM> xTx = DMT::Create(r2,r1);
562 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd(X2,X1,MX1,*xTx);
563 return DMT::NormFrobenius(*xTx);
564 }
565
567 // Find an Op-orthonormal basis for span(X) - span(W)
568 template<class ScalarType, class MV, class OP, class DM>
569 int
570 IMGSOrthoManager<ScalarType, MV, OP, DM>::
571 projectAndNormalizeWithMxImpl(MV &X,
572 Teuchos::RCP<MV> MX,
573 Teuchos::Array<Teuchos::RCP<DM>> C,
574 Teuchos::RCP<DM> B,
575 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
576 {
577 using Teuchos::Array;
578 using Teuchos::null;
579 using Teuchos::is_null;
580 using Teuchos::RCP;
581 using Teuchos::rcp;
582 typedef typename Array< RCP< const MV > >::size_type size_type;
583
584#ifdef BELOS_TEUCHOS_TIME_MONITOR
585 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
586#endif
587
588 ScalarType ONE = SCT::one();
589 const MagnitudeType ZERO = MGT::zero();
590
591 int nq = Q.size();
592 int xc = MVT::GetNumberVecs( X );
593 ptrdiff_t xr = MVT::GetGlobalLength( X );
594 int rank = xc;
595
596 // If the user doesn't want to store the normalization
597 // coefficients, allocate some local memory for them. This will
598 // go away at the end of this method.
599 if (is_null (B)) {
600 B = DMT::Create(xc,xc);
601 }
602 // Likewise, if the user doesn't want to store the projection
603 // coefficients, allocate some local memory for them. Also make
604 // sure that all the entries of C are the right size. We're going
605 // to overwrite them anyway, so we don't have to worry about the
606 // contents (other than to resize them if they are the wrong
607 // size).
608 if (C.size() < nq)
609 C.resize (nq);
610 for (size_type k = 0; k < nq; ++k)
611 {
612 const int numRows = MVT::GetNumberVecs (*Q[k]);
613 const int numCols = xc; // Number of vectors in X
614
615 if (is_null (C[k]))
616 C[k] = DMT::Create(numRows,numCols);
617 else if (DMT::GetNumRows(*C[k]) != numRows || DMT::GetNumCols(*C[k]) != numCols)
618 {
619 DMT::Reshape(*C[k],numRows,numCols);
620 }
621 }
622
623 /****** DO NOT MODIFY *MX IF _hasOp == false ******/
624 if (this->_hasOp) {
625 if (MX == Teuchos::null) {
626 // we need to allocate space for MX
627 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
628 OPT::Apply(*(this->_Op),X,*MX);
629 }
630 }
631 else {
632 // Op == I --> MX = X (ignore it if the user passed it in)
633 MX = Teuchos::rcp( &X, false );
634 }
635
636 int mxc = MVT::GetNumberVecs( *MX );
637 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
638
639 // short-circuit
640 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
641
642 int numbas = 0;
643 for (int i=0; i<nq; i++) {
644 numbas += MVT::GetNumberVecs( *Q[i] );
645 }
646
647 // check size of B
648 TEUCHOS_TEST_FOR_EXCEPTION( DMT::GetNumRows(*B) != xc || DMT::GetNumCols(*B) != xc, std::invalid_argument,
649 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
650 // check size of X and MX
651 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
652 "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
653 // check size of X w.r.t. MX
654 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
655 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
656 // check feasibility
657 //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
658 // "Belos::IMGSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
659
660 // Some flags for checking dependency returns from the internal orthogonalization methods
661 bool dep_flg = false;
662
663 // Make a temporary copy of X and MX, just in case a block dependency is detected.
664 Teuchos::RCP<MV> tmpX, tmpMX;
665 tmpX = MVT::CloneCopy(X);
666 if (this->_hasOp) {
667 tmpMX = MVT::CloneCopy(*MX);
668 }
669
670 if (xc == 1) {
671
672 // Use the cheaper block orthogonalization.
673 // NOTE: Don't check for dependencies because the update has one vector.
674 dep_flg = blkOrtho1( X, MX, C, Q );
675
676 // Normalize the new block X
677 if ( B == Teuchos::null ) {
678 B = DMT::Create(xc,xc);
679 }
680 std::vector<ScalarType> dot(xc);
681 {
682#ifdef BELOS_TEUCHOS_TIME_MONITOR
683 Teuchos::TimeMonitor normTimer( *timerNorm_ );
684#endif
685 MVT::MvDot( X, *MX, dot );
686 }
687
688 ScalarType diag = SCT::squareroot(SCT::magnitude(dot[0]));
689 if (SCT::magnitude(diag) > ZERO) {
690 rank = 1;
691 MVT::MvScale( X, ONE/diag );
692 if (this->_hasOp) {
693 // Update MXj.
694 MVT::MvScale( *MX, ONE/diag );
695 }
696 }
697
698 Teuchos::RCP<DM> B00 = DMT::Subview(*B,1,1);
699 DMT::PutScalar(*B00, diag);
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, class DM>
749 int IMGSOrthoManager<ScalarType, MV, OP, DM>::normalize(
750 MV &X, Teuchos::RCP<MV> MX,
751 Teuchos::RCP<DM> 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, class DM>
765 void IMGSOrthoManager<ScalarType, MV, OP, DM>::project(
766 MV &X, Teuchos::RCP<MV> MX,
767 Teuchos::Array<Teuchos::RCP<DM>> 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] = DMT::Create(qcs[i],xc);
836 }
837 else {
838 TEUCHOS_TEST_FOR_EXCEPTION( DMT::GetNumRows(*C[i]) != qcs[i] || DMT::GetNumCols(*C[i]) != 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, class DM>
852 int IMGSOrthoManager<ScalarType, MV, OP, DM>::findBasis(
853 MV &X, Teuchos::RCP<MV> MX,
854 Teuchos::RCP<DM> 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 = DMT::Create(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( DMT::GetNumRows(*B) != xc || DMT::GetNumCols(*B) != 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::RCP<DM> product = DMT::Create(numX,1);
952 Teuchos::RCP<DM> P2 = DMT::Create(1,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,DM>::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 Teuchos::RCP<DM> product_ii = DMT::Subview(*product,1,1,ii,0);
1010 if ( i==0 )
1011 DMT::Assign(*product_ii, *P2);
1012 else
1013 DMT::Add(*product_ii, *P2);
1014
1015 } // for (int i=0; i<max_ortho_steps_; ++i)
1016
1017 } // for (int ii=0; ii<numX; ++ii)
1018
1019 // Compute Op-norm with old MXj
1020 if (numX > 0) {
1021#ifdef BELOS_TEUCHOS_TIME_MONITOR
1022 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1023#endif
1024 MVT::MvDot( *Xj, *oldMXj, newDot );
1025 }
1026 else {
1027 newDot[0] = oldDot[0];
1028 }
1029
1030 // Check to see if the new vector is dependent.
1031 if (completeBasis) {
1032 //
1033 // We need a complete basis, so add random vectors if necessary
1034 //
1035 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1036
1037 // Add a random vector and orthogonalize it against previous vectors in block.
1038 addVec = true;
1039#ifdef ORTHO_DEBUG
1040 std::cout << "Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1041#endif
1042 //
1043 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1044 Teuchos::RCP<MV> tempMXj;
1045 MVT::MvRandom( *tempXj );
1046 if (this->_hasOp) {
1047 tempMXj = MVT::Clone( X, 1 );
1048 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1049 }
1050 else {
1051 tempMXj = tempXj;
1052 }
1053 {
1054#ifdef BELOS_TEUCHOS_TIME_MONITOR
1055 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1056#endif
1057 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1058 }
1059 //
1060 // Perform MGS one vector at a time
1061 for (int ii=0; ii<numX; ii++) {
1062
1063 index[0] = ii;
1064 prevX = MVT::CloneView( X, index );
1065 if (this->_hasOp) {
1066 prevMX = MVT::CloneView( *MX, index );
1067 }
1068
1069 for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1070 {
1071#ifdef BELOS_TEUCHOS_TIME_MONITOR
1072 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1073#endif
1074 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd(*prevX,*tempXj,tempMXj,*P2);
1075 }
1076 {
1077#ifdef BELOS_TEUCHOS_TIME_MONITOR
1078 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1079#endif
1080 MVT::MvTimesMatAddMv( -ONE, *prevX, *P2, ONE, *tempXj );
1081 }
1082 if (this->_hasOp) {
1083#ifdef BELOS_TEUCHOS_TIME_MONITOR
1084 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1085#endif
1086 MVT::MvTimesMatAddMv( -ONE, *prevMX, *P2, ONE, *tempMXj );
1087 }
1088
1089 // Set coefficients
1090 Teuchos::RCP<DM> product_ii = DMT::Subview(*product,1,1,ii,0);
1091 if ( num_orth==0 )
1092 DMT::Assign(*product_ii, *P2);
1093 else
1094 DMT::Add(*product_ii, *P2);
1095 }
1096 }
1097
1098 // Compute new Op-norm
1099 {
1100#ifdef BELOS_TEUCHOS_TIME_MONITOR
1101 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1102#endif
1103 MVT::MvDot( *tempXj, *tempMXj, newDot );
1104 }
1105 //
1106 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1107 // Copy vector into current column of _basisvecs
1108 MVT::Assign( *tempXj, *Xj );
1109 if (this->_hasOp) {
1110 MVT::Assign( *tempMXj, *MXj );
1111 }
1112 }
1113 else {
1114 return numX;
1115 }
1116 }
1117 }
1118 else {
1119 //
1120 // We only need to detect dependencies.
1121 //
1122 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1123 return numX;
1124 }
1125 }
1126
1127
1128 // If we haven't left this method yet, then we can normalize the new vector Xj.
1129 // Normalize Xj.
1130 // Xj <- Xj / std::sqrt(newDot)
1131 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1132 if (SCT::magnitude(diag) > ZERO) {
1133 MVT::MvScale( *Xj, ONE/diag );
1134 if (this->_hasOp) {
1135 // Update MXj.
1136 MVT::MvScale( *MXj, ONE/diag );
1137 }
1138 }
1139
1140 // If we've added a random vector, enter a zero in the j'th diagonal element.
1141 Teuchos::RCP<DM> Bjj = DMT::Subview(*B,1,1,j,j);
1142 if (addVec) {
1143 DMT::PutScalar(*Bjj, ZERO);
1144 }
1145 else {
1146 DMT::PutScalar(*Bjj, diag);
1147 }
1148
1149 // Save the coefficients, if we are working on the original vector and not a randomly generated one
1150 if (!addVec) {
1151 Teuchos::RCP<DM> Bcolj = DMT::Subview(*B,numX,1,0,j);
1152 DMT::Assign(*Bcolj,*product);
1153 }
1154
1155 } // for (j = 0; j < xc; ++j)
1156
1157 return xc;
1158 }
1159
1161 // Routine to compute the block orthogonalization
1162 template<class ScalarType, class MV, class OP, class DM>
1163 bool
1164 IMGSOrthoManager<ScalarType, MV, OP, DM>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1165 Teuchos::Array<Teuchos::RCP<DM>> C,
1166 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1167 {
1168 int nq = Q.size();
1169 int xc = MVT::GetNumberVecs( X );
1170 const ScalarType ONE = SCT::one();
1171
1172 std::vector<int> qcs( nq );
1173 for (int i=0; i<nq; i++) {
1174 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1175 }
1176
1177 // Perform the Gram-Schmidt transformation for a block of vectors
1178 std::vector<int> index(1);
1179 Teuchos::RCP<const MV> tempQ;
1180
1181 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1182 // Define the product Q^T * (Op*X)
1183 for (int i=0; i<nq; i++) {
1184
1185 // Perform MGS one vector at a time
1186 for (int ii=0; ii<qcs[i]; ii++) {
1187
1188 index[0] = ii;
1189 tempQ = MVT::CloneView( *Q[i], index );
1190 Teuchos::RCP tempC = DMT::Subview(*C[i], 1, 1, ii, 0);
1191
1192 // Multiply Q' with MX
1193 {
1194#ifdef BELOS_TEUCHOS_TIME_MONITOR
1195 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1196#endif
1197 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd(*tempQ,X,MX,*tempC);
1198 }
1199 // Multiply by Q and subtract the result in X
1200 {
1201#ifdef BELOS_TEUCHOS_TIME_MONITOR
1202 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1203#endif
1204 MVT::MvTimesMatAddMv( -ONE, *tempQ, *tempC, ONE, X );
1205 }
1206 }
1207
1208 // Update MX, with the least number of applications of Op as possible
1209 if (this->_hasOp) {
1210 if (xc <= qcs[i]) {
1211 OPT::Apply( *(this->_Op), X, *MX);
1212 }
1213 else {
1214 // this will possibly be used again below; don't delete it
1215 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1216 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1217 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1218 }
1219 }
1220 }
1221
1222 // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1223 for (int j = 1; j < max_ortho_steps_; ++j) {
1224
1225 for (int i=0; i<nq; i++) {
1226
1227 Teuchos::RCP<DM> C2 = DMT::Create(qcs[i],1);
1228
1229 // Perform MGS one vector at a time
1230 for (int ii=0; ii<qcs[i]; ii++) {
1231
1232 index[0] = ii;
1233 tempQ = MVT::CloneView( *Q[i], index );
1234 Teuchos::RCP tempC = DMT::Subview(*C[i], 1, 1, ii, 0);
1235 Teuchos::RCP tempC2 = DMT::Subview(*C2, 1, 1, ii, 0);
1236
1237 // Apply another step of modified Gram-Schmidt
1238 {
1239#ifdef BELOS_TEUCHOS_TIME_MONITOR
1240 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1241#endif
1242 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd( *tempQ, X, MX, *tempC2 );
1243 }
1244 DMT::Add(*tempC,*tempC2);
1245 {
1246#ifdef BELOS_TEUCHOS_TIME_MONITOR
1247 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1248#endif
1249 MVT::MvTimesMatAddMv( -ONE, *tempQ, *tempC2, ONE, X );
1250 }
1251
1252 }
1253
1254 // Update MX, with the least number of applications of Op as possible
1255 if (this->_hasOp) {
1256 if (MQ[i].get()) {
1257 // MQ was allocated and computed above; use it
1258 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C2, ONE, *MX );
1259 }
1260 else if (xc <= qcs[i]) {
1261 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1262 OPT::Apply( *(this->_Op), X, *MX);
1263 }
1264 }
1265 } // for (int i=0; i<nq; i++)
1266 } // for (int j = 0; j < max_ortho_steps; ++j)
1267
1268 return false;
1269 }
1270
1272 // Routine to compute the block orthogonalization
1273 template<class ScalarType, class MV, class OP, class DM>
1274 bool
1275 IMGSOrthoManager<ScalarType, MV, OP, DM>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1276 Teuchos::Array<Teuchos::RCP<DM> > C,
1277 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1278 {
1279 int nq = Q.size();
1280 int xc = MVT::GetNumberVecs( X );
1281 bool dep_flg = false;
1282 const ScalarType ONE = SCT::one();
1283
1284 std::vector<int> qcs( nq );
1285 for (int i=0; i<nq; i++) {
1286 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1287 }
1288
1289 // Perform the Gram-Schmidt transformation for a block of vectors
1290
1291 // Compute the initial Op-norms
1292 std::vector<ScalarType> oldDot( xc );
1293 {
1294#ifdef BELOS_TEUCHOS_TIME_MONITOR
1295 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1296#endif
1297 MVT::MvDot( X, *MX, oldDot );
1298 }
1299
1300 std::vector<int> index(1);
1301 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1302 Teuchos::RCP<const MV> tempQ;
1303
1304 // Define the product Q^T * (Op*X)
1305 for (int i=0; i<nq; i++) {
1306
1307 // Perform MGS one vector at a time
1308 for (int ii=0; ii<qcs[i]; ii++) {
1309
1310 index[0] = ii;
1311 tempQ = MVT::CloneView( *Q[i], index );
1312 Teuchos::RCP tempC = DMT::Subview(*C[i], 1, xc, ii, 0);
1313
1314 // Multiply Q' with MX
1315 {
1316#ifdef BELOS_TEUCHOS_TIME_MONITOR
1317 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1318#endif
1319 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd( *tempQ, X, MX, *tempC);
1320 }
1321 // Multiply by Q and subtract the result in X
1322 {
1323#ifdef BELOS_TEUCHOS_TIME_MONITOR
1324 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1325#endif
1326 MVT::MvTimesMatAddMv( -ONE, *tempQ, *tempC, ONE, X );
1327 }
1328 }
1329
1330 // Update MX, with the least number of applications of Op as possible
1331 if (this->_hasOp) {
1332 if (xc <= qcs[i]) {
1333 OPT::Apply( *(this->_Op), X, *MX);
1334 }
1335 else {
1336 // this will possibly be used again below; don't delete it
1337 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1338 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1339 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1340 }
1341 }
1342 }
1343
1344 // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1345 for (int j = 1; j < max_ortho_steps_; ++j) {
1346
1347 for (int i=0; i<nq; i++) {
1348 Teuchos::RCP C2 = DMT::Create(qcs[i],xc);
1349
1350 // Perform MGS one vector at a time
1351 for (int ii=0; ii<qcs[i]; ii++) {
1352
1353 index[0] = ii;
1354 tempQ = MVT::CloneView( *Q[i], index );
1355 Teuchos::RCP tempC = DMT::Subview(*C[i], 1, xc, ii, 0);
1356 Teuchos::RCP tempC2 = DMT::Subview(*C2, 1, xc, ii, 0);
1357
1358 // Apply another step of modified Gram-Schmidt
1359 {
1360#ifdef BELOS_TEUCHOS_TIME_MONITOR
1361 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1362#endif
1363 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd( *tempQ, X, MX, *tempC2 );
1364 }
1365 DMT::Add(*tempC,*tempC2);
1366 {
1367#ifdef BELOS_TEUCHOS_TIME_MONITOR
1368 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1369#endif
1370 MVT::MvTimesMatAddMv( -ONE, *tempQ, *tempC2, ONE, X );
1371 }
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 // MQ was allocated and computed above; use it
1378 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C2, ONE, *MX );
1379 }
1380 else if (xc <= qcs[i]) {
1381 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1382 OPT::Apply( *(this->_Op), X, *MX);
1383 }
1384 }
1385 } // for (int i=0; i<nq; i++)
1386 } // for (int j = 0; j < max_ortho_steps; ++j)
1387
1388 // Compute new Op-norms
1389 std::vector<ScalarType> newDot(xc);
1390 {
1391#ifdef BELOS_TEUCHOS_TIME_MONITOR
1392 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1393#endif
1394 MVT::MvDot( X, *MX, newDot );
1395 }
1396
1397 // Check to make sure the new block of vectors are not dependent on previous vectors
1398 for (int i=0; i<xc; i++){
1399 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1400 dep_flg = true;
1401 break;
1402 }
1403 } // end for (i=0;...)
1404
1405 return dep_flg;
1406 }
1407
1408 template<class ScalarType, class MV, class OP, class DM>
1409 int
1410 IMGSOrthoManager<ScalarType, MV, OP, DM>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1411 Teuchos::Array<Teuchos::RCP<DM> > C,
1412 Teuchos::RCP<DM> B,
1413 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
1414 {
1415 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1416
1417 const ScalarType ONE = SCT::one();
1418 const ScalarType ZERO = SCT::zero();
1419
1420 int nq = Q.size();
1421 int xc = MVT::GetNumberVecs( X );
1422 std::vector<int> indX( 1 );
1423 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1424
1425 std::vector<int> qcs( nq );
1426 for (int i=0; i<nq; i++) {
1427 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1428 }
1429
1430 // Create pointers for the previous vectors of X that have already been orthonormalized.
1431 Teuchos::RCP<const MV> lastQ;
1432 Teuchos::RCP<MV> Xj, MXj;
1433
1434 // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
1435 for (int j=0; j<xc; j++) {
1436
1437 bool dep_flg = false;
1438
1439 // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
1440 if (j > 0) {
1441 std::vector<int> index( j );
1442 for (int ind=0; ind<j; ind++) {
1443 index[ind] = ind;
1444 }
1445 lastQ = MVT::CloneView( X, index );
1446
1447 // Add these views to the Q and C arrays.
1448 Q.push_back( lastQ );
1449 C.push_back( B );
1450 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1451 }
1452
1453 // Get a view of the current vector in X to orthogonalize.
1454 indX[0] = j;
1455 Xj = MVT::CloneViewNonConst( X, indX );
1456 if (this->_hasOp) {
1457 MXj = MVT::CloneViewNonConst( *MX, indX );
1458 }
1459 else {
1460 MXj = Xj;
1461 }
1462
1463 // Compute the initial Op-norms
1464 {
1465#ifdef BELOS_TEUCHOS_TIME_MONITOR
1466 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1467#endif
1468 MVT::MvDot( *Xj, *MXj, oldDot );
1469 }
1470
1471 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1472 Teuchos::RCP<const MV> tempQ;
1473
1474 // Define the product Q^T * (Op*X)
1475 for (int i=0; i<Q.size(); i++) {
1476
1477 // Perform MGS one vector at a time
1478 for (int ii=0; ii<qcs[i]; ii++) {
1479
1480 indX[0] = ii;
1481 tempQ = MVT::CloneView( *Q[i], indX );
1482 // Get a view of the current serial dense matrix
1483 Teuchos::RCP tempC = DMT::Subview(*C[i], 1, 1, ii, j);
1484
1485 // Multiply Q' with MX
1486 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd(*tempQ,*Xj,MXj,*tempC);
1487
1488 // Multiply by Q and subtract the result in Xj
1489 MVT::MvTimesMatAddMv( -ONE, *tempQ, *tempC, ONE, *Xj );
1490 }
1491
1492 // Update MXj, with the least number of applications of Op as possible
1493 if (this->_hasOp) {
1494 if (xc <= qcs[i]) {
1495 OPT::Apply( *(this->_Op), *Xj, *MXj);
1496 }
1497 else {
1498 // this will possibly be used again below; don't delete it
1499 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1500 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1501 Teuchos::RCP tempC = DMT::Subview(*C[i], qcs[i], 1, 0, j);
1502 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *tempC, ONE, *MXj );
1503 }
1504 }
1505 }
1506
1507 // Do any additional steps of modified Gram-Schmidt orthogonalization
1508 for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1509
1510 for (int i=0; i<Q.size(); i++) {
1511 Teuchos::RCP<DM> C2 = DMT::Create(qcs[i],1);
1512
1513 // Perform MGS one vector at a time
1514 for (int ii=0; ii<qcs[i]; ii++) {
1515
1516 indX[0] = ii;
1517 tempQ = MVT::CloneView( *Q[i], indX );
1518 // Get a view of the current serial dense matrix
1519 Teuchos::RCP tempC2 = DMT::Subview(*C2, 1, 1, ii);
1520
1521 // Apply another step of modified Gram-Schmidt
1522 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd( *tempQ, *Xj, MXj, *tempC2);
1523 MVT::MvTimesMatAddMv( -ONE, *tempQ, *tempC2, ONE, *Xj );
1524 }
1525
1526 // Add the coefficients into C[i]
1527 Teuchos::RCP tempC = DMT::Subview(*C[i], qcs[i], 1, 0, j);
1528 DMT::Add(*tempC,*C2);
1529
1530 // Update MXj, with the least number of applications of Op as possible
1531 if (this->_hasOp) {
1532 if (MQ[i].get()) {
1533 // MQ was allocated and computed above; use it
1534 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C2, ONE, *MXj );
1535 }
1536 else if (xc <= qcs[i]) {
1537 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1538 OPT::Apply( *(this->_Op), *Xj, *MXj);
1539 }
1540 }
1541 } // for (int i=0; i<Q.size(); i++)
1542
1543 } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
1544
1545 // Compute the Op-norms after the correction step.
1546 {
1547#ifdef BELOS_TEUCHOS_TIME_MONITOR
1548 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1549#endif
1550 MVT::MvDot( *Xj, *MXj, newDot );
1551 }
1552
1553 // Check for linear dependence.
1554 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1555 dep_flg = true;
1556 }
1557
1558 // Normalize the new vector if it's not dependent
1559 if (!dep_flg) {
1560 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
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 Teuchos::RCP<DM> Bjj = DMT::Subview(*B,1,1,j,j);
1569 DMT::PutScalar(*Bjj, diag);
1570 }
1571 else {
1572 // Create a random vector and orthogonalize it against all previous columns of Q.
1573 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1574 Teuchos::RCP<MV> tempMXj;
1575 MVT::MvRandom( *tempXj );
1576 if (this->_hasOp) {
1577 tempMXj = MVT::Clone( X, 1 );
1578 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1579 }
1580 else {
1581 tempMXj = tempXj;
1582 }
1583 {
1584#ifdef BELOS_TEUCHOS_TIME_MONITOR
1585 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1586#endif
1587 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1588 }
1589 //
1590 for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1591
1592 for (int i=0; i<Q.size(); i++) {
1593 Teuchos::RCP<DM> product = DMT::Create(qcs[i], 1);
1594
1595 // Perform MGS one vector at a time
1596 for (int ii=0; ii<qcs[i]; ii++) {
1597
1598 indX[0] = ii;
1599 tempQ = MVT::CloneView( *Q[i], indX );
1600 Teuchos::RCP tempC = DMT::Subview(*product, 1, 1, ii);
1601
1602 // Apply another step of modified Gram-Schmidt
1603 MatOrthoManager<ScalarType,MV,OP,DM>::innerProd( *tempQ, *tempXj, tempMXj, *tempC );
1604 MVT::MvTimesMatAddMv( -ONE, *tempQ, *tempC, ONE, *tempXj );
1605
1606 }
1607
1608 // Update MXj, with the least number of applications of Op as possible
1609 if (this->_hasOp) {
1610 if (MQ[i].get()) {
1611 // MQ was allocated and computed above; use it
1612 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *product, ONE, *tempMXj );
1613 }
1614 else if (xc <= qcs[i]) {
1615 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1616 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1617 }
1618 }
1619 } // for (int i=0; i<nq; i++)
1620 } // for (int num_orth=0; num_orth<max_orth_steps_; num_orth++)
1621
1622 // Compute the Op-norms after the correction step.
1623 {
1624#ifdef BELOS_TEUCHOS_TIME_MONITOR
1625 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1626#endif
1627 MVT::MvDot( *tempXj, *tempMXj, newDot );
1628 }
1629
1630 // Copy vector into current column of Xj
1631 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1632 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1633
1634 // Enter value on diagonal of B.
1635 Teuchos::RCP<DM> Bjj = DMT::Subview(*B,1,1,j,j);
1636 DMT::PutScalar(*Bjj, ZERO);
1637
1638 // Copy vector into current column of _basisvecs
1639 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1640 if (this->_hasOp) {
1641 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1642 }
1643 }
1644 else {
1645 return j;
1646 }
1647 } // if (!dep_flg)
1648
1649 // Remove the vectors from array
1650 if (j > 0) {
1651 Q.resize( nq );
1652 C.resize( nq );
1653 qcs.resize( nq );
1654 }
1655
1656 } // for (int j=0; j<xc; j++)
1657
1658 return xc;
1659 }
1660
1661 template<class ScalarType, class MV, class OP, class DM>
1662 Teuchos::RCP<Teuchos::ParameterList> getIMGSDefaultParameters ()
1663 {
1664 using Teuchos::ParameterList;
1665 using Teuchos::parameterList;
1666 using Teuchos::RCP;
1667
1668 RCP<ParameterList> params = parameterList ("IMGS");
1669
1670 // Default parameter values for IMGS orthogonalization.
1671 // Documentation will be embedded in the parameter list.
1672 params->set ("maxNumOrthogPasses", IMGSOrthoManager<ScalarType, MV, OP, DM>::max_ortho_steps_default_,
1673 "Maximum number of orthogonalization passes (includes the "
1674 "first). Default is 2, since \"twice is enough\" for Krylov "
1675 "methods.");
1676 params->set ("blkTol", IMGSOrthoManager<ScalarType, MV, OP, DM>::blk_tol_default_,
1677 "Block reorthogonalization threshold.");
1678 params->set ("singTol", IMGSOrthoManager<ScalarType, MV, OP, DM>::sing_tol_default_,
1679 "Singular block detection threshold.");
1680
1681 return params;
1682 }
1683
1684 template<class ScalarType, class MV, class OP, class DM>
1685 Teuchos::RCP<Teuchos::ParameterList> getIMGSFastParameters ()
1686 {
1687 using Teuchos::ParameterList;
1688 using Teuchos::RCP;
1689
1690 RCP<ParameterList> params = getIMGSDefaultParameters<ScalarType, MV, OP>();
1691
1692 params->set ("maxNumOrthogPasses",
1693 IMGSOrthoManager<ScalarType, MV, OP, DM>::max_ortho_steps_fast_);
1694 params->set ("blkTol",
1695 IMGSOrthoManager<ScalarType, MV, OP, DM>::blk_tol_fast_);
1696 params->set ("singTol",
1697 IMGSOrthoManager<ScalarType, MV, OP, DM>::sing_tol_fast_);
1698
1699 return params;
1700 }
1701
1702} // namespace Belos
1703
1704#endif // BELOS_IMGS_ORTHOMANAGER_HPP
1705
BelosConfigDefs.hpp
Belos header file which uses auto-configuration information to include necessary C++ headers.
BelosKokkosDenseAdapter.hpp
Full specialization of Belos::DenseMatTraits for Kokkos::DualView with arbitrary scalarType....
BelosMatOrthoManager.hpp
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
BelosMultiVecTraits.hpp
Declaration of basic traits for the multivector type.
BelosOperatorTraits.hpp
Class which defines basic traits for the operator type.
BelosTeuchosDenseAdapter.hpp
Full specialization of Belos::DenseMatTraits for Teuchos::SerialDenseMatrix with ordinal type int and...
Belos::IMGSOrthoManager
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
Definition BelosIMGSOrthoManager.hpp:53
Belos::IMGSOrthoManager::projectAndNormalizeWithMxImpl
virtual int projectAndNormalizeWithMxImpl(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< DM > > C, Teuchos::RCP< DM > B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for .
Definition BelosIMGSOrthoManager.hpp:571
Belos::IMGSOrthoManager::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Definition BelosIMGSOrthoManager.hpp:205
Belos::IMGSOrthoManager::getSingTol
MagnitudeType getSingTol() const
Return parameter for singular block detection.
Definition BelosIMGSOrthoManager.hpp:229
Belos::IMGSOrthoManager::getBlkTol
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
Definition BelosIMGSOrthoManager.hpp:226
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:416
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:371
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:423
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:386
Belos::IMGSOrthoManager::normalize
int normalize(MV &X, Teuchos::RCP< MV > MX, Teuchos::RCP< DM > B) const
This method takes a multivector X and attempts to compute an orthonormal basis for ,...
Definition BelosIMGSOrthoManager.hpp:749
Belos::IMGSOrthoManager::setParameterList
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
Definition BelosIMGSOrthoManager.hpp:133
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:514
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:97
Belos::IMGSOrthoManager::project
void project(MV &X, Teuchos::Array< Teuchos::RCP< DM > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
This method calls project(X,Teuchos::null,C,Q); see documentation for that function.
Definition BelosIMGSOrthoManager.hpp:271
Belos::IMGSOrthoManager::project
void project(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< DM > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a list of (mutually and internally) orthonormal bases Q, this method takes a multivector X and ...
Definition BelosIMGSOrthoManager.hpp:765
Belos::IMGSOrthoManager::setBlkTol
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
Definition BelosIMGSOrthoManager.hpp:220
Belos::IMGSOrthoManager::setSingTol
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
Definition BelosIMGSOrthoManager.hpp:223
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:408
Belos::IMGSOrthoManager::sing_tol_default_
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
Definition BelosIMGSOrthoManager.hpp:420
Belos::IMGSOrthoManager::normalize
int normalize(MV &X, Teuchos::RCP< DM > B) const
This method calls normalize(X,Teuchos::null,B); see documentation for that function.
Definition BelosIMGSOrthoManager.hpp:309
Belos::IMGSOrthoManager::blk_tol_fast_
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
Definition BelosIMGSOrthoManager.hpp:425
Belos::IMGSOrthoManager::sing_tol_fast_
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
Definition BelosIMGSOrthoManager.hpp:427
Belos::IMGSOrthoManager::blk_tol_default_
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
Definition BelosIMGSOrthoManager.hpp:418
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:67
Belos::MatOrthoManager
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Definition BelosMatOrthoManager.hpp:48
Belos::MatOrthoManager::innerProd
void innerProd(const MV &X, const MV &Y, DM &Z) const
Provides the inner product defining the orthogonality concepts, using the provided operator.
Definition BelosMatOrthoManager.hpp:85
Belos::Operator
Alternative run-time polymorphic interface for operators.
Definition BelosOperator.hpp:48
Belos::Operator::Operator
Operator()
Default constructor (does nothing).
Definition BelosOperator.hpp:54
Belos::getIMGSFastParameters
Teuchos::RCP< Teuchos::ParameterList > getIMGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
Definition BelosIMGSOrthoManager.hpp:1685
Belos::getIMGSDefaultParameters
Teuchos::RCP< Teuchos::ParameterList > getIMGSDefaultParameters()
"Default" parameters for robustness and accuracy.
Definition BelosIMGSOrthoManager.hpp:1662
Generated for Belos by doxygen 1.9.8