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