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