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