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