51#ifndef BELOS_ICGS_ORTHOMANAGER_MP_VECTOR_HPP
52#define BELOS_ICGS_ORTHOMANAGER_MP_VECTOR_HPP
55#include "BelosICGSOrthoManager.hpp"
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>
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;
77 Teuchos::RCP<const OP> Op = Teuchos::null,
78 const int max_ortho_steps = max_ortho_steps_default_,
82 max_ortho_steps_( max_ortho_steps ),
84 sing_tol_( sing_tol ),
87#ifdef BELOS_TEUCHOS_TIME_MONITOR
89 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
91 std::string orthoLabel = ss.str() +
": Orthogonalization";
92 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
94 std::string updateLabel = ss.str() +
": Ortho (Update)";
95 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
97 std::string normLabel = ss.str() +
": Ortho (Norm)";
98 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
100 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
101 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
107 const std::string& label =
"Belos",
108 Teuchos::RCP<const OP> Op = Teuchos::null) :
110 max_ortho_steps_ (max_ortho_steps_default_),
111 blk_tol_ (blk_tol_default_),
112 sing_tol_ (sing_tol_default_),
115 setParameterList (plist);
117#ifdef BELOS_TEUCHOS_TIME_MONITOR
118 std::stringstream ss;
119 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
121 std::string orthoLabel = ss.str() +
": Orthogonalization";
122 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
124 std::string updateLabel = ss.str() +
": Ortho (Update)";
125 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
127 std::string normLabel = ss.str() +
": Ortho (Norm)";
128 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
130 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
131 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
145 using Teuchos::Exceptions::InvalidParameterName;
146 using Teuchos::ParameterList;
147 using Teuchos::parameterList;
150 RCP<const ParameterList> defaultParams = getValidParameters();
151 RCP<ParameterList> params;
152 if (plist.is_null()) {
153 params = parameterList (*defaultParams);
168 int maxNumOrthogPasses;
173 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
174 }
catch (InvalidParameterName&) {
175 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
176 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
188 }
catch (InvalidParameterName&) {
193 params->remove (
"depTol");
194 }
catch (InvalidParameterName&) {
197 params->set (
"blkTol", blkTol);
202 }
catch (InvalidParameterName&) {
204 params->set (
"singTol", singTol);
207 max_ortho_steps_ = maxNumOrthogPasses;
211 this->setMyParamList (params);
214 Teuchos::RCP<const Teuchos::ParameterList>
217 if (defaultParams_.is_null()) {
218 defaultParams_ = Belos::getICGSDefaultParameters<ScalarType, MV, OP>();
221 return defaultParams_;
230 Teuchos::RCP<const Teuchos::ParameterList>
234 using Teuchos::ParameterList;
235 using Teuchos::parameterList;
238 RCP<const ParameterList> defaultParams = getValidParameters ();
240 RCP<ParameterList> params = parameterList (*defaultParams);
242 params->set (
"maxNumOrthogPasses", max_ortho_steps_fast_);
243 params->set (
"blkTol", blk_tol_fast_);
244 params->set (
"singTol", sing_tol_fast_);
255 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
256 if (! params.is_null()) {
261 params->set (
"blkTol", blk_tol);
269 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
270 if (! params.is_null()) {
275 params->set (
"singTol", sing_tol);
277 sing_tol_ = sing_tol;
320 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
321 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
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);
359 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
364 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
365 return normalize(X,Teuchos::null,B);
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;
426 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
428 return orthonormError(X,Teuchos::null);
435 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
441 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
443 return orthogError(X1,Teuchos::null,X2);
450 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
451 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
464 const std::string&
getLabel()
const {
return label_; }
498#ifdef BELOS_TEUCHOS_TIME_MONITOR
499 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
507 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
508 bool completeBasis,
int howMany = -1 )
const;
512 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
513 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
517 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
518 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
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;
540 template<
class StorageType,
class MV,
class OP>
541 const int ICGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::max_ortho_steps_default_ = 2;
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(
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();
554 template<
class StorageType,
class MV,
class OP>
555 const int ICGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::max_ortho_steps_fast_ = 1;
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();
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();
569 template<
class StorageType,
class MV,
class OP>
570 void ICGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::setLabel(
const std::string& label)
572 if (label != label_) {
574#ifdef BELOS_TEUCHOS_TIME_MONITOR
575 std::stringstream ss;
576 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
578 std::string orthoLabel = ss.str() +
": Orthogonalization";
579 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
581 std::string updateLabel = ss.str() +
": Ortho (Update)";
582 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
584 std::string normLabel = ss.str() +
": Ortho (Norm)";
585 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
587 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
588 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
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++) {
605 return xTx.normFrobenius();
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();
622 template<
class StorageType,
class MV,
class OP>
624 ICGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::
625 projectAndNormalizeWithMxImpl (MV &X,
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
631 using Teuchos::Array;
633 using Teuchos::is_null;
636 using Teuchos::SerialDenseMatrix;
637 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
638 typedef typename Array< RCP< const MV > >::size_type size_type;
640#ifdef BELOS_TEUCHOS_TIME_MONITOR
641 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
644 ScalarType ONE = SCT::one();
645 const MagnitudeType ZERO = MGT::zero();
648 int xc = MVT::GetNumberVecs( X );
649 ptrdiff_t xr = MVT::GetGlobalLength( X );
656 B = rcp (
new serial_dense_matrix_type (xc, xc));
666 for (size_type k = 0; k < nq; ++k)
668 const int numRows = MVT::GetNumberVecs (*Q[k]);
669 const int numCols = xc;
672 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
673 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
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 <<
"]).");
686 if (MX == Teuchos::null) {
688 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
689 OPT::Apply(*(this->_Op),X,*MX);
694 MX = Teuchos::rcp( &X,
false );
697 int mxc = MVT::GetNumberVecs( *MX );
698 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
701 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::ICGSOrthoManager::projectAndNormalize(): X must be non-empty" );
704 for (
int i=0; i<nq; i++) {
705 numbas += MVT::GetNumberVecs( *Q[i] );
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" );
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" );
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" );
722 bool dep_flg =
false;
728 dep_flg = blkOrtho1( X, MX, C, Q );
731 if ( B == Teuchos::null ) {
732 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
734 std::vector<ScalarType>
diag(xc);
736#ifdef BELOS_TEUCHOS_TIME_MONITOR
737 Teuchos::TimeMonitor normTimer( *timerNorm_ );
739 MVT::MvDot( X, *MX,
diag );
741 (*B)(0,0) = SCT::squareroot(SCT::magnitude(
diag[0]));
743 ScalarType scale = ONE;
748 MVT::MvScale( X, scale );
751 MVT::MvScale( *MX, scale );
757 Teuchos::RCP<MV> tmpX, tmpMX;
758 tmpX = MVT::CloneCopy(X);
760 tmpMX = MVT::CloneCopy(*MX);
764 dep_flg = blkOrtho( X, MX, C, Q );
770 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
773 MVT::Assign( *tmpX, X );
775 MVT::Assign( *tmpMX, *MX );
780 rank = findBasis( X, MX, B,
false );
785 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
788 MVT::Assign( *tmpX, X );
790 MVT::Assign( *tmpMX, *MX );
797 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
798 "Belos::ICGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
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 {
813#ifdef BELOS_TEUCHOS_TIME_MONITOR
814 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
818 return findBasis(X, MX, B,
true);
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 {
845#ifdef BELOS_TEUCHOS_TIME_MONITOR
846 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
849 int xc = MVT::GetNumberVecs( X );
850 ptrdiff_t xr = MVT::GetGlobalLength( X );
852 std::vector<int> qcs(nq);
854 if (nq == 0 || xc == 0 || xr == 0) {
857 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
866 if (MX == Teuchos::null) {
868 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
869 OPT::Apply(*(this->_Op),X,*MX);
874 MX = Teuchos::rcp( &X,
false );
876 int mxc = MVT::GetNumberVecs( *MX );
877 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
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" );
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" );
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" );
897 if ( C[i] == Teuchos::null ) {
898 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
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" );
907 blkOrtho( X, MX, C, Q );
914 template<
class StorageType,
class MV,
class OP>
916 ICGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::
919 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
937 const ScalarType ONE = SCT::one ();
938 const MagnitudeType ZERO = SCT::magnitude (SCT::zero ());
940 const int xc = MVT::GetNumberVecs (X);
941 const ptrdiff_t xr = MVT::GetGlobalLength (X);
954 if (MX == Teuchos::null) {
956 MX = MVT::Clone(X,xc);
957 OPT::Apply(*(this->_Op),X,*MX);
964 if ( B == Teuchos::null ) {
965 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
968 const int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
969 const ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
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" );
986 int xstart = xc - howMany;
988 for (
int j = xstart;
j < xc;
j++) {
997 std::vector<int> index(1);
999 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1000 Teuchos::RCP<MV> MXj;
1003 MXj = MVT::CloneViewNonConst( *MX, index );
1011 std::vector<int> prev_idx( numX );
1012 Teuchos::RCP<const MV> prevX, prevMX;
1013 Teuchos::RCP<MV> oldMXj;
1016 for (
int i=0; i<numX; i++) {
1019 prevX = MVT::CloneView( X, prev_idx );
1021 prevMX = MVT::CloneView( *MX, prev_idx );
1024 oldMXj = MVT::CloneCopy( *MXj );
1028 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1029 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1034#ifdef BELOS_TEUCHOS_TIME_MONITOR
1035 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1037 MVT::MvDot( *Xj, *MXj, oldDot );
1040 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
1041 "Belos::ICGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1045 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1047 for (
int i=0; i<max_ortho_steps_; ++i) {
1051#ifdef BELOS_TEUCHOS_TIME_MONITOR
1052 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1054 MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
1060#ifdef BELOS_TEUCHOS_TIME_MONITOR
1061 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1063 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1071#ifdef BELOS_TEUCHOS_TIME_MONITOR
1072 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1074 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1088#ifdef BELOS_TEUCHOS_TIME_MONITOR
1089 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1091 MVT::MvDot( *Xj, *oldMXj, newDot );
1094 newDot[0] = oldDot[0];
1098 if (completeBasis) {
1102 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1107 std::cout <<
"Belos::ICGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1110 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1111 Teuchos::RCP<MV> tempMXj;
1112 MVT::MvRandom( *tempXj );
1114 tempMXj = MVT::Clone( X, 1 );
1115 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1121#ifdef BELOS_TEUCHOS_TIME_MONITOR
1122 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1124 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1127 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1129#ifdef BELOS_TEUCHOS_TIME_MONITOR
1130 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1132 MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,product);
1135#ifdef BELOS_TEUCHOS_TIME_MONITOR
1136 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1138 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1141#ifdef BELOS_TEUCHOS_TIME_MONITOR
1142 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1144 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1149#ifdef BELOS_TEUCHOS_TIME_MONITOR
1150 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1152 MVT::MvDot( *tempXj, *tempMXj, newDot );
1155 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1157 MVT::Assign( *tempXj, *Xj );
1159 MVT::Assign( *tempMXj, *MXj );
1171 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1179 ScalarType
diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1180 ScalarType scale = ONE;
1182 MVT::MvScale( *Xj, scale );
1185 MVT::MvScale( *MXj, scale );
1198 for (
int i=0; i<numX; i++) {
1199 (*B)(i,
j) = product(i,0);
1210 template<
class StorageType,
class MV,
class OP>
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
1217 int xc = MVT::GetNumberVecs( X );
1218 const ScalarType ONE = SCT::one();
1220 std::vector<int> qcs( nq );
1221 for (
int i=0; i<nq; i++) {
1222 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1227 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1229 for (
int i=0; i<nq; i++) {
1232#ifdef BELOS_TEUCHOS_TIME_MONITOR
1233 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1235 MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
1239#ifdef BELOS_TEUCHOS_TIME_MONITOR
1240 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1242 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1248 OPT::Apply( *(this->_Op), X, *MX);
1252 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1253 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1255#ifdef BELOS_TEUCHOS_TIME_MONITOR
1256 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1258 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1265 for (
int j = 1;
j < max_ortho_steps_; ++
j) {
1267 for (
int i=0; i<nq; i++) {
1268 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1272#ifdef BELOS_TEUCHOS_TIME_MONITOR
1273 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1275 MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
1279#ifdef BELOS_TEUCHOS_TIME_MONITOR
1280 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1282 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1288#ifdef BELOS_TEUCHOS_TIME_MONITOR
1289 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1292 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1294 else if (xc <= qcs[i]) {
1296 OPT::Apply( *(this->_Op), X, *MX);
1307 template<
class StorageType,
class MV,
class OP>
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
1314 int xc = MVT::GetNumberVecs( X );
1315 bool dep_flg =
false;
1316 const ScalarType ONE = SCT::one();
1318 std::vector<int> qcs( nq );
1319 for (
int i=0; i<nq; i++) {
1320 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1326 std::vector<ScalarType> oldDot( xc );
1328#ifdef BELOS_TEUCHOS_TIME_MONITOR
1329 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1331 MVT::MvDot( X, *MX, oldDot );
1334 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1336 for (
int i=0; i<nq; i++) {
1339#ifdef BELOS_TEUCHOS_TIME_MONITOR
1340 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1342 MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,*C[i]);
1346#ifdef BELOS_TEUCHOS_TIME_MONITOR
1347 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1349 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1354 OPT::Apply( *(this->_Op), X, *MX);
1358 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1359 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1361#ifdef BELOS_TEUCHOS_TIME_MONITOR
1362 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1364 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1371 for (
int j = 1;
j < max_ortho_steps_; ++
j) {
1373 for (
int i=0; i<nq; i++) {
1374 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1378#ifdef BELOS_TEUCHOS_TIME_MONITOR
1379 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1381 MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],X,MX,C2);
1385#ifdef BELOS_TEUCHOS_TIME_MONITOR
1386 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1388 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1394#ifdef BELOS_TEUCHOS_TIME_MONITOR
1395 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1398 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1400 else if (xc <= qcs[i]) {
1402 OPT::Apply( *(this->_Op), X, *MX);
1409 std::vector<ScalarType> newDot(xc);
1411#ifdef BELOS_TEUCHOS_TIME_MONITOR
1412 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1414 MVT::MvDot( X, *MX, newDot );
1418 for (
int i=0; i<xc; i++){
1419 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1428 template<
class StorageType,
class MV,
class OP>
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
1435 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1437 const ScalarType ONE = SCT::one();
1438 const ScalarType ZERO = SCT::zero();
1441 int xc = MVT::GetNumberVecs( X );
1442 std::vector<int> indX( 1 );
1443 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1445 std::vector<int> qcs( nq );
1446 for (
int i=0; i<nq; i++) {
1447 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1451 Teuchos::RCP<const MV> lastQ;
1452 Teuchos::RCP<MV> Xj, MXj;
1453 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1456 for (
int j=0;
j<xc;
j++) {
1458 bool dep_flg =
false;
1462 std::vector<int> index(
j );
1463 for (
int ind=0; ind<
j; ind++) {
1466 lastQ = MVT::CloneView( X, index );
1469 Q.push_back( lastQ );
1471 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1476 Xj = MVT::CloneViewNonConst( X, indX );
1478 MXj = MVT::CloneViewNonConst( *MX, indX );
1486#ifdef BELOS_TEUCHOS_TIME_MONITOR
1487 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1489 MVT::MvDot( *Xj, *MXj, oldDot );
1492 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1494 for (
int i=0; i<Q.size(); i++) {
1497 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0,
j );
1501#ifdef BELOS_TEUCHOS_TIME_MONITOR
1502 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1504 MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,tempC);
1507#ifdef BELOS_TEUCHOS_TIME_MONITOR
1508 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1511 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1516 OPT::Apply( *(this->_Op), *Xj, *MXj);
1520 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1521 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1523#ifdef BELOS_TEUCHOS_TIME_MONITOR
1524 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1526 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1533 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
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 );
1541#ifdef BELOS_TEUCHOS_TIME_MONITOR
1542 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1544 MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*Xj,MXj,C2);
1548#ifdef BELOS_TEUCHOS_TIME_MONITOR
1549 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1551 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1558#ifdef BELOS_TEUCHOS_TIME_MONITOR
1559 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1561 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1563 else if (xc <= qcs[i]) {
1565 OPT::Apply( *(this->_Op), *Xj, *MXj);
1573 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1579 ScalarType
diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1581 MVT::MvScale( *Xj, ONE/
diag );
1584 MVT::MvScale( *MXj, ONE/
diag );
1592 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1593 Teuchos::RCP<MV> tempMXj;
1594 MVT::MvRandom( *tempXj );
1596 tempMXj = MVT::Clone( X, 1 );
1597 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1603#ifdef BELOS_TEUCHOS_TIME_MONITOR
1604 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1606 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1609 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1611 for (
int i=0; i<Q.size(); i++) {
1612 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1616#ifdef BELOS_TEUCHOS_TIME_MONITOR
1617 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1619 MatOrthoManager<ScalarType,MV,OP>::innerProd(*Q[i],*tempXj,tempMXj,product);
1622#ifdef BELOS_TEUCHOS_TIME_MONITOR
1623 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1625 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1631#ifdef BELOS_TEUCHOS_TIME_MONITOR
1632 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1635 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1637 else if (xc <= qcs[i]) {
1639 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1648#ifdef BELOS_TEUCHOS_TIME_MONITOR
1649 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1651 MVT::MvDot( *tempXj, *tempMXj, newDot );
1655 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1656 ScalarType
diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1662 MVT::MvAddMv( ONE/
diag, *tempXj, ZERO, *tempXj, *Xj );
1664 MVT::MvAddMv( ONE/
diag, *tempMXj, ZERO, *tempMXj, *MXj );
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 MagnitudeType
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.
std::string label_
Label for timers.
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
Teuchos::ScalarTraits< ScalarType > SCT
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.
Sacado::MP::Vector< Storage > ScalarType
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.
~ICGSOrthoManager()
Destructor.
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::ScalarTraits< MagnitudeType > MGT
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
MultiVecTraits< ScalarType, MV > MVT
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
MagnitudeType blk_tol_
Block reorthogonalization threshold.
OperatorTraits< ScalarType, MV, OP > OPT
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 ...
MagnitudeType sing_tol_
Singular block detection threshold.
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)