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
47#ifndef BELOS_ICGS_ORTHOMANAGER_HPP
48#define BELOS_ICGS_ORTHOMANAGER_HPP
49
57// #define ORTHO_DEBUG
58
59#include "BelosConfigDefs.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>
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;
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
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>
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>
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
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>
561
562 template<class ScalarType, class MV, class OP>
563 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
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
571 = 10*Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
572
573 template<class ScalarType, class MV, class OP>
575
576 template<class ScalarType, class MV, class OP>
577 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
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
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);
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);
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
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>
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>
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
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
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
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
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
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
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
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
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.");
1723 "Block reorthogonalization threshold.");
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",
1740 params->set ("blkTol",
1742 params->set ("singTol",
1744
1745 return params;
1746 }
1747
1748} // namespace Belos
1749
1750#endif // BELOS_ICGS_ORTHOMANAGER_HPP
1751
Belos header file which uses auto-configuration information to include necessary C++ headers.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
MagnitudeType getSingTol() const
Return parameter for singular block detection.
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
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.
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
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.
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 ...
static const int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
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.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
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 ...
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 ,...
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
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.
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...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
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 .
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Teuchos::RCP< Teuchos::ParameterList > getICGSDefaultParameters()
"Default" parameters for robustness and accuracy.
Teuchos::RCP< Teuchos::ParameterList > getICGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.

Generated for Belos by doxygen 1.10.0