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