Belos Version of the Day
Loading...
Searching...
No Matches
BelosDGKSOrthoManager.hpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Belos: Block Linear Solvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
42
47#ifndef BELOS_DGKS_ORTHOMANAGER_HPP
48#define BELOS_DGKS_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> getDGKSDefaultParameters ();
74
76 template<class ScalarType, class MV, class OP>
77 Teuchos::RCP<Teuchos::ParameterList> getDGKSFastParameters();
78
79 template<class ScalarType, class MV, class OP>
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 DGKSOrthoManager( const std::string& label = "Belos",
96 Teuchos::RCP<const OP> Op = Teuchos::null,
97 const int max_blk_ortho = max_blk_ortho_default_,
98 const MagnitudeType blk_tol = blk_tol_default_,
99 const MagnitudeType dep_tol = dep_tol_default_,
100 const MagnitudeType sing_tol = sing_tol_default_ )
101 : MatOrthoManager<ScalarType,MV,OP>(Op),
102 max_blk_ortho_( max_blk_ortho ),
103 blk_tol_( blk_tol ),
104 dep_tol_( dep_tol ),
105 sing_tol_( sing_tol ),
106 label_( label )
107 {
108#ifdef BELOS_TEUCHOS_TIME_MONITOR
109 std::stringstream ss;
110 ss << label_ + ": DGKS[" << max_blk_ortho_ << "]";
111
112 std::string orthoLabel = ss.str() + ": Orthogonalization";
113 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
114
115 std::string updateLabel = ss.str() + ": Ortho (Update)";
116 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
117
118 std::string normLabel = ss.str() + ": Ortho (Norm)";
119 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
120
121 std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
122 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
123#endif
124 }
125
127 DGKSOrthoManager (const Teuchos::RCP<Teuchos::ParameterList>& plist,
128 const std::string& label = "Belos",
129 Teuchos::RCP<const OP> Op = Teuchos::null)
130 : MatOrthoManager<ScalarType,MV,OP>(Op),
131 max_blk_ortho_ ( max_blk_ortho_default_ ),
132 blk_tol_ ( blk_tol_default_ ),
133 dep_tol_ ( dep_tol_default_ ),
134 sing_tol_ ( sing_tol_default_ ),
135 label_( label )
136 {
137 setParameterList (plist);
138
139#ifdef BELOS_TEUCHOS_TIME_MONITOR
140 std::stringstream ss;
141 ss << label_ + ": DGKS[" << max_blk_ortho_ << "]";
142
143 std::string orthoLabel = ss.str() + ": Orthogonalization";
144 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
145
146 std::string updateLabel = ss.str() + ": Ortho (Update)";
147 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
148
149 std::string normLabel = ss.str() + ": Ortho (Norm)";
150 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
151
152 std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
153 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
154#endif
155 }
156
160
162
163
164 void
165 setParameterList (const Teuchos::RCP<Teuchos::ParameterList>& plist)
166 {
167 using Teuchos::ParameterList;
168 using Teuchos::parameterList;
169 using Teuchos::RCP;
170
171 RCP<const ParameterList> defaultParams = getValidParameters();
172 RCP<ParameterList> params;
173 if (plist.is_null()) {
174 // No need to validate default parameters.
175 params = parameterList (*defaultParams);
176 } else {
177 params = plist;
178 params->validateParametersAndSetDefaults (*defaultParams);
179 }
180
181 // Using temporary variables and fetching all values before
182 // setting the output arguments ensures the strong exception
183 // guarantee for this function: if an exception is thrown, no
184 // externally visible side effects (in this case, setting the
185 // output arguments) have taken place.
186 const int maxNumOrthogPasses = params->get<int> ("maxNumOrthogPasses");
187 const MagnitudeType blkTol = params->get<MagnitudeType> ("blkTol");
188 const MagnitudeType depTol = params->get<MagnitudeType> ("depTol");
189 const MagnitudeType singTol = params->get<MagnitudeType> ("singTol");
190
191 max_blk_ortho_ = maxNumOrthogPasses;
192 blk_tol_ = blkTol;
193 dep_tol_ = depTol;
194 sing_tol_ = singTol;
195
196 this->setMyParamList (params);
197 }
198
199 Teuchos::RCP<const Teuchos::ParameterList>
201 {
202 if (defaultParams_.is_null()) {
203 defaultParams_ = Belos::getDGKSDefaultParameters<ScalarType, MV, OP>();
204 }
205
206 return defaultParams_;
207 }
208
210
212
213
215 void setBlkTol( const MagnitudeType blk_tol ) {
216 // Update the parameter list as well.
217 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
218 if (! params.is_null()) {
219 // If it's null, then we haven't called setParameterList()
220 // yet. It's entirely possible to construct the parameter
221 // list on demand, so we don't try to create the parameter
222 // list here.
223 params->set ("blkTol", blk_tol);
224 }
225 blk_tol_ = blk_tol;
226 }
227
229 void setDepTol( const MagnitudeType dep_tol ) {
230 // Update the parameter list as well.
231 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
232 if (! params.is_null()) {
233 params->set ("depTol", dep_tol);
234 }
235 dep_tol_ = dep_tol;
236 }
237
239 void setSingTol( const MagnitudeType sing_tol ) {
240 // Update the parameter list as well.
241 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
242 if (! params.is_null()) {
243 params->set ("singTol", sing_tol);
244 }
245 sing_tol_ = sing_tol;
246 }
247
249 MagnitudeType getBlkTol() const { return blk_tol_; }
250
252 MagnitudeType getDepTol() const { return dep_tol_; }
253
255 MagnitudeType getSingTol() const { return sing_tol_; }
256
258
259
261
262
290 void project ( MV &X, Teuchos::RCP<MV> MX,
291 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
292 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
293
294
297 void project ( MV &X,
298 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
299 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
300 project(X,Teuchos::null,C,Q);
301 }
302
303
304
329 int normalize ( MV &X, Teuchos::RCP<MV> MX,
330 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B) const;
331
332
335 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
336 return normalize(X,Teuchos::null,B);
337 }
338
339 protected:
340
396 virtual int
398 Teuchos::RCP<MV> MX,
399 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
400 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
401 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
402
403 public:
405
407
413 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
414 orthonormError(const MV &X) const {
415 return orthonormError(X,Teuchos::null);
416 }
417
424 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
425 orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const;
426
430 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
431 orthogError(const MV &X1, const MV &X2) const {
432 return orthogError(X1,Teuchos::null,X2);
433 }
434
439 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
440 orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const;
441
443
445
446
449 void setLabel(const std::string& label);
450
453 const std::string& getLabel() const { return label_; }
454
456
458
459
461 static const int max_blk_ortho_default_;
463 static const MagnitudeType blk_tol_default_;
465 static const MagnitudeType dep_tol_default_;
467 static const MagnitudeType sing_tol_default_;
468
470 static const int max_blk_ortho_fast_;
472 static const MagnitudeType blk_tol_fast_;
474 static const MagnitudeType dep_tol_fast_;
476 static const MagnitudeType sing_tol_fast_;
477
479
480 private:
481
483 int max_blk_ortho_;
485 MagnitudeType blk_tol_;
487 MagnitudeType dep_tol_;
489 MagnitudeType sing_tol_;
490
492 std::string label_;
493#ifdef BELOS_TEUCHOS_TIME_MONITOR
494 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
495#endif // BELOS_TEUCHOS_TIME_MONITOR
496
498 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
499
501 int findBasis(MV &X, Teuchos::RCP<MV> MX,
502 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
503 bool completeBasis, int howMany = -1 ) const;
504
506 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
507 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
508 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
509
511 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
512 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
513 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const;
514
528 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
529 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
530 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
531 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const;
532 };
533
534 // Set static variables.
535 template<class ScalarType, class MV, class OP>
537
538 template<class ScalarType, class MV, class OP>
539 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
541 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
542 Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps() );
543
544 template<class ScalarType, class MV, class OP>
545 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
547 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one()
548 / Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
549 2*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one() );
550
551 template<class ScalarType, class MV, class OP>
552 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
554 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
555
556 template<class ScalarType, class MV, class OP>
558
559 template<class ScalarType, class MV, class OP>
560 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
562 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
563
564 template<class ScalarType, class MV, class OP>
565 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
567 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
568
569 template<class ScalarType, class MV, class OP>
570 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
572 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
573
575 // Set the label for this orthogonalization manager and create new timers if it's changed
576 template<class ScalarType, class MV, class OP>
577 void DGKSOrthoManager<ScalarType,MV,OP>::setLabel(const std::string& label)
578 {
579 if (label != label_) {
580 label_ = label;
581#ifdef BELOS_TEUCHOS_TIME_MONITOR
582 std::stringstream ss;
583 ss << label_ + ": DGKS[" << max_blk_ortho_ << "]";
584
585 std::string orthoLabel = ss.str() + ": Orthogonalization";
586 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
587
588 std::string updateLabel = ss.str() + ": Ortho (Update)";
589 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
590
591 std::string normLabel = ss.str() + ": Ortho (Norm)";
592 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
593
594 std::string ipLabel = ss.str() + ": Ortho (Inner Product)";
595 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
596#endif
597 }
598 }
599
601 // Compute the distance from orthonormality
602 template<class ScalarType, class MV, class OP>
603 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
604 DGKSOrthoManager<ScalarType,MV,OP>::orthonormError(const MV &X, Teuchos::RCP<const MV> MX) const {
605 const ScalarType ONE = SCT::one();
606 int rank = MVT::GetNumberVecs(X);
607 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
608 {
609#ifdef BELOS_TEUCHOS_TIME_MONITOR
610 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
611#endif
613 }
614 for (int i=0; i<rank; i++) {
615 xTx(i,i) -= ONE;
616 }
617 return xTx.normFrobenius();
618 }
619
621 // Compute the distance from orthogonality
622 template<class ScalarType, class MV, class OP>
623 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
624 DGKSOrthoManager<ScalarType,MV,OP>::orthogError(const MV &X1, Teuchos::RCP<const MV> MX1, const MV &X2) const {
625 int r1 = MVT::GetNumberVecs(X1);
626 int r2 = MVT::GetNumberVecs(X2);
627 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
628 {
629#ifdef BELOS_TEUCHOS_TIME_MONITOR
630 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
631#endif
633 }
634 return xTx.normFrobenius();
635 }
636
638 // Find an Op-orthonormal basis for span(X) - span(W)
639 template<class ScalarType, class MV, class OP>
640 int
643 Teuchos::RCP<MV> MX,
644 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
645 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
646 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
647 {
648 using Teuchos::Array;
649 using Teuchos::null;
650 using Teuchos::is_null;
651 using Teuchos::RCP;
652 using Teuchos::rcp;
653 using Teuchos::SerialDenseMatrix;
654 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
655 typedef typename Array< RCP< const MV > >::size_type size_type;
656
657#ifdef BELOS_TEUCHOS_TIME_MONITOR
658 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
659#endif
660
661 ScalarType ONE = SCT::one();
662 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
663
664 int nq = Q.size();
665 int xc = MVT::GetNumberVecs( X );
666 ptrdiff_t xr = MVT::GetGlobalLength( X );
667 int rank = xc;
668
669 // If the user doesn't want to store the normalization
670 // coefficients, allocate some local memory for them. This will
671 // go away at the end of this method.
672 if (is_null (B)) {
673 B = rcp (new serial_dense_matrix_type (xc, xc));
674 }
675 // Likewise, if the user doesn't want to store the projection
676 // coefficients, allocate some local memory for them. Also make
677 // sure that all the entries of C are the right size. We're going
678 // to overwrite them anyway, so we don't have to worry about the
679 // contents (other than to resize them if they are the wrong
680 // size).
681 if (C.size() < nq)
682 C.resize (nq);
683 for (size_type k = 0; k < nq; ++k)
684 {
685 const int numRows = MVT::GetNumberVecs (*Q[k]);
686 const int numCols = xc; // Number of vectors in X
687
688 if (is_null (C[k]))
689 C[k] = rcp (new serial_dense_matrix_type (numRows, numCols));
690 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
691 {
692 int err = C[k]->reshape (numRows, numCols);
693 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
694 "DGKS orthogonalization: failed to reshape "
695 "C[" << k << "] (the array of block "
696 "coefficients resulting from projecting X "
697 "against Q[1:" << nq << "]).");
698 }
699 }
700
701 /****** DO NO MODIFY *MX IF _hasOp == false ******/
702 if (this->_hasOp) {
703 if (MX == Teuchos::null) {
704 // we need to allocate space for MX
705 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
706 OPT::Apply(*(this->_Op),X,*MX);
707 }
708 }
709 else {
710 // Op == I --> MX = X (ignore it if the user passed it in)
711 MX = Teuchos::rcp( &X, false );
712 }
713
714 int mxc = MVT::GetNumberVecs( *MX );
715 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
716
717 // short-circuit
718 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
719
720 int numbas = 0;
721 for (int i=0; i<nq; i++) {
722 numbas += MVT::GetNumberVecs( *Q[i] );
723 }
724
725 // check size of B
726 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
727 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
728 // check size of X and MX
729 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
730 "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
731 // check size of X w.r.t. MX
732 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
733 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
734 // check feasibility
735 //TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
736 // "Belos::DGKSOrthoManager::projectAndNormalize(): Orthogonality constraints not feasible" );
737
738 // Some flags for checking dependency returns from the internal orthogonalization methods
739 bool dep_flg = false;
740
741 if (xc == 1) {
742
743 // Use the cheaper block orthogonalization.
744 // NOTE: Don't check for dependencies because the update has one vector.
745 dep_flg = blkOrtho1( X, MX, C, Q );
746
747 // Normalize the new block X
748 if ( B == Teuchos::null ) {
749 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
750 }
751 std::vector<ScalarType> diag(1);
752 {
753#ifdef BELOS_TEUCHOS_TIME_MONITOR
754 Teuchos::TimeMonitor normTimer( *timerNorm_ );
755#endif
756 MVT::MvDot( X, *MX, diag );
757 }
758 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
759
760 if (SCT::magnitude((*B)(0,0)) > ZERO) {
761 rank = 1;
762 MVT::MvScale( X, ONE/(*B)(0,0) );
763 if (this->_hasOp) {
764 // Update MXj.
765 MVT::MvScale( *MX, ONE/(*B)(0,0) );
766 }
767 }
768 }
769 else {
770
771 // Make a temporary copy of X and MX, just in case a block dependency is detected.
772 Teuchos::RCP<MV> tmpX, tmpMX;
773 tmpX = MVT::CloneCopy(X);
774 if (this->_hasOp) {
775 tmpMX = MVT::CloneCopy(*MX);
776 }
777
778 // Use the cheaper block orthogonalization.
779 dep_flg = blkOrtho( X, MX, C, Q );
780
781 // If a dependency has been detected in this block, then perform
782 // the more expensive single-vector orthogonalization.
783 if (dep_flg) {
784 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
785
786 // Copy tmpX back into X.
787 MVT::Assign( *tmpX, X );
788 if (this->_hasOp) {
789 MVT::Assign( *tmpMX, *MX );
790 }
791 }
792 else {
793 // There is no dependency, so orthonormalize new block X
794 rank = findBasis( X, MX, B, false );
795 if (rank < xc) {
796 // A dependency was found during orthonormalization of X,
797 // rerun orthogonalization using more expensive single-vector orthogonalization.
798 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
799
800 // Copy tmpX back into X.
801 MVT::Assign( *tmpX, X );
802 if (this->_hasOp) {
803 MVT::Assign( *tmpMX, *MX );
804 }
805 }
806 }
807 } // if (xc == 1)
808
809 // this should not raise an std::exception; but our post-conditions oblige us to check
810 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
811 "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
812
813 // Return the rank of X.
814 return rank;
815 }
816
817
818
820 // Find an Op-orthonormal basis for span(X), with rank numvectors(X)
821 template<class ScalarType, class MV, class OP>
823 MV &X, Teuchos::RCP<MV> MX,
824 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B ) const {
825
826#ifdef BELOS_TEUCHOS_TIME_MONITOR
827 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
828#endif
829
830 // call findBasis, with the instruction to try to generate a basis of rank numvecs(X)
831 return findBasis(X, MX, B, true);
832
833 }
834
835
836
838 template<class ScalarType, class MV, class OP>
840 MV &X, Teuchos::RCP<MV> MX,
841 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
842 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const {
843 // For the inner product defined by the operator Op or the identity (Op == 0)
844 // -> Orthogonalize X against each Q[i]
845 // Modify MX accordingly
846 //
847 // Note that when Op is 0, MX is not referenced
848 //
849 // Parameter variables
850 //
851 // X : Vectors to be transformed
852 //
853 // MX : Image of the block vector X by the mass matrix
854 //
855 // Q : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently.
856 //
857
858#ifdef BELOS_TEUCHOS_TIME_MONITOR
859 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
860#endif
861
862 int xc = MVT::GetNumberVecs( X );
863 ptrdiff_t xr = MVT::GetGlobalLength( X );
864 int nq = Q.size();
865 std::vector<int> qcs(nq);
866 // short-circuit
867 if (nq == 0 || xc == 0 || xr == 0) {
868 return;
869 }
870 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
871 // if we don't have enough C, expand it with null references
872 // if we have too many, resize to throw away the latter ones
873 // if we have exactly as many as we have Q, this call has no effect
874 C.resize(nq);
875
876
877 /****** DO NO MODIFY *MX IF _hasOp == false ******/
878 if (this->_hasOp) {
879 if (MX == Teuchos::null) {
880 // we need to allocate space for MX
881 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
882 OPT::Apply(*(this->_Op),X,*MX);
883 }
884 }
885 else {
886 // Op == I --> MX = X (ignore it if the user passed it in)
887 MX = Teuchos::rcp( &X, false );
888 }
889 int mxc = MVT::GetNumberVecs( *MX );
890 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
891
892 // check size of X and Q w.r.t. common sense
893 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
894 "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
895 // check size of X w.r.t. MX and Q
896 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
897 "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
898
899 // tally up size of all Q and check/allocate C
900 for (int i=0; i<nq; i++) {
901 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
902 "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
903 qcs[i] = MVT::GetNumberVecs( *Q[i] );
904 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
905 "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
906
907 // check size of C[i]
908 if ( C[i] == Teuchos::null ) {
909 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
910 }
911 else {
912 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
913 "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
914 }
915 }
916
917 // Use the cheaper block orthogonalization, don't check for rank deficiency.
918 blkOrtho( X, MX, C, Q );
919
920 }
921
923 // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that
924 // the rank is numvectors(X)
925 template<class ScalarType, class MV, class OP>
927 MV &X, Teuchos::RCP<MV> MX,
928 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
929 bool completeBasis, int howMany ) const {
930 // For the inner product defined by the operator Op or the identity (Op == 0)
931 // -> Orthonormalize X
932 // Modify MX accordingly
933 //
934 // Note that when Op is 0, MX is not referenced
935 //
936 // Parameter variables
937 //
938 // X : Vectors to be orthonormalized
939 //
940 // MX : Image of the multivector X under the operator Op
941 //
942 // Op : Pointer to the operator for the inner product
943 //
944 //
945
946 const ScalarType ONE = SCT::one();
947 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
948
949 int xc = MVT::GetNumberVecs( X );
950 ptrdiff_t xr = MVT::GetGlobalLength( X );
951
952 if (howMany == -1) {
953 howMany = xc;
954 }
955
956 /*******************************************************
957 * If _hasOp == false, we will not reference MX below *
958 *******************************************************/
959
960 // if Op==null, MX == X (via pointer)
961 // Otherwise, either the user passed in MX or we will allocated and compute it
962 if (this->_hasOp) {
963 if (MX == Teuchos::null) {
964 // we need to allocate space for MX
965 MX = MVT::Clone(X,xc);
966 OPT::Apply(*(this->_Op),X,*MX);
967 }
968 }
969
970 /* if the user doesn't want to store the coefficienets,
971 * allocate some local memory for them
972 */
973 if ( B == Teuchos::null ) {
974 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
975 }
976
977 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
978 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
979
980 // check size of C, B
981 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
982 "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
983 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
984 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
985 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
986 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
987 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
988 "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
989 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
990 "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
991
992 /* xstart is which column we are starting the process with, based on howMany
993 * columns before xstart are assumed to be Op-orthonormal already
994 */
995 int xstart = xc - howMany;
996
997 for (int j = xstart; j < xc; j++) {
998
999 // numX is
1000 // * number of currently orthonormal columns of X
1001 // * the index of the current column of X
1002 int numX = j;
1003 bool addVec = false;
1004
1005 // Get a view of the vector currently being worked on.
1006 std::vector<int> index(1);
1007 index[0] = numX;
1008 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1009 Teuchos::RCP<MV> MXj;
1010 if ((this->_hasOp)) {
1011 // MXj is a view of the current vector in MX
1012 MXj = MVT::CloneViewNonConst( *MX, index );
1013 }
1014 else {
1015 // MXj is a pointer to Xj, and MUST NOT be modified
1016 MXj = Xj;
1017 }
1018
1019 // Get a view of the previous vectors.
1020 std::vector<int> prev_idx( numX );
1021 Teuchos::RCP<const MV> prevX, prevMX;
1022 Teuchos::RCP<MV> oldMXj;
1023
1024 if (numX > 0) {
1025 for (int i=0; i<numX; i++) {
1026 prev_idx[i] = i;
1027 }
1028 prevX = MVT::CloneView( X, prev_idx );
1029 if (this->_hasOp) {
1030 prevMX = MVT::CloneView( *MX, prev_idx );
1031 }
1032
1033 oldMXj = MVT::CloneCopy( *MXj );
1034 }
1035
1036 // Make storage for these Gram-Schmidt iterations.
1037 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1038 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1039 //
1040 // Save old MXj vector and compute Op-norm
1041 //
1042 {
1043#ifdef BELOS_TEUCHOS_TIME_MONITOR
1044 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1045#endif
1046 MVT::MvDot( *Xj, *MXj, oldDot );
1047 }
1048 // Xj^H Op Xj should be real and positive, by the hermitian positive definiteness of Op
1049 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
1050 "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1051
1052 if (numX > 0) {
1053 // Apply the first step of Gram-Schmidt
1054
1055 // product <- prevX^T MXj
1056 {
1057#ifdef BELOS_TEUCHOS_TIME_MONITOR
1058 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1059#endif
1060 MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,product);
1061 }
1062 // Xj <- Xj - prevX prevX^T MXj
1063 // = Xj - prevX product
1064 {
1065#ifdef BELOS_TEUCHOS_TIME_MONITOR
1066 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1067#endif
1068 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1069 }
1070
1071 // Update MXj
1072 if (this->_hasOp) {
1073 // MXj <- Op*Xj_new
1074 // = Op*(Xj_old - prevX prevX^T MXj)
1075 // = MXj - prevMX product
1076#ifdef BELOS_TEUCHOS_TIME_MONITOR
1077 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1078#endif
1079 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1080 }
1081
1082 // Compute new Op-norm
1083 {
1084#ifdef BELOS_TEUCHOS_TIME_MONITOR
1085 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1086#endif
1087 MVT::MvDot( *Xj, *MXj, newDot );
1088 }
1089
1090 // Check if a correction is needed.
1091 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1092 // Apply the second step of Gram-Schmidt
1093 // This is the same as above
1094 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1095 {
1096#ifdef BELOS_TEUCHOS_TIME_MONITOR
1097 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1098#endif
1100 }
1101 product += P2;
1102
1103 {
1104#ifdef BELOS_TEUCHOS_TIME_MONITOR
1105 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1106#endif
1107 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1108 }
1109 if ((this->_hasOp)) {
1110#ifdef BELOS_TEUCHOS_TIME_MONITOR
1111 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1112#endif
1113 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1114 }
1115 } // if (newDot[0] < dep_tol_*oldDot[0])
1116
1117 } // if (numX > 0)
1118
1119 // Compute Op-norm with old MXj
1120 if (numX > 0) {
1121#ifdef BELOS_TEUCHOS_TIME_MONITOR
1122 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1123#endif
1124 MVT::MvDot( *Xj, *oldMXj, newDot );
1125 }
1126 else {
1127 newDot[0] = oldDot[0];
1128 }
1129
1130 // Check to see if the new vector is dependent.
1131 if (completeBasis) {
1132 //
1133 // We need a complete basis, so add random vectors if necessary
1134 //
1135 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1136
1137 // Add a random vector and orthogonalize it against previous vectors in block.
1138 addVec = true;
1139#ifdef ORTHO_DEBUG
1140 std::cout << "Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1141#endif
1142 //
1143 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1144 Teuchos::RCP<MV> tempMXj;
1145 MVT::MvRandom( *tempXj );
1146 if (this->_hasOp) {
1147 tempMXj = MVT::Clone( X, 1 );
1148 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1149 }
1150 else {
1151 tempMXj = tempXj;
1152 }
1153 {
1154#ifdef BELOS_TEUCHOS_TIME_MONITOR
1155 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1156#endif
1157 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1158 }
1159 //
1160 for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
1161 {
1162#ifdef BELOS_TEUCHOS_TIME_MONITOR
1163 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1164#endif
1165 MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,product);
1166 }
1167 {
1168#ifdef BELOS_TEUCHOS_TIME_MONITOR
1169 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1170#endif
1171 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1172 }
1173 if (this->_hasOp) {
1174#ifdef BELOS_TEUCHOS_TIME_MONITOR
1175 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1176#endif
1177 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1178 }
1179 }
1180 // Compute new Op-norm
1181 {
1182#ifdef BELOS_TEUCHOS_TIME_MONITOR
1183 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1184#endif
1185 MVT::MvDot( *tempXj, *tempMXj, newDot );
1186 }
1187 //
1188 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1189 // Copy vector into current column of _basisvecs
1190 MVT::Assign( *tempXj, *Xj );
1191 if (this->_hasOp) {
1192 MVT::Assign( *tempMXj, *MXj );
1193 }
1194 }
1195 else {
1196 return numX;
1197 }
1198 }
1199 }
1200 else {
1201 //
1202 // We only need to detect dependencies.
1203 //
1204 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1205 return numX;
1206 }
1207 }
1208
1209 // If we haven't left this method yet, then we can normalize the new vector Xj.
1210 // Normalize Xj.
1211 // Xj <- Xj / std::sqrt(newDot)
1212 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1213
1214 if (SCT::magnitude(diag) > ZERO) {
1215 MVT::MvScale( *Xj, ONE/diag );
1216 if (this->_hasOp) {
1217 // Update MXj.
1218 MVT::MvScale( *MXj, ONE/diag );
1219 }
1220 }
1221
1222 // If we've added a random vector, enter a zero in the j'th diagonal element.
1223 if (addVec) {
1224 (*B)(j,j) = ZERO;
1225 }
1226 else {
1227 (*B)(j,j) = diag;
1228 }
1229
1230 // Save the coefficients, if we are working on the original vector and not a randomly generated one
1231 if (!addVec) {
1232 for (int i=0; i<numX; i++) {
1233 (*B)(i,j) = product(i,0);
1234 }
1235 }
1236
1237 } // for (j = 0; j < xc; ++j)
1238
1239 return xc;
1240 }
1241
1243 // Routine to compute the block orthogonalization
1244 template<class ScalarType, class MV, class OP>
1245 bool
1246 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1247 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1248 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1249 {
1250 int nq = Q.size();
1251 int xc = MVT::GetNumberVecs( X );
1252 const ScalarType ONE = SCT::one();
1253
1254 std::vector<int> qcs( nq );
1255 for (int i=0; i<nq; i++) {
1256 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1257 }
1258
1259 // Compute the initial Op-norms
1260 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1261 {
1262#ifdef BELOS_TEUCHOS_TIME_MONITOR
1263 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1264#endif
1265 MVT::MvDot( X, *MX, oldDot );
1266 }
1267
1268 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1269 // Define the product Q^T * (Op*X)
1270 for (int i=0; i<nq; i++) {
1271 // Multiply Q' with MX
1272 {
1273#ifdef BELOS_TEUCHOS_TIME_MONITOR
1274 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1275#endif
1277 }
1278 // Multiply by Q and subtract the result in X
1279 {
1280#ifdef BELOS_TEUCHOS_TIME_MONITOR
1281 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1282#endif
1283 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1284 }
1285
1286 // Update MX, with the least number of applications of Op as possible
1287 if (this->_hasOp) {
1288 if (xc <= qcs[i]) {
1289 OPT::Apply( *(this->_Op), X, *MX);
1290 }
1291 else {
1292 // this will possibly be used again below; don't delete it
1293 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1294 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1295 {
1296#ifdef BELOS_TEUCHOS_TIME_MONITOR
1297 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1298#endif
1299 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1300 }
1301 }
1302 }
1303 }
1304
1305 {
1306#ifdef BELOS_TEUCHOS_TIME_MONITOR
1307 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1308#endif
1309 MVT::MvDot( X, *MX, newDot );
1310 }
1311
1312/* // Compute correction bound, compare with PETSc bound.
1313 MagnitudeType hnrm = C[0]->normFrobenius();
1314 for (int i=1; i<nq; i++)
1315 {
1316 hnrm += C[i]->normFrobenius();
1317 }
1318
1319 std::cout << "newDot < 1/sqrt(2)*oldDot < hnrm = " << MGT::squareroot(SCT::magnitude(newDot[0])) << " < " << dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) << " < " << hnrm << std::endl;
1320*/
1321
1322 // Check if a correction is needed.
1323 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1324 // Apply the second step of Gram-Schmidt
1325
1326 for (int i=0; i<nq; i++) {
1327 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1328
1329 // Apply another step of classical Gram-Schmidt
1330 {
1331#ifdef BELOS_TEUCHOS_TIME_MONITOR
1332 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1333#endif
1335 }
1336 *C[i] += C2;
1337
1338 {
1339#ifdef BELOS_TEUCHOS_TIME_MONITOR
1340 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1341#endif
1342 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1343 }
1344
1345 // Update MX, with the least number of applications of Op as possible
1346 if (this->_hasOp) {
1347 if (MQ[i].get()) {
1348#ifdef BELOS_TEUCHOS_TIME_MONITOR
1349 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1350#endif
1351 // MQ was allocated and computed above; use it
1352 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1353 }
1354 else if (xc <= qcs[i]) {
1355 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1356 OPT::Apply( *(this->_Op), X, *MX);
1357 }
1358 }
1359 } // for (int i=0; i<nq; i++)
1360 }
1361
1362 return false;
1363 }
1364
1366 // Routine to compute the block orthogonalization
1367 template<class ScalarType, class MV, class OP>
1368 bool
1369 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1370 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1371 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q) const
1372 {
1373 int nq = Q.size();
1374 int xc = MVT::GetNumberVecs( X );
1375 bool dep_flg = false;
1376 const ScalarType ONE = SCT::one();
1377
1378 std::vector<int> qcs( nq );
1379 for (int i=0; i<nq; i++) {
1380 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1381 }
1382
1383 // Perform the Gram-Schmidt transformation for a block of vectors
1384
1385 // Compute the initial Op-norms
1386 std::vector<ScalarType> oldDot( xc );
1387 {
1388#ifdef BELOS_TEUCHOS_TIME_MONITOR
1389 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1390#endif
1391 MVT::MvDot( X, *MX, oldDot );
1392 }
1393
1394 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1395 // Define the product Q^T * (Op*X)
1396 for (int i=0; i<nq; i++) {
1397 // Multiply Q' with MX
1398 {
1399#ifdef BELOS_TEUCHOS_TIME_MONITOR
1400 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1401#endif
1403 }
1404 // Multiply by Q and subtract the result in X
1405 {
1406#ifdef BELOS_TEUCHOS_TIME_MONITOR
1407 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1408#endif
1409 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1410 }
1411
1412 // Update MX, with the least number of applications of Op as possible
1413 if (this->_hasOp) {
1414 if (xc <= qcs[i]) {
1415 OPT::Apply( *(this->_Op), X, *MX);
1416 }
1417 else {
1418 // this will possibly be used again below; don't delete it
1419 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1420 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1421 {
1422#ifdef BELOS_TEUCHOS_TIME_MONITOR
1423 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1424#endif
1425 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1426 }
1427 }
1428 }
1429 }
1430
1431 // Do as many steps of classical Gram-Schmidt as required by max_blk_ortho_
1432 for (int j = 1; j < max_blk_ortho_; ++j) {
1433
1434 for (int i=0; i<nq; i++) {
1435 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1436
1437 // Apply another step of classical Gram-Schmidt
1438 {
1439#ifdef BELOS_TEUCHOS_TIME_MONITOR
1440 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1441#endif
1443 }
1444 *C[i] += C2;
1445
1446 {
1447#ifdef BELOS_TEUCHOS_TIME_MONITOR
1448 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1449#endif
1450 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1451 }
1452
1453 // Update MX, with the least number of applications of Op as possible
1454 if (this->_hasOp) {
1455 if (MQ[i].get()) {
1456#ifdef BELOS_TEUCHOS_TIME_MONITOR
1457 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1458#endif
1459 // MQ was allocated and computed above; use it
1460 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1461 }
1462 else if (xc <= qcs[i]) {
1463 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1464 OPT::Apply( *(this->_Op), X, *MX);
1465 }
1466 }
1467 } // for (int i=0; i<nq; i++)
1468 } // for (int j = 0; j < max_blk_ortho; ++j)
1469
1470 // Compute new Op-norms
1471 std::vector<ScalarType> newDot(xc);
1472 {
1473#ifdef BELOS_TEUCHOS_TIME_MONITOR
1474 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1475#endif
1476 MVT::MvDot( X, *MX, newDot );
1477 }
1478
1479 // Check to make sure the new block of vectors are not dependent on previous vectors
1480 for (int i=0; i<xc; i++){
1481 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1482 dep_flg = true;
1483 break;
1484 }
1485 } // end for (i=0;...)
1486
1487 return dep_flg;
1488 }
1489
1490
1491 template<class ScalarType, class MV, class OP>
1492 int
1493 DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1494 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1495 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1496 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ) const
1497 {
1498 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1499
1500 const ScalarType ONE = SCT::one();
1501 const ScalarType ZERO = SCT::zero();
1502
1503 int nq = Q.size();
1504 int xc = MVT::GetNumberVecs( X );
1505 std::vector<int> indX( 1 );
1506 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1507
1508 std::vector<int> qcs( nq );
1509 for (int i=0; i<nq; i++) {
1510 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1511 }
1512
1513 // Create pointers for the previous vectors of X that have already been orthonormalized.
1514 Teuchos::RCP<const MV> lastQ;
1515 Teuchos::RCP<MV> Xj, MXj;
1516 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1517
1518 // Perform the Gram-Schmidt transformation for each vector in the block of vectors.
1519 for (int j=0; j<xc; j++) {
1520
1521 bool dep_flg = false;
1522
1523 // Get a view of the previously orthogonalized vectors and B, add it to the arrays.
1524 if (j > 0) {
1525 std::vector<int> index( j );
1526 for (int ind=0; ind<j; ind++) {
1527 index[ind] = ind;
1528 }
1529 lastQ = MVT::CloneView( X, index );
1530
1531 // Add these views to the Q and C arrays.
1532 Q.push_back( lastQ );
1533 C.push_back( B );
1534 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1535 }
1536
1537 // Get a view of the current vector in X to orthogonalize.
1538 indX[0] = j;
1539 Xj = MVT::CloneViewNonConst( X, indX );
1540 if (this->_hasOp) {
1541 MXj = MVT::CloneViewNonConst( *MX, indX );
1542 }
1543 else {
1544 MXj = Xj;
1545 }
1546
1547 // Compute the initial Op-norms
1548 {
1549#ifdef BELOS_TEUCHOS_TIME_MONITOR
1550 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1551#endif
1552 MVT::MvDot( *Xj, *MXj, oldDot );
1553 }
1554
1555 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1556 // Define the product Q^T * (Op*X)
1557 for (int i=0; i<Q.size(); i++) {
1558
1559 // Get a view of the current serial dense matrix
1560 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1561
1562 // Multiply Q' with MX
1563 {
1564#ifdef BELOS_TEUCHOS_TIME_MONITOR
1565 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1566#endif
1568 }
1569 // Multiply by Q and subtract the result in Xj
1570 {
1571#ifdef BELOS_TEUCHOS_TIME_MONITOR
1572 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1573#endif
1574 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1575 }
1576
1577 // Update MXj, with the least number of applications of Op as possible
1578 if (this->_hasOp) {
1579 if (xc <= qcs[i]) {
1580 OPT::Apply( *(this->_Op), *Xj, *MXj);
1581 }
1582 else {
1583 // this will possibly be used again below; don't delete it
1584 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1585 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1586 {
1587#ifdef BELOS_TEUCHOS_TIME_MONITOR
1588 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1589#endif
1590 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1591 }
1592 }
1593 }
1594 }
1595
1596 // Compute the Op-norms
1597 {
1598#ifdef BELOS_TEUCHOS_TIME_MONITOR
1599 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1600#endif
1601 MVT::MvDot( *Xj, *MXj, newDot );
1602 }
1603
1604 // Do one step of classical Gram-Schmidt orthogonalization
1605 // with a second correction step if needed.
1606
1607 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
1608
1609 for (int i=0; i<Q.size(); i++) {
1610 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1611 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1612
1613 // Apply another step of classical Gram-Schmidt
1614 {
1615#ifdef BELOS_TEUCHOS_TIME_MONITOR
1616 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1617#endif
1619 }
1620 tempC += C2;
1621 {
1622#ifdef BELOS_TEUCHOS_TIME_MONITOR
1623 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1624#endif
1625 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1626 }
1627
1628 // Update MXj, with the least number of applications of Op as possible
1629 if (this->_hasOp) {
1630 if (MQ[i].get()) {
1631#ifdef BELOS_TEUCHOS_TIME_MONITOR
1632 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1633#endif
1634 // MQ was allocated and computed above; use it
1635 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1636 }
1637 else if (xc <= qcs[i]) {
1638 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1639 OPT::Apply( *(this->_Op), *Xj, *MXj);
1640 }
1641 }
1642 } // for (int i=0; i<Q.size(); i++)
1643
1644 // Compute the Op-norms after the correction step.
1645 {
1646#ifdef BELOS_TEUCHOS_TIME_MONITOR
1647 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1648#endif
1649 MVT::MvDot( *Xj, *MXj, newDot );
1650 }
1651 } // if ()
1652
1653 // Check for linear dependence.
1654 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1655 dep_flg = true;
1656 }
1657
1658 // Normalize the new vector if it's not dependent
1659 if (!dep_flg) {
1660 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1661
1662 MVT::MvScale( *Xj, ONE/diag );
1663 if (this->_hasOp) {
1664 // Update MXj.
1665 MVT::MvScale( *MXj, ONE/diag );
1666 }
1667
1668 // Enter value on diagonal of B.
1669 (*B)(j,j) = diag;
1670 }
1671 else {
1672 // Create a random vector and orthogonalize it against all previous columns of Q.
1673 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1674 Teuchos::RCP<MV> tempMXj;
1675 MVT::MvRandom( *tempXj );
1676 if (this->_hasOp) {
1677 tempMXj = MVT::Clone( X, 1 );
1678 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1679 }
1680 else {
1681 tempMXj = tempXj;
1682 }
1683 {
1684#ifdef BELOS_TEUCHOS_TIME_MONITOR
1685 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1686#endif
1687 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1688 }
1689 //
1690 for (int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
1691
1692 for (int i=0; i<Q.size(); i++) {
1693 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1694
1695 // Apply another step of classical Gram-Schmidt
1696 {
1697#ifdef BELOS_TEUCHOS_TIME_MONITOR
1698 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1699#endif
1700 MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*tempXj,tempMXj,product);
1701 }
1702 {
1703#ifdef BELOS_TEUCHOS_TIME_MONITOR
1704 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1705#endif
1706 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1707 }
1708
1709 // Update MXj, with the least number of applications of Op as possible
1710 if (this->_hasOp) {
1711 if (MQ[i].get()) {
1712#ifdef BELOS_TEUCHOS_TIME_MONITOR
1713 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1714#endif
1715 // MQ was allocated and computed above; use it
1716 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1717 }
1718 else if (xc <= qcs[i]) {
1719 // MQ was not allocated and computed above; it was cheaper to use X before and it still is
1720 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1721 }
1722 }
1723 } // for (int i=0; i<nq; i++)
1724
1725 }
1726
1727 // Compute the Op-norms after the correction step.
1728 {
1729#ifdef BELOS_TEUCHOS_TIME_MONITOR
1730 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1731#endif
1732 MVT::MvDot( *tempXj, *tempMXj, newDot );
1733 }
1734
1735 // Copy vector into current column of Xj
1736 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1737 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1738
1739 // Enter value on diagonal of B.
1740 (*B)(j,j) = ZERO;
1741
1742 // Copy vector into current column of _basisvecs
1743 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1744 if (this->_hasOp) {
1745 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1746 }
1747 }
1748 else {
1749 return j;
1750 }
1751 } // if (!dep_flg)
1752
1753 // Remove the vectors from array
1754 if (j > 0) {
1755 Q.resize( nq );
1756 C.resize( nq );
1757 qcs.resize( nq );
1758 }
1759
1760 } // for (int j=0; j<xc; j++)
1761
1762 return xc;
1763 }
1764
1765 template<class ScalarType, class MV, class OP>
1766 Teuchos::RCP<Teuchos::ParameterList> getDGKSDefaultParameters ()
1767 {
1768 using Teuchos::ParameterList;
1769 using Teuchos::parameterList;
1770 using Teuchos::RCP;
1771
1772 RCP<ParameterList> params = parameterList ("DGKS");
1773
1774 // Default parameter values for DGKS orthogonalization.
1775 // Documentation will be embedded in the parameter list.
1776 params->set ("maxNumOrthogPasses", DGKSOrthoManager<ScalarType, MV, OP>::max_blk_ortho_default_,
1777 "Maximum number of orthogonalization passes (includes the "
1778 "first). Default is 2, since \"twice is enough\" for Krylov "
1779 "methods.");
1781 "Block reorthogonalization threshold.");
1783 "(Non-block) reorthogonalization threshold.");
1785 "Singular block detection threshold.");
1786
1787 return params;
1788 }
1789
1790 template<class ScalarType, class MV, class OP>
1791 Teuchos::RCP<Teuchos::ParameterList> getDGKSFastParameters ()
1792 {
1793 using Teuchos::ParameterList;
1794 using Teuchos::RCP;
1795
1796 RCP<ParameterList> params = getDGKSDefaultParameters<ScalarType, MV, OP>();
1797
1798 params->set ("maxNumOrthogPasses",
1800 params->set ("blkTol",
1802 params->set ("depTol",
1804 params->set ("singTol",
1806
1807 return params;
1808 }
1809
1810} // namespace Belos
1811
1812#endif // BELOS_DGKS_ORTHOMANAGER_HPP
1813
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 (potentially) m...
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
DGKSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_blk_ortho=max_blk_ortho_default_, const MagnitudeType blk_tol=blk_tol_default_, const MagnitudeType dep_tol=dep_tol_default_, const MagnitudeType sing_tol=sing_tol_default_)
Constructor specifying re-orthogonalization tolerance.
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_blk_ortho_default_
Max number of (re)orthogonalization steps, including the first (default).
MagnitudeType getSingTol() const
Return parameter for singular block detection.
DGKSOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &plist, const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor that takes a list of parameters.
static const MagnitudeType dep_tol_default_
(Non-block) reorthogonalization threshold (default).
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.
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
static const int max_blk_ortho_fast_
Max number of (re)orthogonalization steps, including the first (fast).
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
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.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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 .
MagnitudeType getDepTol() const
Return parameter for re-orthogonalization threshhold.
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
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...
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
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 dep_tol_fast_
(Non-block) reorthogonalization threshold (fast).
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 setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
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 > getDGKSDefaultParameters()
"Default" parameters for robustness and accuracy.
Teuchos::RCP< Teuchos::ParameterList > getDGKSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.

Generated for Belos by doxygen 1.9.6