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