47 #ifndef BELOS_DGKS_ORTHOMANAGER_HPP 48 #define BELOS_DGKS_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,
110 #ifdef BELOS_TEUCHOS_TIME_MONITOR 111 std::stringstream ss;
114 std::string orthoLabel = ss.str() +
": Orthogonalization";
115 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
117 std::string updateLabel = ss.str() +
": Ortho (Update)";
118 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
120 std::string normLabel = ss.str() +
": Ortho (Norm)";
121 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
123 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
124 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
130 const std::string& label =
"Belos",
131 Teuchos::RCP<const OP> Op = Teuchos::null)
141 #ifdef BELOS_TEUCHOS_TIME_MONITOR 142 std::stringstream ss;
145 std::string orthoLabel = ss.str() +
": Orthogonalization";
146 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
148 std::string updateLabel = ss.str() +
": Ortho (Update)";
149 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
151 std::string normLabel = ss.str() +
": Ortho (Norm)";
152 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
154 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
155 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
169 using Teuchos::ParameterList;
170 using Teuchos::parameterList;
174 RCP<ParameterList> params;
175 if (plist.is_null()) {
177 params = parameterList (*defaultParams);
180 params->validateParametersAndSetDefaults (*defaultParams);
188 const int maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
198 setMyParamList (params);
201 Teuchos::RCP<const Teuchos::ParameterList>
205 defaultParams_ = Belos::getDGKSDefaultParameters<ScalarType, MV, OP>();
219 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
220 if (! params.is_null()) {
225 params->set (
"blkTol", blk_tol);
233 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
234 if (! params.is_null()) {
235 params->set (
"depTol", dep_tol);
243 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
244 if (! params.is_null()) {
245 params->set (
"singTol", sing_tol);
292 void project ( MV &X, Teuchos::RCP<MV> MX,
293 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
294 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
300 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
301 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
331 int normalize ( MV &X, Teuchos::RCP<MV> MX,
332 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
337 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
401 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
402 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
403 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
415 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
426 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
432 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
441 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
442 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
451 void setLabel(
const std::string& label);
495 #ifdef BELOS_TEUCHOS_TIME_MONITOR 496 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
497 #endif // BELOS_TEUCHOS_TIME_MONITOR 503 int findBasis(MV &X, Teuchos::RCP<MV> MX,
504 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
505 bool completeBasis,
int howMany = -1 )
const;
508 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
509 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
510 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
513 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
514 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
515 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
531 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
532 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
533 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
537 template<
class ScalarType,
class MV,
class OP>
540 template<
class ScalarType,
class MV,
class OP>
543 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
546 template<
class ScalarType,
class MV,
class OP>
549 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one()
550 / Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
553 template<
class ScalarType,
class MV,
class OP>
556 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
558 template<
class ScalarType,
class MV,
class OP>
561 template<
class ScalarType,
class MV,
class OP>
564 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
566 template<
class ScalarType,
class MV,
class OP>
569 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
571 template<
class ScalarType,
class MV,
class OP>
574 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
578 template<
class ScalarType,
class MV,
class OP>
581 if (label != label_) {
583 #ifdef BELOS_TEUCHOS_TIME_MONITOR 584 std::stringstream ss;
585 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
587 std::string orthoLabel = ss.str() +
": Orthogonalization";
588 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
590 std::string updateLabel = ss.str() +
": Ortho (Update)";
591 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
593 std::string normLabel = ss.str() +
": Ortho (Norm)";
594 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
596 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
597 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
604 template<
class ScalarType,
class MV,
class OP>
605 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
607 const ScalarType ONE = SCT::one();
608 int rank = MVT::GetNumberVecs(X);
609 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
611 #ifdef BELOS_TEUCHOS_TIME_MONITOR 612 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
616 for (
int i=0; i<rank; i++) {
619 return xTx.normFrobenius();
624 template<
class ScalarType,
class MV,
class OP>
625 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
627 int r1 = MVT::GetNumberVecs(X1);
628 int r2 = MVT::GetNumberVecs(X2);
629 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
631 #ifdef BELOS_TEUCHOS_TIME_MONITOR 632 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
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();
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 "DGKS 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::DGKSOrthoManager::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::DGKSOrthoManager::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::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
734 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
735 "Belos::DGKSOrthoManager::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(1);
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 );
786 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
789 MVT::Assign( *tmpX, X );
791 MVT::Assign( *tmpMX, *MX );
796 rank = findBasis( X, MX, B,
false );
800 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
803 MVT::Assign( *tmpX, X );
805 MVT::Assign( *tmpMX, *MX );
812 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
813 "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
823 template<
class ScalarType,
class MV,
class OP>
825 MV &X, Teuchos::RCP<MV> MX,
826 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
828 #ifdef BELOS_TEUCHOS_TIME_MONITOR 829 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
833 return findBasis(X, MX, B,
true);
840 template<
class ScalarType,
class MV,
class OP>
842 MV &X, Teuchos::RCP<MV> MX,
843 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
844 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
860 #ifdef BELOS_TEUCHOS_TIME_MONITOR 861 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
864 int xc = MVT::GetNumberVecs( X );
865 ptrdiff_t xr = MVT::GetGlobalLength( X );
867 std::vector<int> qcs(nq);
869 if (nq == 0 || xc == 0 || xr == 0) {
872 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
881 if (MX == Teuchos::null) {
883 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
884 OPT::Apply(*(this->_Op),X,*MX);
889 MX = Teuchos::rcp( &X,
false );
891 int mxc = MVT::GetNumberVecs( *MX );
892 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
895 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
896 "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
898 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
899 "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
903 for (
int i=0; i<nq; i++) {
904 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
905 "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
906 qcs[i] = MVT::GetNumberVecs( *Q[i] );
907 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
908 "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
912 if ( C[i] == Teuchos::null ) {
913 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
916 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
917 "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
922 blkOrtho( X, MX, C, Q );
929 template<
class ScalarType,
class MV,
class OP>
931 MV &X, Teuchos::RCP<MV> MX,
932 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
933 bool completeBasis,
int howMany )
const {
950 const ScalarType ONE = SCT::one();
953 int xc = MVT::GetNumberVecs( X );
954 ptrdiff_t xr = MVT::GetGlobalLength( X );
967 if (MX == Teuchos::null) {
969 MX = MVT::Clone(X,xc);
970 OPT::Apply(*(this->_Op),X,*MX);
977 if ( B == Teuchos::null ) {
978 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
981 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
982 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
985 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
986 "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
987 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
988 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
989 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
990 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
991 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
992 "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
993 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
994 "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
999 int xstart = xc - howMany;
1001 for (
int j = xstart; j < xc; j++) {
1007 bool addVec =
false;
1010 std::vector<int> index(1);
1012 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1013 Teuchos::RCP<MV> MXj;
1014 if ((this->_hasOp)) {
1016 MXj = MVT::CloneViewNonConst( *MX, index );
1024 std::vector<int> prev_idx( numX );
1025 Teuchos::RCP<const MV> prevX, prevMX;
1026 Teuchos::RCP<MV> oldMXj;
1029 for (
int i=0; i<numX; i++) {
1032 prevX = MVT::CloneView( X, prev_idx );
1034 prevMX = MVT::CloneView( *MX, prev_idx );
1037 oldMXj = MVT::CloneCopy( *MXj );
1041 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1042 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1047 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1048 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1050 MVT::MvDot( *Xj, *MXj, oldDot );
1053 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO,
OrthoError,
1054 "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1061 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1062 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1069 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1070 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1072 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1080 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1081 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1083 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1088 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1089 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1091 MVT::MvDot( *Xj, *MXj, newDot );
1095 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1098 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1100 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1101 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1108 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1109 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1111 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1113 if ((this->_hasOp)) {
1114 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1115 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1117 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1125 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1126 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1128 MVT::MvDot( *Xj, *oldMXj, newDot );
1131 newDot[0] = oldDot[0];
1135 if (completeBasis) {
1139 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1144 std::cout <<
"Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1147 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1148 Teuchos::RCP<MV> tempMXj;
1149 MVT::MvRandom( *tempXj );
1151 tempMXj = MVT::Clone( X, 1 );
1152 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1158 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1159 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1161 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1164 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
1166 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1167 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1172 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1173 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1175 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1178 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1179 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1181 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1186 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1187 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1189 MVT::MvDot( *tempXj, *tempMXj, newDot );
1192 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1194 MVT::Assign( *tempXj, *Xj );
1196 MVT::Assign( *tempMXj, *MXj );
1208 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1216 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1218 if (SCT::magnitude(diag) > ZERO) {
1219 MVT::MvScale( *Xj, ONE/diag );
1222 MVT::MvScale( *MXj, ONE/diag );
1236 for (
int i=0; i<numX; i++) {
1237 (*B)(i,j) = product(i,0);
1248 template<
class ScalarType,
class MV,
class OP>
1251 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1252 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1255 int xc = MVT::GetNumberVecs( X );
1256 const ScalarType ONE = SCT::one();
1258 std::vector<int> qcs( nq );
1259 for (
int i=0; i<nq; i++) {
1260 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1264 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1266 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1267 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1269 MVT::MvDot( X, *MX, oldDot );
1272 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1274 for (
int i=0; i<nq; i++) {
1277 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1278 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1284 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1285 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1287 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1293 OPT::Apply( *(this->_Op), X, *MX);
1297 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1298 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1300 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1301 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1303 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1310 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1311 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1313 MVT::MvDot( X, *MX, newDot );
1327 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1330 for (
int i=0; i<nq; i++) {
1331 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1335 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1336 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1343 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1344 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1346 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1352 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1353 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1356 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1358 else if (xc <= qcs[i]) {
1360 OPT::Apply( *(this->_Op), X, *MX);
1371 template<
class ScalarType,
class MV,
class OP>
1374 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1375 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1378 int xc = MVT::GetNumberVecs( X );
1379 bool dep_flg =
false;
1380 const ScalarType ONE = SCT::one();
1382 std::vector<int> qcs( nq );
1383 for (
int i=0; i<nq; i++) {
1384 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1390 std::vector<ScalarType> oldDot( xc );
1392 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1393 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1395 MVT::MvDot( X, *MX, oldDot );
1398 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1400 for (
int i=0; i<nq; i++) {
1403 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1404 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1410 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1411 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1413 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1419 OPT::Apply( *(this->_Op), X, *MX);
1423 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1424 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1426 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1427 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1429 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1436 for (
int j = 1; j < max_blk_ortho_; ++j) {
1438 for (
int i=0; i<nq; i++) {
1439 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1443 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1444 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1451 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1452 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1454 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1460 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1461 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1464 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1466 else if (xc <= qcs[i]) {
1468 OPT::Apply( *(this->_Op), X, *MX);
1475 std::vector<ScalarType> newDot(xc);
1477 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1478 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1480 MVT::MvDot( X, *MX, newDot );
1484 for (
int i=0; i<xc; i++){
1485 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1495 template<
class ScalarType,
class MV,
class OP>
1498 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1499 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1500 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1502 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1504 const ScalarType ONE = SCT::one();
1505 const ScalarType ZERO = SCT::zero();
1508 int xc = MVT::GetNumberVecs( X );
1509 std::vector<int> indX( 1 );
1510 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1512 std::vector<int> qcs( nq );
1513 for (
int i=0; i<nq; i++) {
1514 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1518 Teuchos::RCP<const MV> lastQ;
1519 Teuchos::RCP<MV> Xj, MXj;
1520 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1523 for (
int j=0; j<xc; j++) {
1525 bool dep_flg =
false;
1529 std::vector<int> index( j );
1530 for (
int ind=0; ind<j; ind++) {
1533 lastQ = MVT::CloneView( X, index );
1536 Q.push_back( lastQ );
1538 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1543 Xj = MVT::CloneViewNonConst( X, indX );
1545 MXj = MVT::CloneViewNonConst( *MX, indX );
1553 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1554 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1556 MVT::MvDot( *Xj, *MXj, oldDot );
1559 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1561 for (
int i=0; i<Q.size(); i++) {
1564 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1568 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1569 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1575 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1576 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1578 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1584 OPT::Apply( *(this->_Op), *Xj, *MXj);
1588 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1589 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1591 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1592 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1594 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1602 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1603 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1605 MVT::MvDot( *Xj, *MXj, newDot );
1611 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
1613 for (
int i=0; i<Q.size(); i++) {
1614 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1615 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1619 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1620 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1626 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1627 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1629 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1635 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1636 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1639 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1641 else if (xc <= qcs[i]) {
1643 OPT::Apply( *(this->_Op), *Xj, *MXj);
1650 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1651 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1653 MVT::MvDot( *Xj, *MXj, newDot );
1658 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1664 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1666 MVT::MvScale( *Xj, ONE/diag );
1669 MVT::MvScale( *MXj, ONE/diag );
1677 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1678 Teuchos::RCP<MV> tempMXj;
1679 MVT::MvRandom( *tempXj );
1681 tempMXj = MVT::Clone( X, 1 );
1682 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1688 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1689 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1691 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1694 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
1696 for (
int i=0; i<Q.size(); i++) {
1697 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1701 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1702 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1707 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1708 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1710 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1716 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1717 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1720 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1722 else if (xc <= qcs[i]) {
1724 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1733 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1734 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1736 MVT::MvDot( *tempXj, *tempMXj, newDot );
1740 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1741 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1747 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1749 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1769 template<
class ScalarType,
class MV,
class OP>
1772 using Teuchos::ParameterList;
1773 using Teuchos::parameterList;
1776 RCP<ParameterList> params = parameterList (
"DGKS");
1781 "Maximum number of orthogonalization passes (includes the " 1782 "first). Default is 2, since \"twice is enough\" for Krylov " 1785 "Block reorthogonalization threshold.");
1787 "(Non-block) reorthogonalization threshold.");
1789 "Singular block detection threshold.");
1794 template<
class ScalarType,
class MV,
class OP>
1797 using Teuchos::ParameterList;
1800 RCP<ParameterList> params = getDGKSDefaultParameters<ScalarType, MV, OP>();
1802 params->set (
"maxNumOrthogPasses",
1804 params->set (
"blkTol",
1806 params->set (
"depTol",
1808 params->set (
"singTol",
1816 #endif // BELOS_DGKS_ORTHOMANAGER_HPP void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
static const MagnitudeType dep_tol_fast_
(Non-block) reorthogonalization threshold (fast).
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 sing_tol_default_
Singular block detection threshold (default).
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
OperatorTraits< ScalarType, MV, OP > OPT
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
MagnitudeType sing_tol_
Singular block detection threshold.
Teuchos::RCP< Teuchos::ParameterList > getDGKSDefaultParameters()
"Default" parameters for robustness and accuracy.
DGKSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_blk_ortho=max_blk_ortho_default_, const MagnitudeType blk_tol=blk_tol_default_, const MagnitudeType dep_tol=dep_tol_default_, const MagnitudeType sing_tol=sing_tol_default_)
Constructor specifying re-orthogonalization tolerance.
Teuchos::ScalarTraits< ScalarType > SCT
~DGKSOrthoManager()
Destructor.
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
static const MagnitudeType dep_tol_default_
(Non-block) reorthogonalization threshold (default).
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
Declaration of basic traits for the multivector type.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
DGKSOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &plist, const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor that takes a list of parameters.
Class which defines basic traits for the operator type.
static const int max_blk_ortho_default_
Max number of (re)orthogonalization steps, including the first (default).
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Traits class which defines basic operations on multivectors.
MagnitudeType dep_tol_
(Non-block) reorthogonalization threshold.
Teuchos::ScalarTraits< MagnitudeType > MGT
std::string label_
Label for timer(s).
Exception thrown to signal error in an orthogonalization manager method.
MagnitudeType blk_tol_
Block reorthogonalization threshold.
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.
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.
static const int max_blk_ortho_fast_
Max number of (re)orthogonalization steps, including the first (fast).
MagnitudeType getDepTol() const
Return parameter for re-orthogonalization threshhold.
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
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.
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.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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 setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
int max_blk_ortho_
Max number of (re)orthogonalization steps, including the first.
Class which defines basic traits for the operator type.
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.
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...
Teuchos::RCP< Teuchos::ParameterList > getDGKSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
Belos header file which uses auto-configuration information to include necessary C++ headers...
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.
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...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
MultiVecTraits< ScalarType, MV > MVT
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 ...