47 #ifndef BELOS_ICGS_ORTHOMANAGER_HPP 48 #define BELOS_ICGS_ORTHOMANAGER_HPP 64 #include "Teuchos_as.hpp" 65 #ifdef BELOS_TEUCHOS_TIME_MONITOR 66 #include "Teuchos_TimeMonitor.hpp" 67 #endif // BELOS_TEUCHOS_TIME_MONITOR 72 template<
class ScalarType,
class MV,
class OP>
76 template<
class ScalarType,
class MV,
class OP>
79 template<
class ScalarType,
class MV,
class OP>
84 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
85 typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
86 typedef Teuchos::ScalarTraits<ScalarType> SCT;
96 Teuchos::RCP<const OP> Op = Teuchos::null,
101 max_ortho_steps_( max_ortho_steps ),
103 sing_tol_( sing_tol ),
106 #ifdef BELOS_TEUCHOS_TIME_MONITOR 107 std::stringstream ss;
108 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
110 std::string orthoLabel = ss.str() +
": Orthogonalization";
111 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
113 std::string updateLabel = ss.str() +
": Ortho (Update)";
114 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
116 std::string normLabel = ss.str() +
": Ortho (Norm)";
117 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
119 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
120 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
126 const std::string& label =
"Belos",
127 Teuchos::RCP<const OP> Op = Teuchos::null) :
136 #ifdef BELOS_TEUCHOS_TIME_MONITOR 137 std::stringstream ss;
138 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
140 std::string orthoLabel = ss.str() +
": Orthogonalization";
141 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
143 std::string updateLabel = ss.str() +
": Ortho (Update)";
144 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
146 std::string normLabel = ss.str() +
": Ortho (Norm)";
147 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
149 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
150 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
164 using Teuchos::Exceptions::InvalidParameterName;
165 using Teuchos::ParameterList;
166 using Teuchos::parameterList;
170 RCP<ParameterList> params;
171 if (plist.is_null()) {
172 params = parameterList (*defaultParams);
187 int maxNumOrthogPasses;
188 MagnitudeType blkTol;
189 MagnitudeType singTol;
192 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
193 }
catch (InvalidParameterName&) {
194 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
195 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
206 blkTol = params->get<MagnitudeType> (
"blkTol");
207 }
catch (InvalidParameterName&) {
209 blkTol = params->get<MagnitudeType> (
"depTol");
212 params->remove (
"depTol");
213 }
catch (InvalidParameterName&) {
214 blkTol = defaultParams->get<MagnitudeType> (
"blkTol");
216 params->set (
"blkTol", blkTol);
220 singTol = params->get<MagnitudeType> (
"singTol");
221 }
catch (InvalidParameterName&) {
222 singTol = defaultParams->get<MagnitudeType> (
"singTol");
223 params->set (
"singTol", singTol);
226 max_ortho_steps_ = maxNumOrthogPasses;
230 this->setMyParamList (params);
233 Teuchos::RCP<const Teuchos::ParameterList>
236 if (defaultParams_.is_null()) {
237 defaultParams_ = Belos::getICGSDefaultParameters<ScalarType, MV, OP>();
240 return defaultParams_;
249 Teuchos::RCP<const Teuchos::ParameterList>
253 using Teuchos::ParameterList;
254 using Teuchos::parameterList;
259 RCP<ParameterList> params = parameterList (*defaultParams);
274 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
275 if (! params.is_null()) {
280 params->set (
"blkTol", blk_tol);
288 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
289 if (! params.is_null()) {
294 params->set (
"singTol", sing_tol);
296 sing_tol_ = sing_tol;
338 void project ( MV &X, Teuchos::RCP<MV> MX,
339 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
340 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
346 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
347 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
377 int normalize ( MV &X, Teuchos::RCP<MV> MX,
378 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
383 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
433 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
434 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
435 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
445 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
454 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
460 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
469 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
470 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
479 void setLabel(
const std::string& label);
483 const std::string&
getLabel()
const {
return label_; }
509 int max_ortho_steps_;
511 MagnitudeType blk_tol_;
513 MagnitudeType sing_tol_;
517 #ifdef BELOS_TEUCHOS_TIME_MONITOR 518 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
519 #endif // BELOS_TEUCHOS_TIME_MONITOR 522 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
525 int findBasis(MV &X, Teuchos::RCP<MV> MX,
526 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
527 bool completeBasis,
int howMany = -1 )
const;
530 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
531 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
532 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
535 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
536 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
537 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
552 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
553 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
554 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
555 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
559 template<
class ScalarType,
class MV,
class OP>
562 template<
class ScalarType,
class MV,
class OP>
563 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
565 = 10*Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
566 Teuchos::ScalarTraits<
typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps() );
568 template<
class ScalarType,
class MV,
class OP>
569 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
571 = 10*Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
573 template<
class ScalarType,
class MV,
class OP>
576 template<
class ScalarType,
class MV,
class OP>
577 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
579 = Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
581 template<
class ScalarType,
class MV,
class OP>
582 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
584 = Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
588 template<
class ScalarType,
class MV,
class OP>
591 if (label != label_) {
593 #ifdef BELOS_TEUCHOS_TIME_MONITOR 594 std::stringstream ss;
595 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
597 std::string orthoLabel = ss.str() +
": Orthogonalization";
598 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
600 std::string updateLabel = ss.str() +
": Ortho (Update)";
601 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
603 std::string normLabel = ss.str() +
": Ortho (Norm)";
604 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
606 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
607 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
614 template<
class ScalarType,
class MV,
class OP>
615 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
617 const ScalarType ONE = SCT::one();
618 int rank = MVT::GetNumberVecs(X);
619 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
621 for (
int i=0; i<rank; i++) {
624 return xTx.normFrobenius();
629 template<
class ScalarType,
class MV,
class OP>
630 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
632 int r1 = MVT::GetNumberVecs(X1);
633 int r2 = MVT::GetNumberVecs(X2);
634 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
636 return xTx.normFrobenius();
641 template<
class ScalarType,
class MV,
class OP>
646 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
647 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
648 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 650 using Teuchos::Array;
652 using Teuchos::is_null;
655 using Teuchos::SerialDenseMatrix;
656 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
657 typedef typename Array< RCP< const MV > >::size_type size_type;
659 #ifdef BELOS_TEUCHOS_TIME_MONITOR 660 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
663 ScalarType ONE = SCT::one();
664 const MagnitudeType ZERO = MGT::zero();
667 int xc = MVT::GetNumberVecs( X );
668 ptrdiff_t xr = MVT::GetGlobalLength( X );
675 B = rcp (
new serial_dense_matrix_type (xc, xc));
685 for (size_type k = 0; k < nq; ++k)
687 const int numRows = MVT::GetNumberVecs (*Q[k]);
688 const int numCols = xc;
691 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
692 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
694 int err = C[k]->reshape (numRows, numCols);
695 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
696 "IMGS orthogonalization: failed to reshape " 697 "C[" << k <<
"] (the array of block " 698 "coefficients resulting from projecting X " 699 "against Q[1:" << nq <<
"]).");
705 if (MX == Teuchos::null) {
707 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
708 OPT::Apply(*(this->_Op),X,*MX);
713 MX = Teuchos::rcp( &X,
false );
716 int mxc = MVT::GetNumberVecs( *MX );
717 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
720 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::ICGSOrthoManager::projectAndNormalize(): X must be non-empty" );
723 for (
int i=0; i<nq; i++) {
724 numbas += MVT::GetNumberVecs( *Q[i] );
728 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
729 "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
731 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
732 "Belos::ICGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
734 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
735 "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
741 bool dep_flg =
false;
747 dep_flg = blkOrtho1( X, MX, C, Q );
750 if ( B == Teuchos::null ) {
751 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
753 std::vector<ScalarType> diag(xc);
755 #ifdef BELOS_TEUCHOS_TIME_MONITOR 756 Teuchos::TimeMonitor normTimer( *timerNorm_ );
758 MVT::MvDot( X, *MX, diag );
760 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
762 if (SCT::magnitude((*B)(0,0)) > ZERO) {
764 MVT::MvScale( X, ONE/(*B)(0,0) );
767 MVT::MvScale( *MX, ONE/(*B)(0,0) );
774 Teuchos::RCP<MV> tmpX, tmpMX;
775 tmpX = MVT::CloneCopy(X);
777 tmpMX = MVT::CloneCopy(*MX);
781 dep_flg = blkOrtho( X, MX, C, Q );
787 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
790 MVT::Assign( *tmpX, X );
792 MVT::Assign( *tmpMX, *MX );
797 rank = findBasis( X, MX, B,
false );
802 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
805 MVT::Assign( *tmpX, X );
807 MVT::Assign( *tmpMX, *MX );
814 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
815 "Belos::ICGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
825 template<
class ScalarType,
class MV,
class OP>
827 MV &X, Teuchos::RCP<MV> MX,
828 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
830 #ifdef BELOS_TEUCHOS_TIME_MONITOR 831 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
835 return findBasis(X, MX, B,
true);
842 template<
class ScalarType,
class MV,
class OP>
844 MV &X, Teuchos::RCP<MV> MX,
845 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
846 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
862 #ifdef BELOS_TEUCHOS_TIME_MONITOR 863 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
866 int xc = MVT::GetNumberVecs( X );
867 ptrdiff_t xr = MVT::GetGlobalLength( X );
869 std::vector<int> qcs(nq);
871 if (nq == 0 || xc == 0 || xr == 0) {
874 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
883 if (MX == Teuchos::null) {
885 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
886 OPT::Apply(*(this->_Op),X,*MX);
891 MX = Teuchos::rcp( &X,
false );
893 int mxc = MVT::GetNumberVecs( *MX );
894 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
897 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
898 "Belos::ICGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
900 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
901 "Belos::ICGSOrthoManager::project(): Size of X not consistant with MX,Q" );
905 for (
int i=0; i<nq; i++) {
906 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
907 "Belos::ICGSOrthoManager::project(): Q lengths not mutually consistant" );
908 qcs[i] = MVT::GetNumberVecs( *Q[i] );
909 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
910 "Belos::ICGSOrthoManager::project(): Q has less rows than columns" );
914 if ( C[i] == Teuchos::null ) {
915 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
918 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
919 "Belos::ICGSOrthoManager::project(): Size of Q not consistant with size of C" );
924 blkOrtho( X, MX, C, Q );
931 template<
class ScalarType,
class MV,
class OP>
936 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
954 const ScalarType ONE = SCT::one ();
955 const MagnitudeType ZERO = SCT::magnitude (SCT::zero ());
957 const int xc = MVT::GetNumberVecs (X);
958 const ptrdiff_t xr = MVT::GetGlobalLength (X);
971 if (MX == Teuchos::null) {
973 MX = MVT::Clone(X,xc);
974 OPT::Apply(*(this->_Op),X,*MX);
981 if ( B == Teuchos::null ) {
982 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
985 const int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
986 const ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
989 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
990 "Belos::ICGSOrthoManager::findBasis(): X must be non-empty" );
991 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
992 "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
993 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
994 "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
995 TEUCHOS_TEST_FOR_EXCEPTION( as<ptrdiff_t> (xc) > xr, std::invalid_argument,
996 "Belos::ICGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
997 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
998 "Belos::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
1003 int xstart = xc - howMany;
1005 for (
int j = xstart; j < xc; j++) {
1011 bool addVec =
false;
1014 std::vector<int> index(1);
1016 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1017 Teuchos::RCP<MV> MXj;
1020 MXj = MVT::CloneViewNonConst( *MX, index );
1028 std::vector<int> prev_idx( numX );
1029 Teuchos::RCP<const MV> prevX, prevMX;
1030 Teuchos::RCP<MV> oldMXj;
1033 for (
int i=0; i<numX; i++) {
1036 prevX = MVT::CloneView( X, prev_idx );
1038 prevMX = MVT::CloneView( *MX, prev_idx );
1041 oldMXj = MVT::CloneCopy( *MXj );
1045 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1046 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1051 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1052 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1054 MVT::MvDot( *Xj, *MXj, oldDot );
1057 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
1058 "Belos::ICGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1062 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1064 for (
int i=0; i<max_ortho_steps_; ++i) {
1068 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1069 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1077 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1078 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1080 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1088 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1089 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1091 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1105 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1106 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1108 MVT::MvDot( *Xj, *oldMXj, newDot );
1111 newDot[0] = oldDot[0];
1115 if (completeBasis) {
1119 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1124 std::cout <<
"Belos::ICGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1127 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1128 Teuchos::RCP<MV> tempMXj;
1129 MVT::MvRandom( *tempXj );
1131 tempMXj = MVT::Clone( X, 1 );
1132 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1138 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1139 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1141 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1144 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1146 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1147 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1152 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1153 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1155 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1158 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1159 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1161 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1166 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1167 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1169 MVT::MvDot( *tempXj, *tempMXj, newDot );
1172 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1174 MVT::Assign( *tempXj, *Xj );
1176 MVT::Assign( *tempMXj, *MXj );
1188 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1196 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1197 if (SCT::magnitude(diag) > ZERO) {
1198 MVT::MvScale( *Xj, ONE/diag );
1201 MVT::MvScale( *MXj, ONE/diag );
1215 for (
int i=0; i<numX; i++) {
1216 (*B)(i,j) = product(i,0);
1227 template<
class ScalarType,
class MV,
class OP>
1229 ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1230 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1231 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1234 int xc = MVT::GetNumberVecs( X );
1235 const ScalarType ONE = SCT::one();
1237 std::vector<int> qcs( nq );
1238 for (
int i=0; i<nq; i++) {
1239 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1244 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1246 for (
int i=0; i<nq; i++) {
1249 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1250 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1256 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1257 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1259 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1265 OPT::Apply( *(this->_Op), X, *MX);
1269 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1270 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1272 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1273 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1275 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1282 for (
int j = 1; j < max_ortho_steps_; ++j) {
1284 for (
int i=0; i<nq; i++) {
1285 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1289 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1290 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1296 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1297 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1299 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1305 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1306 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1309 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1311 else if (xc <= qcs[i]) {
1313 OPT::Apply( *(this->_Op), X, *MX);
1324 template<
class ScalarType,
class MV,
class OP>
1326 ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1327 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1328 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1331 int xc = MVT::GetNumberVecs( X );
1332 bool dep_flg =
false;
1333 const ScalarType ONE = SCT::one();
1335 std::vector<int> qcs( nq );
1336 for (
int i=0; i<nq; i++) {
1337 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1343 std::vector<ScalarType> oldDot( xc );
1345 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1346 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1348 MVT::MvDot( X, *MX, oldDot );
1351 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1353 for (
int i=0; i<nq; i++) {
1356 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1357 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1363 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1364 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1366 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1371 OPT::Apply( *(this->_Op), X, *MX);
1375 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1376 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1378 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1379 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1381 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1388 for (
int j = 1; j < max_ortho_steps_; ++j) {
1390 for (
int i=0; i<nq; i++) {
1391 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1395 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1396 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1402 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1403 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1405 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1411 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1412 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1415 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1417 else if (xc <= qcs[i]) {
1419 OPT::Apply( *(this->_Op), X, *MX);
1426 std::vector<ScalarType> newDot(xc);
1428 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1429 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1431 MVT::MvDot( X, *MX, newDot );
1435 for (
int i=0; i<xc; i++){
1436 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1445 template<
class ScalarType,
class MV,
class OP>
1447 ICGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1448 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1449 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1450 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1452 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1454 const ScalarType ONE = SCT::one();
1455 const ScalarType ZERO = SCT::zero();
1458 int xc = MVT::GetNumberVecs( X );
1459 std::vector<int> indX( 1 );
1460 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1462 std::vector<int> qcs( nq );
1463 for (
int i=0; i<nq; i++) {
1464 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1468 Teuchos::RCP<const MV> lastQ;
1469 Teuchos::RCP<MV> Xj, MXj;
1470 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1473 for (
int j=0; j<xc; j++) {
1475 bool dep_flg =
false;
1479 std::vector<int> index( j );
1480 for (
int ind=0; ind<j; ind++) {
1483 lastQ = MVT::CloneView( X, index );
1486 Q.push_back( lastQ );
1488 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1493 Xj = MVT::CloneViewNonConst( X, indX );
1495 MXj = MVT::CloneViewNonConst( *MX, indX );
1503 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1504 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1506 MVT::MvDot( *Xj, *MXj, oldDot );
1509 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1511 for (
int i=0; i<Q.size(); i++) {
1514 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1518 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1519 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1524 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1525 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1528 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1533 OPT::Apply( *(this->_Op), *Xj, *MXj);
1537 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1538 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1540 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1541 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1543 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1550 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1552 for (
int i=0; i<Q.size(); i++) {
1553 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1554 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1558 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1559 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1565 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1566 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1568 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1575 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1576 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1578 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1580 else if (xc <= qcs[i]) {
1582 OPT::Apply( *(this->_Op), *Xj, *MXj);
1591 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1592 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1594 MVT::MvDot( *Xj, *MXj, newDot );
1598 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1604 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1606 MVT::MvScale( *Xj, ONE/diag );
1609 MVT::MvScale( *MXj, ONE/diag );
1617 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1618 Teuchos::RCP<MV> tempMXj;
1619 MVT::MvRandom( *tempXj );
1621 tempMXj = MVT::Clone( X, 1 );
1622 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1628 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1629 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1631 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1634 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1636 for (
int i=0; i<Q.size(); i++) {
1637 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1641 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1642 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1647 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1648 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1650 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1656 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1657 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1660 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1662 else if (xc <= qcs[i]) {
1664 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1673 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1674 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1676 MVT::MvDot( *tempXj, *tempMXj, newDot );
1680 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1681 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1687 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1689 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1709 template<
class ScalarType,
class MV,
class OP>
1712 using Teuchos::ParameterList;
1713 using Teuchos::parameterList;
1716 RCP<ParameterList> params = parameterList (
"ICGS");
1721 "Maximum number of orthogonalization passes (includes the " 1722 "first). Default is 2, since \"twice is enough\" for Krylov " 1725 "Block reorthogonalization threshold.");
1727 "Singular block detection threshold.");
1732 template<
class ScalarType,
class MV,
class OP>
1735 using Teuchos::ParameterList;
1738 RCP<ParameterList> params = getICGSDefaultParameters<ScalarType, MV, OP>();
1740 params->set (
"maxNumOrthogPasses",
1742 params->set (
"blkTol",
1744 params->set (
"singTol",
1752 #endif // BELOS_ICGS_ORTHOMANAGER_HPP void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
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.
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...
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
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.
static const int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
Teuchos::RCP< Teuchos::ParameterList > getICGSDefaultParameters()
"Default" parameters for robustness and accuracy.
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.
Declaration of basic traits for the multivector type.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
int normalize(MV &X, Teuchos::RCP< MV > MX, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method takes a multivector X and attempts to compute an orthonormal basis for ...
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
Class which defines basic traits for the operator type.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
Traits class which defines basic operations on multivectors.
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
void project(MV &X, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
This method calls project(X,Teuchos::null,C,Q); see documentation for that function.
~ICGSOrthoManager()
Destructor.
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 ...
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
MagnitudeType getSingTol() const
Return parameter for singular block detection.
Class which defines basic traits for the operator type.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Provides the inner product defining the orthogonality concepts, using the provided operator...
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (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 ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
void project(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a list of (mutually and internally) orthonormal bases Q, this method takes a multivector X and ...
Teuchos::RCP< Teuchos::ParameterList > getICGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...