47 #ifndef BELOS_ICGS_ORTHOMANAGER_HPP 48 #define BELOS_ICGS_ORTHOMANAGER_HPP 64 #include "Teuchos_as.hpp" 65 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp" 66 #ifdef BELOS_TEUCHOS_TIME_MONITOR 67 #include "Teuchos_TimeMonitor.hpp" 68 #endif // BELOS_TEUCHOS_TIME_MONITOR 73 template<
class ScalarType,
class MV,
class OP>
77 template<
class ScalarType,
class MV,
class OP>
80 template<
class ScalarType,
class MV,
class OP>
83 public Teuchos::ParameterListAcceptorDefaultBase
86 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
87 typedef typename Teuchos::ScalarTraits<MagnitudeType>
MGT;
88 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
98 Teuchos::RCP<const OP> Op = Teuchos::null,
108 #ifdef BELOS_TEUCHOS_TIME_MONITOR 109 std::stringstream ss;
112 std::string orthoLabel = ss.str() +
": Orthogonalization";
113 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
115 std::string updateLabel = ss.str() +
": Ortho (Update)";
116 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
118 std::string normLabel = ss.str() +
": Ortho (Norm)";
119 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
121 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
122 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
128 const std::string& label =
"Belos",
129 Teuchos::RCP<const OP> Op = Teuchos::null) :
138 #ifdef BELOS_TEUCHOS_TIME_MONITOR 139 std::stringstream ss;
142 std::string orthoLabel = ss.str() +
": Orthogonalization";
143 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
145 std::string updateLabel = ss.str() +
": Ortho (Update)";
146 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
148 std::string normLabel = ss.str() +
": Ortho (Norm)";
149 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
151 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
152 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
166 using Teuchos::Exceptions::InvalidParameterName;
167 using Teuchos::ParameterList;
168 using Teuchos::parameterList;
172 RCP<ParameterList> params;
173 if (plist.is_null()) {
174 params = parameterList (*defaultParams);
189 int maxNumOrthogPasses;
194 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
195 }
catch (InvalidParameterName&) {
196 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
197 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
209 }
catch (InvalidParameterName&) {
214 params->remove (
"depTol");
215 }
catch (InvalidParameterName&) {
218 params->set (
"blkTol", blkTol);
223 }
catch (InvalidParameterName&) {
225 params->set (
"singTol", singTol);
232 setMyParamList (params);
235 Teuchos::RCP<const Teuchos::ParameterList>
239 defaultParams_ = Belos::getICGSDefaultParameters<ScalarType, MV, OP>();
251 Teuchos::RCP<const Teuchos::ParameterList>
255 using Teuchos::ParameterList;
256 using Teuchos::parameterList;
261 RCP<ParameterList> params = parameterList (*defaultParams);
276 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
277 if (! params.is_null()) {
282 params->set (
"blkTol", blk_tol);
290 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
291 if (! params.is_null()) {
296 params->set (
"singTol", sing_tol);
340 void project ( MV &X, Teuchos::RCP<MV> MX,
341 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
342 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
348 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
349 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
379 int normalize ( MV &X, Teuchos::RCP<MV> MX,
380 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
385 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
435 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
436 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
437 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
447 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
456 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
462 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
471 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
472 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
481 void setLabel(
const std::string& label);
519 #ifdef BELOS_TEUCHOS_TIME_MONITOR 520 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
521 #endif // BELOS_TEUCHOS_TIME_MONITOR 527 int findBasis(MV &X, Teuchos::RCP<MV> MX,
528 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
529 bool completeBasis,
int howMany = -1 )
const;
532 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
533 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
534 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
537 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
538 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
539 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
555 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
556 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
557 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
561 template<
class ScalarType,
class MV,
class OP>
564 template<
class ScalarType,
class MV,
class OP>
567 = 10*Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
570 template<
class ScalarType,
class MV,
class OP>
573 = 10*Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
575 template<
class ScalarType,
class MV,
class OP>
578 template<
class ScalarType,
class MV,
class OP>
581 = Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
583 template<
class ScalarType,
class MV,
class OP>
586 = Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
590 template<
class ScalarType,
class MV,
class OP>
593 if (label != label_) {
595 #ifdef BELOS_TEUCHOS_TIME_MONITOR 596 std::stringstream ss;
597 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
599 std::string orthoLabel = ss.str() +
": Orthogonalization";
600 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
602 std::string updateLabel = ss.str() +
": Ortho (Update)";
603 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
605 std::string normLabel = ss.str() +
": Ortho (Norm)";
606 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
608 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
609 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
616 template<
class ScalarType,
class MV,
class OP>
617 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
619 const ScalarType ONE = SCT::one();
620 int rank = MVT::GetNumberVecs(X);
621 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
623 for (
int i=0; i<rank; i++) {
626 return xTx.normFrobenius();
631 template<
class ScalarType,
class MV,
class OP>
632 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
634 int r1 = MVT::GetNumberVecs(X1);
635 int r2 = MVT::GetNumberVecs(X2);
636 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
638 return xTx.normFrobenius();
643 template<
class ScalarType,
class MV,
class OP>
648 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
649 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
650 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 652 using Teuchos::Array;
654 using Teuchos::is_null;
657 using Teuchos::SerialDenseMatrix;
658 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
659 typedef typename Array< RCP< const MV > >::size_type size_type;
661 #ifdef BELOS_TEUCHOS_TIME_MONITOR 662 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
665 ScalarType ONE = SCT::one();
669 int xc = MVT::GetNumberVecs( X );
670 ptrdiff_t xr = MVT::GetGlobalLength( X );
677 B = rcp (
new serial_dense_matrix_type (xc, xc));
687 for (size_type k = 0; k < nq; ++k)
689 const int numRows = MVT::GetNumberVecs (*Q[k]);
690 const int numCols = xc;
693 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
694 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
696 int err = C[k]->reshape (numRows, numCols);
697 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
698 "IMGS orthogonalization: failed to reshape " 699 "C[" << k <<
"] (the array of block " 700 "coefficients resulting from projecting X " 701 "against Q[1:" << nq <<
"]).");
707 if (MX == Teuchos::null) {
709 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
710 OPT::Apply(*(this->_Op),X,*MX);
715 MX = Teuchos::rcp( &X,
false );
718 int mxc = MVT::GetNumberVecs( *MX );
719 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
722 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::ICGSOrthoManager::projectAndNormalize(): X must be non-empty" );
725 for (
int i=0; i<nq; i++) {
726 numbas += MVT::GetNumberVecs( *Q[i] );
730 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
731 "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
733 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
734 "Belos::ICGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
736 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
737 "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
743 bool dep_flg =
false;
749 dep_flg = blkOrtho1( X, MX, C, Q );
752 if ( B == Teuchos::null ) {
753 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
755 std::vector<ScalarType> diag(xc);
757 #ifdef BELOS_TEUCHOS_TIME_MONITOR 758 Teuchos::TimeMonitor normTimer( *timerNorm_ );
760 MVT::MvDot( X, *MX, diag );
762 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
764 if (SCT::magnitude((*B)(0,0)) > ZERO) {
766 MVT::MvScale( X, ONE/(*B)(0,0) );
769 MVT::MvScale( *MX, ONE/(*B)(0,0) );
776 Teuchos::RCP<MV> tmpX, tmpMX;
777 tmpX = MVT::CloneCopy(X);
779 tmpMX = MVT::CloneCopy(*MX);
783 dep_flg = blkOrtho( X, MX, C, Q );
789 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
792 MVT::Assign( *tmpX, X );
794 MVT::Assign( *tmpMX, *MX );
799 rank = findBasis( X, MX, B,
false );
804 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
807 MVT::Assign( *tmpX, X );
809 MVT::Assign( *tmpMX, *MX );
816 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
817 "Belos::ICGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
827 template<
class ScalarType,
class MV,
class OP>
829 MV &X, Teuchos::RCP<MV> MX,
830 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
832 #ifdef BELOS_TEUCHOS_TIME_MONITOR 833 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
837 return findBasis(X, MX, B,
true);
844 template<
class ScalarType,
class MV,
class OP>
846 MV &X, Teuchos::RCP<MV> MX,
847 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
848 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
864 #ifdef BELOS_TEUCHOS_TIME_MONITOR 865 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
868 int xc = MVT::GetNumberVecs( X );
869 ptrdiff_t xr = MVT::GetGlobalLength( X );
871 std::vector<int> qcs(nq);
873 if (nq == 0 || xc == 0 || xr == 0) {
876 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
885 if (MX == Teuchos::null) {
887 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
888 OPT::Apply(*(this->_Op),X,*MX);
893 MX = Teuchos::rcp( &X,
false );
895 int mxc = MVT::GetNumberVecs( *MX );
896 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
899 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
900 "Belos::ICGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
902 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
903 "Belos::ICGSOrthoManager::project(): Size of X not consistant with MX,Q" );
907 for (
int i=0; i<nq; i++) {
908 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
909 "Belos::ICGSOrthoManager::project(): Q lengths not mutually consistant" );
910 qcs[i] = MVT::GetNumberVecs( *Q[i] );
911 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
912 "Belos::ICGSOrthoManager::project(): Q has less rows than columns" );
916 if ( C[i] == Teuchos::null ) {
917 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
920 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
921 "Belos::ICGSOrthoManager::project(): Size of Q not consistant with size of C" );
926 blkOrtho( X, MX, C, Q );
933 template<
class ScalarType,
class MV,
class OP>
938 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
956 const ScalarType ONE = SCT::one ();
959 const int xc = MVT::GetNumberVecs (X);
960 const ptrdiff_t xr = MVT::GetGlobalLength (X);
973 if (MX == Teuchos::null) {
975 MX = MVT::Clone(X,xc);
976 OPT::Apply(*(this->_Op),X,*MX);
983 if ( B == Teuchos::null ) {
984 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
987 const int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
988 const ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
991 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
992 "Belos::ICGSOrthoManager::findBasis(): X must be non-empty" );
993 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
994 "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
995 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
996 "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
997 TEUCHOS_TEST_FOR_EXCEPTION( as<ptrdiff_t> (xc) > xr, std::invalid_argument,
998 "Belos::ICGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
999 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
1000 "Belos::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
1005 int xstart = xc - howMany;
1007 for (
int j = xstart; j < xc; j++) {
1013 bool addVec =
false;
1016 std::vector<int> index(1);
1018 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1019 Teuchos::RCP<MV> MXj;
1022 MXj = MVT::CloneViewNonConst( *MX, index );
1030 std::vector<int> prev_idx( numX );
1031 Teuchos::RCP<const MV> prevX, prevMX;
1032 Teuchos::RCP<MV> oldMXj;
1035 for (
int i=0; i<numX; i++) {
1038 prevX = MVT::CloneView( X, prev_idx );
1040 prevMX = MVT::CloneView( *MX, prev_idx );
1043 oldMXj = MVT::CloneCopy( *MXj );
1047 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1048 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1053 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1054 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1056 MVT::MvDot( *Xj, *MXj, oldDot );
1059 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO,
OrthoError,
1060 "Belos::ICGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1064 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1066 for (
int i=0; i<max_ortho_steps_; ++i) {
1070 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1071 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1079 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1080 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1082 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1090 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1091 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1093 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1107 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1108 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1110 MVT::MvDot( *Xj, *oldMXj, newDot );
1113 newDot[0] = oldDot[0];
1117 if (completeBasis) {
1121 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1126 std::cout <<
"Belos::ICGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1129 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1130 Teuchos::RCP<MV> tempMXj;
1131 MVT::MvRandom( *tempXj );
1133 tempMXj = MVT::Clone( X, 1 );
1134 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1140 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1141 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1143 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1146 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1148 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1149 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1154 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1155 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1157 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1160 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1161 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1163 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1168 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1169 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1171 MVT::MvDot( *tempXj, *tempMXj, newDot );
1174 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1176 MVT::Assign( *tempXj, *Xj );
1178 MVT::Assign( *tempMXj, *MXj );
1190 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1198 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1200 MVT::MvScale( *Xj, ONE/diag );
1203 MVT::MvScale( *MXj, ONE/diag );
1217 for (
int i=0; i<numX; i++) {
1218 (*B)(i,j) = product(i,0);
1229 template<
class ScalarType,
class MV,
class OP>
1232 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1233 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1236 int xc = MVT::GetNumberVecs( X );
1237 const ScalarType ONE = SCT::one();
1239 std::vector<int> qcs( nq );
1240 for (
int i=0; i<nq; i++) {
1241 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1246 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1248 for (
int i=0; i<nq; i++) {
1251 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1252 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1258 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1259 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1261 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1267 OPT::Apply( *(this->_Op), X, *MX);
1271 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1272 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1274 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1275 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1277 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1284 for (
int j = 1; j < max_ortho_steps_; ++j) {
1286 for (
int i=0; i<nq; i++) {
1287 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1291 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1292 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1298 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1299 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1301 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1307 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1308 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1311 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1313 else if (xc <= qcs[i]) {
1315 OPT::Apply( *(this->_Op), X, *MX);
1326 template<
class ScalarType,
class MV,
class OP>
1329 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1330 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1333 int xc = MVT::GetNumberVecs( X );
1334 bool dep_flg =
false;
1335 const ScalarType ONE = SCT::one();
1337 std::vector<int> qcs( nq );
1338 for (
int i=0; i<nq; i++) {
1339 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1345 std::vector<ScalarType> oldDot( xc );
1347 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1348 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1350 MVT::MvDot( X, *MX, oldDot );
1353 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1355 for (
int i=0; i<nq; i++) {
1358 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1359 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1365 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1366 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1368 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1373 OPT::Apply( *(this->_Op), X, *MX);
1377 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1378 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1380 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1381 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1383 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1390 for (
int j = 1; j < max_ortho_steps_; ++j) {
1392 for (
int i=0; i<nq; i++) {
1393 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1397 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1398 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1404 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1405 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1407 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1413 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1414 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1417 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1419 else if (xc <= qcs[i]) {
1421 OPT::Apply( *(this->_Op), X, *MX);
1428 std::vector<ScalarType> newDot(xc);
1430 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1431 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1433 MVT::MvDot( X, *MX, newDot );
1437 for (
int i=0; i<xc; i++){
1438 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1447 template<
class ScalarType,
class MV,
class OP>
1450 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1451 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1452 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1454 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1456 const ScalarType ONE = SCT::one();
1457 const ScalarType ZERO = SCT::zero();
1460 int xc = MVT::GetNumberVecs( X );
1461 std::vector<int> indX( 1 );
1462 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1464 std::vector<int> qcs( nq );
1465 for (
int i=0; i<nq; i++) {
1466 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1470 Teuchos::RCP<const MV> lastQ;
1471 Teuchos::RCP<MV> Xj, MXj;
1472 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1475 for (
int j=0; j<xc; j++) {
1477 bool dep_flg =
false;
1481 std::vector<int> index( j );
1482 for (
int ind=0; ind<j; ind++) {
1485 lastQ = MVT::CloneView( X, index );
1488 Q.push_back( lastQ );
1490 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1495 Xj = MVT::CloneViewNonConst( X, indX );
1497 MXj = MVT::CloneViewNonConst( *MX, indX );
1505 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1506 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1508 MVT::MvDot( *Xj, *MXj, oldDot );
1511 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1513 for (
int i=0; i<Q.size(); i++) {
1516 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1520 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1521 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1526 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1527 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1530 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1535 OPT::Apply( *(this->_Op), *Xj, *MXj);
1539 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1540 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1542 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1543 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1545 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1552 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1554 for (
int i=0; i<Q.size(); i++) {
1555 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1556 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1560 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1561 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1567 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1568 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1570 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1577 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1578 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1580 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1582 else if (xc <= qcs[i]) {
1584 OPT::Apply( *(this->_Op), *Xj, *MXj);
1592 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1598 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1600 MVT::MvScale( *Xj, ONE/diag );
1603 MVT::MvScale( *MXj, ONE/diag );
1611 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1612 Teuchos::RCP<MV> tempMXj;
1613 MVT::MvRandom( *tempXj );
1615 tempMXj = MVT::Clone( X, 1 );
1616 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1622 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1623 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1625 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1628 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1630 for (
int i=0; i<Q.size(); i++) {
1631 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1635 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1636 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1641 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1642 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1644 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1650 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1651 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1654 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1656 else if (xc <= qcs[i]) {
1658 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1667 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1668 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1670 MVT::MvDot( *tempXj, *tempMXj, newDot );
1674 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1675 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1681 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1683 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1703 template<
class ScalarType,
class MV,
class OP>
1706 using Teuchos::ParameterList;
1707 using Teuchos::parameterList;
1710 RCP<ParameterList> params = parameterList (
"ICGS");
1715 "Maximum number of orthogonalization passes (includes the " 1716 "first). Default is 2, since \"twice is enough\" for Krylov " 1719 "Block reorthogonalization threshold.");
1721 "Singular block detection threshold.");
1726 template<
class ScalarType,
class MV,
class OP>
1729 using Teuchos::ParameterList;
1732 RCP<ParameterList> params = getICGSDefaultParameters<ScalarType, MV, OP>();
1734 params->set (
"maxNumOrthogPasses",
1736 params->set (
"blkTol",
1738 params->set (
"singTol",
1746 #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).
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
Project X against QQ and normalize X, one vector at a time.
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.
MultiVecTraits< ScalarType, MV > MVT
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
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
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).
int max_ortho_steps_
Max number of (re)orthogonalization steps, including the first.
Class which defines basic traits for the operator type.
MagnitudeType sing_tol_
Singular block detection threshold.
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.
MagnitudeType blk_tol_
Block reorthogonalization threshold.
Exception thrown to signal error in an orthogonalization manager method.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::ScalarTraits< MagnitudeType > MGT
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 ...
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
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.
OperatorTraits< ScalarType, MV, OP > OPT
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
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.
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).
std::string label_
Label for timers.
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.
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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...