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

Generated for Belos by doxygen 1.10.0