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