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