Belos Version of the Day
Loading...
Searching...
No Matches
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
46
47#ifndef BELOS_IMGS_ORTHOMANAGER_HPP
48#define BELOS_IMGS_ORTHOMANAGER_HPP
49
56
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 for (int i=0; i<nq; i++) {
859 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
860 "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
861 qcs[i] = MVT::GetNumberVecs( *Q[i] );
862 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
863 "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
864
865 // check size of C[i]
866 if ( C[i] == Teuchos::null ) {
867 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
868 }
869 else {
870 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
871 "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
872 }
873 }
874
875 // Use the cheaper block orthogonalization, don't check for rank deficiency.
876 blkOrtho( X, MX, C, Q );
877
878 }
879
881 // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
882 // the rank is numvectors(X)
883 template<class ScalarType, class MV, class OP>
884 int IMGSOrthoManager<ScalarType, MV, OP>::findBasis(
885 MV &X, Teuchos::RCP<MV> MX,
886 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
887 bool completeBasis, int howMany ) const {
888 // For the inner product defined by the operator Op or the identity (Op == 0)
889 // -> Orthonormalize X
890 // Modify MX accordingly
891 //
892 // Note that when Op is 0, MX is not referenced
893 //
894 // Parameter variables
895 //
896 // X : Vectors to be orthonormalized
897 //
898 // MX : Image of the multivector X under the operator Op
899 //
900 // Op : Pointer to the operator for the inner product
901 //
902 //
903
904 const ScalarType ONE = SCT::one();
905 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
906
907 int xc = MVT::GetNumberVecs( X );
908 ptrdiff_t xr = MVT::GetGlobalLength( X );
909
910 if (howMany == -1) {
911 howMany = xc;
912 }
913
914 /*******************************************************
915 * If _hasOp == false, we will not reference MX below *
916 *******************************************************/
917
918 // if Op==null, MX == X (via pointer)
919 // Otherwise, either the user passed in MX or we will allocated and compute it
920 if (this->_hasOp) {
921 if (MX == Teuchos::null) {
922 // we need to allocate space for MX
923 MX = MVT::Clone(X,xc);
924 OPT::Apply(*(this->_Op),X,*MX);
925 }
926 }
927
928 /* if the user doesn't want to store the coefficienets,
929 * allocate some local memory for them
930 */
931 if ( B == Teuchos::null ) {
932 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
933 }
934
935 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
936 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
937
938 // check size of C, B
939 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
940 "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
941 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
942 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
943 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
944 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
945 TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
946 "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
947 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
948 "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
949
950 /* xstart is which column we are starting the process with, based on howMany
951 * columns before xstart are assumed to be Op-orthonormal already
952 */
953 int xstart = xc - howMany;
954
955 for (int j = xstart; j < xc; j++) {
956
957 // numX is
958 // * number of currently orthonormal columns of X
959 // * the index of the current column of X
960 int numX = j;
961 bool addVec = false;
962
963 // Get a view of the vector currently being worked on.
964 std::vector<int> index(1);
965 index[0] = numX;
966 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
967 Teuchos::RCP<MV> MXj;
968 if ((this->_hasOp)) {
969 // MXj is a view of the current vector in MX
970 MXj = MVT::CloneViewNonConst( *MX, index );
971 }
972 else {
973 // MXj is a pointer to Xj, and MUST NOT be modified
974 MXj = Xj;
975 }
976
977 Teuchos::RCP<MV> oldMXj;
978 if (numX > 0) {
979 oldMXj = MVT::CloneCopy( *MXj );
980 }
981
982 // Make storage for these Gram-Schmidt iterations.
983 Teuchos::SerialDenseVector<int,ScalarType> product(numX);
984 Teuchos::SerialDenseVector<int,ScalarType> P2(1);
985 Teuchos::RCP<const MV> prevX, prevMX;
986
987 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
988 //
989 // Save old MXj vector and compute Op-norm
990 //
991 {
992#ifdef BELOS_TEUCHOS_TIME_MONITOR
993 Teuchos::TimeMonitor normTimer( *timerNorm_ );
994#endif
995 MVT::MvDot( *Xj, *MXj, oldDot );
996 }
997 // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
998 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
999 "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1000
1001 // Perform MGS one vector at a time
1002 for (int ii=0; ii<numX; ii++) {
1003
1004 index[0] = ii;
1005 prevX = MVT::CloneView( X, index );
1006 if (this->_hasOp) {
1007 prevMX = MVT::CloneView( *MX, index );
1008 }
1009
1010 for (int i=0; i<max_ortho_steps_; ++i) {
1011
1012 // product <- prevX^T MXj
1013 {
1014#ifdef BELOS_TEUCHOS_TIME_MONITOR
1015 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1016#endif
1017 MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
1018 }
1019
1020 // Xj <- Xj - prevX prevX^T MXj
1021 // = Xj - prevX product
1022 {
1023#ifdef BELOS_TEUCHOS_TIME_MONITOR
1024 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1025#endif
1026 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1027 }
1028
1029 // Update MXj
1030 if (this->_hasOp) {
1031 // MXj <- Op*Xj_new
1032 // = Op*(Xj_old - prevX prevX^T MXj)
1033 // = MXj - prevMX product
1034#ifdef BELOS_TEUCHOS_TIME_MONITOR
1035 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1036#endif
1037 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1038 }
1039
1040 // Set coefficients
1041 if ( i==0 )
1042 product[ii] = P2[0];
1043 else
1044 product[ii] += P2[0];
1045
1046 } // for (int i=0; i<max_ortho_steps_; ++i)
1047
1048 } // for (int ii=0; ii<numX; ++ii)
1049
1050 // Compute Op-norm with old MXj
1051 if (numX > 0) {
1052#ifdef BELOS_TEUCHOS_TIME_MONITOR
1053 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1054#endif
1055 MVT::MvDot( *Xj, *oldMXj, newDot );
1056 }
1057 else {
1058 newDot[0] = oldDot[0];
1059 }
1060
1061 // Check to see if the new vector is dependent.
1062 if (completeBasis) {
1063 //
1064 // We need a complete basis, so add random vectors if necessary
1065 //
1066 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1067
1068 // Add a random vector and orthogonalize it against previous vectors in block.
1069 addVec = true;
1070#ifdef ORTHO_DEBUG
1071 std::cout << "Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1072#endif
1073 //
1074 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1075 Teuchos::RCP<MV> tempMXj;
1076 MVT::MvRandom( *tempXj );
1077 if (this->_hasOp) {
1078 tempMXj = MVT::Clone( X, 1 );
1079 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1080 }
1081 else {
1082 tempMXj = tempXj;
1083 }
1084 {
1085#ifdef BELOS_TEUCHOS_TIME_MONITOR
1086 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1087#endif
1088 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1089 }
1090 //
1091 // Perform MGS one vector at a time
1092 for (int ii=0; ii<numX; ii++) {
1093
1094 index[0] = ii;
1095 prevX = MVT::CloneView( X, index );
1096 if (this->_hasOp) {
1097 prevMX = MVT::CloneView( *MX, index );
1098 }
1099
1100 for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1101 {
1102#ifdef BELOS_TEUCHOS_TIME_MONITOR
1103 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1104#endif
1105 MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,P2);
1106 }
1107 {
1108#ifdef BELOS_TEUCHOS_TIME_MONITOR
1109 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1110#endif
1111 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
1112 }
1113 if (this->_hasOp) {
1114#ifdef BELOS_TEUCHOS_TIME_MONITOR
1115 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1116#endif
1117 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
1118 }
1119
1120 // Set coefficients
1121 if ( num_orth==0 )
1122 product[ii] = P2[0];
1123 else
1124 product[ii] += P2[0];
1125 }
1126 }
1127
1128 // Compute new Op-norm
1129 {
1130#ifdef BELOS_TEUCHOS_TIME_MONITOR
1131 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1132#endif
1133 MVT::MvDot( *tempXj, *tempMXj, newDot );
1134 }
1135 //
1136 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1137 // Copy vector into current column of _basisvecs
1138 MVT::Assign( *tempXj, *Xj );
1139 if (this->_hasOp) {
1140 MVT::Assign( *tempMXj, *MXj );
1141 }
1142 }
1143 else {
1144 return numX;
1145 }
1146 }
1147 }
1148 else {
1149 //
1150 // We only need to detect dependencies.
1151 //
1152 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1153 return numX;
1154 }
1155 }
1156
1157
1158 // If we haven't left this method yet, then we can normalize the new vector Xj.
1159 // Normalize Xj.
1160 // Xj <- Xj / std::sqrt(newDot)
1161 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1162 if (SCT::magnitude(diag) > ZERO) {
1163 MVT::MvScale( *Xj, ONE/diag );
1164 if (this->_hasOp) {
1165 // Update MXj.
1166 MVT::MvScale( *MXj, ONE/diag );
1167 }
1168 }
1169
1170 // If we've added a random vector, enter a zero in the j'th diagonal element.
1171 if (addVec) {
1172 (*B)(j,j) = ZERO;
1173 }
1174 else {
1175 (*B)(j,j) = diag;
1176 }
1177
1178 // Save the coefficients, if we are working on the original vector and not a randomly generated one
1179 if (!addVec) {
1180 for (int i=0; i<numX; i++) {
1181 (*B)(i,j) = product(i);
1182 }
1183 }
1184
1185 } // for (j = 0; j < xc; ++j)
1186
1187 return xc;
1188 }
1189
1191 // Routine to compute the block orthogonalization
1192 template<class ScalarType, class MV, class OP>
1193 bool
1194 IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1195 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1196 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1197 {
1198 int nq = Q.size();
1199 int xc = MVT::GetNumberVecs( X );
1200 const ScalarType ONE = SCT::one();
1201
1202 std::vector<int> qcs( nq );
1203 for (int i=0; i<nq; i++) {
1204 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1205 }
1206
1207 // Perform the Gram-Schmidt transformation for a block of vectors
1208 std::vector<int> index(1);
1209 Teuchos::RCP<const MV> tempQ;
1210
1211 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1212 // Define the product Q^T * (Op*X)
1213 for (int i=0; i<nq; i++) {
1214
1215 // Perform MGS one vector at a time
1216 for (int ii=0; ii<qcs[i]; ii++) {
1217
1218 index[0] = ii;
1219 tempQ = MVT::CloneView( *Q[i], index );
1220 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1221
1222 // Multiply Q' with MX
1223 {
1224#ifdef BELOS_TEUCHOS_TIME_MONITOR
1225 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1226#endif
1227 MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,X,MX,tempC);
1228 }
1229 // Multiply by Q and subtract the result in X
1230 {
1231#ifdef BELOS_TEUCHOS_TIME_MONITOR
1232 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1233#endif
1234 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1235 }
1236 }
1237
1238 // Update MX, with the least number of applications of Op as possible
1239 if (this->_hasOp) {
1240 if (xc <= qcs[i]) {
1241 OPT::Apply( *(this->_Op), X, *MX);
1242 }
1243 else {
1244 // this will possibly be used again below; don't delete it
1245 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1246 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1247 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1248 }
1249 }
1250 }
1251
1252 // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1253 for (int j = 1; j < max_ortho_steps_; ++j) {
1254
1255 for (int i=0; i<nq; i++) {
1256
1257 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1258
1259 // Perform MGS one vector at a time
1260 for (int ii=0; ii<qcs[i]; ii++) {
1261
1262 index[0] = ii;
1263 tempQ = MVT::CloneView( *Q[i], index );
1264 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1265 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1266
1267 // Apply another step of modified Gram-Schmidt
1268 {
1269#ifdef BELOS_TEUCHOS_TIME_MONITOR
1270 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1271#endif
1272 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1273 }
1274 tempC += tempC2;
1275 {
1276#ifdef BELOS_TEUCHOS_TIME_MONITOR
1277 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1278#endif
1279 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1280 }
1281
1282 }
1283
1284 // Update MX, with the least number of applications of Op as possible
1285 if (this->_hasOp) {
1286 if (MQ[i].get()) {
1287 // MQ was allocated and computed above; use it
1288 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1289 }
1290 else if (xc <= qcs[i]) {
1291 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1292 OPT::Apply( *(this->_Op), X, *MX);
1293 }
1294 }
1295 } // for (int i=0; i<nq; i++)
1296 } // for (int j = 0; j < max_ortho_steps; ++j)
1297
1298 return false;
1299 }
1300
1302 // Routine to compute the block orthogonalization
1303 template<class ScalarType, class MV, class OP>
1304 bool
1305 IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1306 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1307 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1308 {
1309 int nq = Q.size();
1310 int xc = MVT::GetNumberVecs( X );
1311 bool dep_flg = false;
1312 const ScalarType ONE = SCT::one();
1313
1314 std::vector<int> qcs( nq );
1315 for (int i=0; i<nq; i++) {
1316 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1317 }
1318
1319 // Perform the Gram-Schmidt transformation for a block of vectors
1320
1321 // Compute the initial Op-norms
1322 std::vector<ScalarType> oldDot( xc );
1323 {
1324#ifdef BELOS_TEUCHOS_TIME_MONITOR
1325 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1326#endif
1327 MVT::MvDot( X, *MX, oldDot );
1328 }
1329
1330 std::vector<int> index(1);
1331 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1332 Teuchos::RCP<const MV> tempQ;
1333
1334 // Define the product Q^T * (Op*X)
1335 for (int i=0; i<nq; i++) {
1336
1337 // Perform MGS one vector at a time
1338 for (int ii=0; ii<qcs[i]; ii++) {
1339
1340 index[0] = ii;
1341 tempQ = MVT::CloneView( *Q[i], index );
1342 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1343
1344 // Multiply Q' with MX
1345 {
1346#ifdef BELOS_TEUCHOS_TIME_MONITOR
1347 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1348#endif
1349 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC);
1350 }
1351 // Multiply by Q and subtract the result in X
1352 {
1353#ifdef BELOS_TEUCHOS_TIME_MONITOR
1354 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1355#endif
1356 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1357 }
1358 }
1359
1360 // Update MX, with the least number of applications of Op as possible
1361 if (this->_hasOp) {
1362 if (xc <= qcs[i]) {
1363 OPT::Apply( *(this->_Op), X, *MX);
1364 }
1365 else {
1366 // this will possibly be used again below; don't delete it
1367 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1368 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1369 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1370 }
1371 }
1372 }
1373
1374 // Do as many steps of modified Gram-Schmidt as required by max_ortho_steps_
1375 for (int j = 1; j < max_ortho_steps_; ++j) {
1376
1377 for (int i=0; i<nq; i++) {
1378 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1379
1380 // Perform MGS one vector at a time
1381 for (int ii=0; ii<qcs[i]; ii++) {
1382
1383 index[0] = ii;
1384 tempQ = MVT::CloneView( *Q[i], index );
1385 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1386 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1387
1388 // Apply another step of modified Gram-Schmidt
1389 {
1390#ifdef BELOS_TEUCHOS_TIME_MONITOR
1391 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1392#endif
1393 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1394 }
1395 tempC += tempC2;
1396 {
1397#ifdef BELOS_TEUCHOS_TIME_MONITOR
1398 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1399#endif
1400 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1401 }
1402 }
1403
1404 // Update MX, with the least number of applications of Op as possible
1405 if (this->_hasOp) {
1406 if (MQ[i].get()) {
1407 // MQ was allocated and computed above; use it
1408 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1409 }
1410 else if (xc <= qcs[i]) {
1411 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1412 OPT::Apply( *(this->_Op), X, *MX);
1413 }
1414 }
1415 } // for (int i=0; i<nq; i++)
1416 } // for (int j = 0; j < max_ortho_steps; ++j)
1417
1418 // Compute new Op-norms
1419 std::vector<ScalarType> newDot(xc);
1420 {
1421#ifdef BELOS_TEUCHOS_TIME_MONITOR
1422 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1423#endif
1424 MVT::MvDot( X, *MX, newDot );
1425 }
1426
1427 // Check to make sure the new block of vectors are not dependent on previous vectors
1428 for (int i=0; i<xc; i++){
1429 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1430 dep_flg = true;
1431 break;
1432 }
1433 } // end for (i=0;...)
1434
1435 return dep_flg;
1436 }
1437
1438 template<class ScalarType, class MV, class OP>
1439 int
1440 IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1441 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1442 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1443 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
1444 {
1445 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1446
1447 const ScalarType ONE = SCT::one();
1448 const ScalarType ZERO = SCT::zero();
1449
1450 int nq = Q.size();
1451 int xc = MVT::GetNumberVecs( X );
1452 std::vector<int> indX( 1 );
1453 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1454
1455 std::vector<int> qcs( nq );
1456 for (int i=0; i<nq; i++) {
1457 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1458 }
1459
1460 // Create pointers for the previous vectors of X that have already been orthonormalized.
1461 Teuchos::RCP<const MV> lastQ;
1462 Teuchos::RCP<MV> Xj, MXj;
1463 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1464
1465 // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
1466 for (int j=0; j<xc; j++) {
1467
1468 bool dep_flg = false;
1469
1470 // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
1471 if (j > 0) {
1472 std::vector<int> index( j );
1473 for (int ind=0; ind<j; ind++) {
1474 index[ind] = ind;
1475 }
1476 lastQ = MVT::CloneView( X, index );
1477
1478 // Add these views to the Q and C arrays.
1479 Q.push_back( lastQ );
1480 C.push_back( B );
1481 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1482 }
1483
1484 // Get a view of the current vector in X to orthogonalize.
1485 indX[0] = j;
1486 Xj = MVT::CloneViewNonConst( X, indX );
1487 if (this->_hasOp) {
1488 MXj = MVT::CloneViewNonConst( *MX, indX );
1489 }
1490 else {
1491 MXj = Xj;
1492 }
1493
1494 // Compute the initial Op-norms
1495 {
1496#ifdef BELOS_TEUCHOS_TIME_MONITOR
1497 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1498#endif
1499 MVT::MvDot( *Xj, *MXj, oldDot );
1500 }
1501
1502 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1503 Teuchos::RCP<const MV> tempQ;
1504
1505 // Define the product Q^T * (Op*X)
1506 for (int i=0; i<Q.size(); i++) {
1507
1508 // Perform MGS one vector at a time
1509 for (int ii=0; ii<qcs[i]; ii++) {
1510
1511 indX[0] = ii;
1512 tempQ = MVT::CloneView( *Q[i], indX );
1513 // Get a view of the current serial dense matrix
1514 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1515
1516 // Multiply Q' with MX
1517 MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,*Xj,MXj,tempC);
1518
1519 // Multiply by Q and subtract the result in Xj
1520 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1521 }
1522
1523 // Update MXj, with the least number of applications of Op as possible
1524 if (this->_hasOp) {
1525 if (xc <= qcs[i]) {
1526 OPT::Apply( *(this->_Op), *Xj, *MXj);
1527 }
1528 else {
1529 // this will possibly be used again below; don't delete it
1530 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1531 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1532 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1533 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1534 }
1535 }
1536 }
1537
1538 // Do any additional steps of modified Gram-Schmidt orthogonalization
1539 for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1540
1541 for (int i=0; i<Q.size(); i++) {
1542 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1543
1544 // Perform MGS one vector at a time
1545 for (int ii=0; ii<qcs[i]; ii++) {
1546
1547 indX[0] = ii;
1548 tempQ = MVT::CloneView( *Q[i], indX );
1549 // Get a view of the current serial dense matrix
1550 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1551
1552 // Apply another step of modified Gram-Schmidt
1553 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *Xj, MXj, tempC2);
1554 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1555 }
1556
1557 // Add the coefficients into C[i]
1558 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1559 tempC += C2;
1560
1561 // Update MXj, with the least number of applications of Op as possible
1562 if (this->_hasOp) {
1563 if (MQ[i].get()) {
1564 // MQ was allocated and computed above; use it
1565 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1566 }
1567 else if (xc <= qcs[i]) {
1568 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1569 OPT::Apply( *(this->_Op), *Xj, *MXj);
1570 }
1571 }
1572 } // for (int i=0; i<Q.size(); i++)
1573
1574 } // for (int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps)
1575
1576 // Compute the Op-norms after the correction step.
1577 {
1578#ifdef BELOS_TEUCHOS_TIME_MONITOR
1579 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1580#endif
1581 MVT::MvDot( *Xj, *MXj, newDot );
1582 }
1583
1584 // Check for linear dependence.
1585 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1586 dep_flg = true;
1587 }
1588
1589 // Normalize the new vector if it's not dependent
1590 if (!dep_flg) {
1591 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1592
1593 MVT::MvScale( *Xj, ONE/diag );
1594 if (this->_hasOp) {
1595 // Update MXj.
1596 MVT::MvScale( *MXj, ONE/diag );
1597 }
1598
1599 // Enter value on diagonal of B.
1600 (*B)(j,j) = diag;
1601 }
1602 else {
1603 // Create a random vector and orthogonalize it against all previous columns of Q.
1604 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1605 Teuchos::RCP<MV> tempMXj;
1606 MVT::MvRandom( *tempXj );
1607 if (this->_hasOp) {
1608 tempMXj = MVT::Clone( X, 1 );
1609 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1610 }
1611 else {
1612 tempMXj = tempXj;
1613 }
1614 {
1615#ifdef BELOS_TEUCHOS_TIME_MONITOR
1616 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1617#endif
1618 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1619 }
1620 //
1621 for (int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1622
1623 for (int i=0; i<Q.size(); i++) {
1624 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1625
1626 // Perform MGS one vector at a time
1627 for (int ii=0; ii<qcs[i]; ii++) {
1628
1629 indX[0] = ii;
1630 tempQ = MVT::CloneView( *Q[i], indX );
1631 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1632
1633 // Apply another step of modified Gram-Schmidt
1634 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *tempXj, tempMXj, tempC );
1635 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1636
1637 }
1638
1639 // Update MXj, with the least number of applications of Op as possible
1640 if (this->_hasOp) {
1641 if (MQ[i].get()) {
1642 // MQ was allocated and computed above; use it
1643 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1644 }
1645 else if (xc <= qcs[i]) {
1646 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1647 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1648 }
1649 }
1650 } // for (int i=0; i<nq; i++)
1651 } // for (int num_orth=0; num_orth<max_orth_steps_; num_orth++)
1652
1653 // Compute the Op-norms after the correction step.
1654 {
1655#ifdef BELOS_TEUCHOS_TIME_MONITOR
1656 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1657#endif
1658 MVT::MvDot( *tempXj, *tempMXj, newDot );
1659 }
1660
1661 // Copy vector into current column of Xj
1662 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1663 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1664
1665 // Enter value on diagonal of B.
1666 (*B)(j,j) = ZERO;
1667
1668 // Copy vector into current column of _basisvecs
1669 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1670 if (this->_hasOp) {
1671 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1672 }
1673 }
1674 else {
1675 return j;
1676 }
1677 } // if (!dep_flg)
1678
1679 // Remove the vectors from array
1680 if (j > 0) {
1681 Q.resize( nq );
1682 C.resize( nq );
1683 qcs.resize( nq );
1684 }
1685
1686 } // for (int j=0; j<xc; j++)
1687
1688 return xc;
1689 }
1690
1691 template<class ScalarType, class MV, class OP>
1692 Teuchos::RCP<Teuchos::ParameterList> getIMGSDefaultParameters ()
1693 {
1694 using Teuchos::ParameterList;
1695 using Teuchos::parameterList;
1696 using Teuchos::RCP;
1697
1698 RCP<ParameterList> params = parameterList ("IMGS");
1699
1700 // Default parameter values for IMGS orthogonalization.
1701 // Documentation will be embedded in the parameter list.
1702 params->set ("maxNumOrthogPasses", IMGSOrthoManager<ScalarType, MV, OP>::max_ortho_steps_default_,
1703 "Maximum number of orthogonalization passes (includes the "
1704 "first). Default is 2, since \"twice is enough\" for Krylov "
1705 "methods.");
1706 params->set ("blkTol", IMGSOrthoManager<ScalarType, MV, OP>::blk_tol_default_,
1707 "Block reorthogonalization threshold.");
1708 params->set ("singTol", IMGSOrthoManager<ScalarType, MV, OP>::sing_tol_default_,
1709 "Singular block detection threshold.");
1710
1711 return params;
1712 }
1713
1714 template<class ScalarType, class MV, class OP>
1715 Teuchos::RCP<Teuchos::ParameterList> getIMGSFastParameters ()
1716 {
1717 using Teuchos::ParameterList;
1718 using Teuchos::RCP;
1719
1720 RCP<ParameterList> params = getIMGSDefaultParameters<ScalarType, MV, OP>();
1721
1722 params->set ("maxNumOrthogPasses",
1723 IMGSOrthoManager<ScalarType, MV, OP>::max_ortho_steps_fast_);
1724 params->set ("blkTol",
1725 IMGSOrthoManager<ScalarType, MV, OP>::blk_tol_fast_);
1726 params->set ("singTol",
1727 IMGSOrthoManager<ScalarType, MV, OP>::sing_tol_fast_);
1728
1729 return params;
1730 }
1731
1732} // namespace Belos
1733
1734#endif // BELOS_IMGS_ORTHOMANAGER_HPP
1735
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::setParameterList
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
Definition BelosIMGSOrthoManager.hpp:163
Belos::IMGSOrthoManager::sing_tol_fast_
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
Definition BelosIMGSOrthoManager.hpp:457
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::sing_tol_default_
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
Definition BelosIMGSOrthoManager.hpp:450
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::blk_tol_fast_
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
Definition BelosIMGSOrthoManager.hpp:455
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::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::blk_tol_default_
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
Definition BelosIMGSOrthoManager.hpp:448
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::MatOrthoManager
MatOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null)
Default constructor.
Definition BelosMatOrthoManager.hpp:87
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::MatOrthoManager::_Op
Teuchos::RCP< const OP > _Op
Definition BelosMatOrthoManager.hpp:80
Belos::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition BelosMultiVecTraits.hpp:129
Belos::MultiVecTraits::MvScale
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Definition BelosMultiVecTraits.hpp:265
Belos::MultiVecTraits::CloneCopy
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Definition BelosMultiVecTraits.hpp:145
Belos::MultiVecTraits::GetGlobalLength
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Definition BelosMultiVecTraits.hpp:214
Belos::MultiVecTraits::MvDot
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
Definition BelosMultiVecTraits.hpp:280
Belos::MultiVecTraits::Clone
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
Definition BelosMultiVecTraits.hpp:138
Belos::MultiVecTraits::GetNumberVecs
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
Definition BelosMultiVecTraits.hpp:218
Belos::MultiVecTraits::Assign
static void Assign(const MV &A, MV &mv)
mv := A
Definition BelosMultiVecTraits.hpp:327
Belos::OperatorTraits
Class which defines basic traits for the operator type.
Definition BelosOperatorTraits.hpp:110
Belos::OperatorTraits::Apply
static void Apply(const OP &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Apply Op to x, putting the result into y.
Definition BelosOperatorTraits.hpp:123
Belos::OrthoError
Exception thrown to signal error in an orthogonalization manager method.
Definition BelosOrthoManager.hpp:84
Belos::getIMGSFastParameters
Teuchos::RCP< Teuchos::ParameterList > getIMGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
Definition BelosIMGSOrthoManager.hpp:1715
Belos::getIMGSDefaultParameters
Teuchos::RCP< Teuchos::ParameterList > getIMGSDefaultParameters()
"Default" parameters for robustness and accuracy.
Definition BelosIMGSOrthoManager.hpp:1692
Generated for Belos by doxygen 1.13.2