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