47 #ifndef BELOS_DGKS_ORTHOMANAGER_HPP 48 #define BELOS_DGKS_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,
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)
139 #ifdef BELOS_TEUCHOS_TIME_MONITOR 140 std::stringstream ss;
143 std::string orthoLabel = ss.str() +
": Orthogonalization";
144 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
146 std::string updateLabel = ss.str() +
": Ortho (Update)";
147 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
149 std::string normLabel = ss.str() +
": Ortho (Norm)";
150 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
152 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
153 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
167 using Teuchos::ParameterList;
168 using Teuchos::parameterList;
172 RCP<ParameterList> params;
173 if (plist.is_null()) {
175 params = parameterList (*defaultParams);
178 params->validateParametersAndSetDefaults (*defaultParams);
186 const int maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
196 this->setMyParamList (params);
199 Teuchos::RCP<const Teuchos::ParameterList>
203 defaultParams_ = Belos::getDGKSDefaultParameters<ScalarType, MV, OP>();
217 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
218 if (! params.is_null()) {
223 params->set (
"blkTol", blk_tol);
231 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
232 if (! params.is_null()) {
233 params->set (
"depTol", dep_tol);
241 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
242 if (! params.is_null()) {
243 params->set (
"singTol", sing_tol);
290 void project ( MV &X, Teuchos::RCP<MV> MX,
291 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
292 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
298 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
299 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
329 int normalize ( MV &X, Teuchos::RCP<MV> MX,
330 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
335 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
399 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
400 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
401 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
413 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
424 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
430 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
439 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
440 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
449 void setLabel(
const std::string& label);
493 #ifdef BELOS_TEUCHOS_TIME_MONITOR 494 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
495 #endif // BELOS_TEUCHOS_TIME_MONITOR 501 int findBasis(MV &X, Teuchos::RCP<MV> MX,
502 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
503 bool completeBasis,
int howMany = -1 )
const;
506 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
507 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
508 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
511 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
512 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
513 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
529 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
530 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
531 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
535 template<
class ScalarType,
class MV,
class OP>
538 template<
class ScalarType,
class MV,
class OP>
541 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
544 template<
class ScalarType,
class MV,
class OP>
547 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one()
548 / Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
551 template<
class ScalarType,
class MV,
class OP>
554 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
556 template<
class ScalarType,
class MV,
class OP>
559 template<
class ScalarType,
class MV,
class OP>
562 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
564 template<
class ScalarType,
class MV,
class OP>
567 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
569 template<
class ScalarType,
class MV,
class OP>
572 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
576 template<
class ScalarType,
class MV,
class OP>
579 if (label != label_) {
581 #ifdef BELOS_TEUCHOS_TIME_MONITOR 582 std::stringstream ss;
583 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
585 std::string orthoLabel = ss.str() +
": Orthogonalization";
586 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
588 std::string updateLabel = ss.str() +
": Ortho (Update)";
589 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
591 std::string normLabel = ss.str() +
": Ortho (Norm)";
592 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
594 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
595 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
602 template<
class ScalarType,
class MV,
class OP>
603 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
605 const ScalarType ONE = SCT::one();
606 int rank = MVT::GetNumberVecs(X);
607 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
609 #ifdef BELOS_TEUCHOS_TIME_MONITOR 610 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
614 for (
int i=0; i<rank; i++) {
617 return xTx.normFrobenius();
622 template<
class ScalarType,
class MV,
class OP>
623 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
625 int r1 = MVT::GetNumberVecs(X1);
626 int r2 = MVT::GetNumberVecs(X2);
627 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
629 #ifdef BELOS_TEUCHOS_TIME_MONITOR 630 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
634 return xTx.normFrobenius();
639 template<
class ScalarType,
class MV,
class OP>
644 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
645 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
646 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 648 using Teuchos::Array;
650 using Teuchos::is_null;
653 using Teuchos::SerialDenseMatrix;
654 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
655 typedef typename Array< RCP< const MV > >::size_type size_type;
657 #ifdef BELOS_TEUCHOS_TIME_MONITOR 658 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
661 ScalarType ONE = SCT::one();
665 int xc = MVT::GetNumberVecs( X );
666 ptrdiff_t xr = MVT::GetGlobalLength( X );
673 B = rcp (
new serial_dense_matrix_type (xc, xc));
683 for (size_type k = 0; k < nq; ++k)
685 const int numRows = MVT::GetNumberVecs (*Q[k]);
686 const int numCols = xc;
689 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
690 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
692 int err = C[k]->reshape (numRows, numCols);
693 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
694 "DGKS orthogonalization: failed to reshape " 695 "C[" << k <<
"] (the array of block " 696 "coefficients resulting from projecting X " 697 "against Q[1:" << nq <<
"]).");
703 if (MX == Teuchos::null) {
705 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
706 OPT::Apply(*(this->_Op),X,*MX);
711 MX = Teuchos::rcp( &X,
false );
714 int mxc = MVT::GetNumberVecs( *MX );
715 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
718 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
721 for (
int i=0; i<nq; i++) {
722 numbas += MVT::GetNumberVecs( *Q[i] );
726 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
727 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
729 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
730 "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
732 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
733 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
739 bool dep_flg =
false;
745 dep_flg = blkOrtho1( X, MX, C, Q );
748 if ( B == Teuchos::null ) {
749 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
751 std::vector<ScalarType> diag(1);
753 #ifdef BELOS_TEUCHOS_TIME_MONITOR 754 Teuchos::TimeMonitor normTimer( *timerNorm_ );
756 MVT::MvDot( X, *MX, diag );
758 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
760 if (SCT::magnitude((*B)(0,0)) > ZERO) {
762 MVT::MvScale( X, ONE/(*B)(0,0) );
765 MVT::MvScale( *MX, ONE/(*B)(0,0) );
772 Teuchos::RCP<MV> tmpX, tmpMX;
773 tmpX = MVT::CloneCopy(X);
775 tmpMX = MVT::CloneCopy(*MX);
779 dep_flg = blkOrtho( X, MX, C, Q );
784 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
787 MVT::Assign( *tmpX, X );
789 MVT::Assign( *tmpMX, *MX );
794 rank = findBasis( X, MX, B,
false );
798 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
801 MVT::Assign( *tmpX, X );
803 MVT::Assign( *tmpMX, *MX );
810 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
811 "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
821 template<
class ScalarType,
class MV,
class OP>
823 MV &X, Teuchos::RCP<MV> MX,
824 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
826 #ifdef BELOS_TEUCHOS_TIME_MONITOR 827 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
831 return findBasis(X, MX, B,
true);
838 template<
class ScalarType,
class MV,
class OP>
840 MV &X, Teuchos::RCP<MV> MX,
841 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
842 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
858 #ifdef BELOS_TEUCHOS_TIME_MONITOR 859 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
862 int xc = MVT::GetNumberVecs( X );
863 ptrdiff_t xr = MVT::GetGlobalLength( X );
865 std::vector<int> qcs(nq);
867 if (nq == 0 || xc == 0 || xr == 0) {
870 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
879 if (MX == Teuchos::null) {
881 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
882 OPT::Apply(*(this->_Op),X,*MX);
887 MX = Teuchos::rcp( &X,
false );
889 int mxc = MVT::GetNumberVecs( *MX );
890 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
893 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
894 "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
896 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
897 "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
901 for (
int i=0; i<nq; i++) {
902 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
903 "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
904 qcs[i] = MVT::GetNumberVecs( *Q[i] );
905 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
906 "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
910 if ( C[i] == Teuchos::null ) {
911 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
914 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
915 "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
920 blkOrtho( X, MX, C, Q );
927 template<
class ScalarType,
class MV,
class OP>
929 MV &X, Teuchos::RCP<MV> MX,
930 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
931 bool completeBasis,
int howMany )
const {
948 const ScalarType ONE = SCT::one();
951 int xc = MVT::GetNumberVecs( X );
952 ptrdiff_t xr = MVT::GetGlobalLength( X );
965 if (MX == Teuchos::null) {
967 MX = MVT::Clone(X,xc);
968 OPT::Apply(*(this->_Op),X,*MX);
975 if ( B == Teuchos::null ) {
976 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
979 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
980 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
983 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
984 "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
985 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
986 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
987 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
988 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
989 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
990 "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
991 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
992 "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
997 int xstart = xc - howMany;
999 for (
int j = xstart; j < xc; j++) {
1005 bool addVec =
false;
1008 std::vector<int> index(1);
1010 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1011 Teuchos::RCP<MV> MXj;
1012 if ((this->_hasOp)) {
1014 MXj = MVT::CloneViewNonConst( *MX, index );
1022 std::vector<int> prev_idx( numX );
1023 Teuchos::RCP<const MV> prevX, prevMX;
1024 Teuchos::RCP<MV> oldMXj;
1027 for (
int i=0; i<numX; i++) {
1030 prevX = MVT::CloneView( X, prev_idx );
1032 prevMX = MVT::CloneView( *MX, prev_idx );
1035 oldMXj = MVT::CloneCopy( *MXj );
1039 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1040 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1045 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1046 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1048 MVT::MvDot( *Xj, *MXj, oldDot );
1051 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO,
OrthoError,
1052 "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1059 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1060 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1067 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1068 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1070 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1078 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1079 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1081 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1086 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1087 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1089 MVT::MvDot( *Xj, *MXj, newDot );
1093 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1096 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1098 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1099 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1106 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1107 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1109 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1111 if ((this->_hasOp)) {
1112 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1113 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1115 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1123 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1124 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1126 MVT::MvDot( *Xj, *oldMXj, newDot );
1129 newDot[0] = oldDot[0];
1133 if (completeBasis) {
1137 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1142 std::cout <<
"Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1145 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1146 Teuchos::RCP<MV> tempMXj;
1147 MVT::MvRandom( *tempXj );
1149 tempMXj = MVT::Clone( X, 1 );
1150 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1156 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1157 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1159 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1162 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
1164 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1165 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1170 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1171 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1173 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1176 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1177 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1179 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1184 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1185 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1187 MVT::MvDot( *tempXj, *tempMXj, newDot );
1190 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1192 MVT::Assign( *tempXj, *Xj );
1194 MVT::Assign( *tempMXj, *MXj );
1206 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1214 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1216 if (SCT::magnitude(diag) > ZERO) {
1217 MVT::MvScale( *Xj, ONE/diag );
1220 MVT::MvScale( *MXj, ONE/diag );
1234 for (
int i=0; i<numX; i++) {
1235 (*B)(i,j) = product(i,0);
1246 template<
class ScalarType,
class MV,
class OP>
1249 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1250 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1253 int xc = MVT::GetNumberVecs( X );
1254 const ScalarType ONE = SCT::one();
1256 std::vector<int> qcs( nq );
1257 for (
int i=0; i<nq; i++) {
1258 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1262 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1264 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1265 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1267 MVT::MvDot( X, *MX, oldDot );
1270 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1272 for (
int i=0; i<nq; i++) {
1275 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1276 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1282 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1283 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1285 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1291 OPT::Apply( *(this->_Op), X, *MX);
1295 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1296 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1298 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1299 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1301 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1308 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1309 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1311 MVT::MvDot( X, *MX, newDot );
1325 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1328 for (
int i=0; i<nq; i++) {
1329 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1333 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1334 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1341 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1342 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1344 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1350 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1351 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1354 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1356 else if (xc <= qcs[i]) {
1358 OPT::Apply( *(this->_Op), X, *MX);
1369 template<
class ScalarType,
class MV,
class OP>
1372 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1373 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1376 int xc = MVT::GetNumberVecs( X );
1377 bool dep_flg =
false;
1378 const ScalarType ONE = SCT::one();
1380 std::vector<int> qcs( nq );
1381 for (
int i=0; i<nq; i++) {
1382 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1388 std::vector<ScalarType> oldDot( xc );
1390 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1391 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1393 MVT::MvDot( X, *MX, oldDot );
1396 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1398 for (
int i=0; i<nq; i++) {
1401 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1402 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1408 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1409 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1411 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1417 OPT::Apply( *(this->_Op), X, *MX);
1421 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1422 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1424 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1425 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1427 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1434 for (
int j = 1; j < max_blk_ortho_; ++j) {
1436 for (
int i=0; i<nq; i++) {
1437 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1441 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1442 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1449 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1450 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1452 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1458 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1459 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1462 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1464 else if (xc <= qcs[i]) {
1466 OPT::Apply( *(this->_Op), X, *MX);
1473 std::vector<ScalarType> newDot(xc);
1475 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1476 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1478 MVT::MvDot( X, *MX, newDot );
1482 for (
int i=0; i<xc; i++){
1483 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1493 template<
class ScalarType,
class MV,
class OP>
1496 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1497 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1498 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1500 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1502 const ScalarType ONE = SCT::one();
1503 const ScalarType ZERO = SCT::zero();
1506 int xc = MVT::GetNumberVecs( X );
1507 std::vector<int> indX( 1 );
1508 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1510 std::vector<int> qcs( nq );
1511 for (
int i=0; i<nq; i++) {
1512 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1516 Teuchos::RCP<const MV> lastQ;
1517 Teuchos::RCP<MV> Xj, MXj;
1518 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1521 for (
int j=0; j<xc; j++) {
1523 bool dep_flg =
false;
1527 std::vector<int> index( j );
1528 for (
int ind=0; ind<j; ind++) {
1531 lastQ = MVT::CloneView( X, index );
1534 Q.push_back( lastQ );
1536 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1541 Xj = MVT::CloneViewNonConst( X, indX );
1543 MXj = MVT::CloneViewNonConst( *MX, indX );
1551 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1552 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1554 MVT::MvDot( *Xj, *MXj, oldDot );
1557 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1559 for (
int i=0; i<Q.size(); i++) {
1562 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1566 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1567 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1573 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1574 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1576 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1582 OPT::Apply( *(this->_Op), *Xj, *MXj);
1586 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1587 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1589 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1590 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1592 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1600 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1601 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1603 MVT::MvDot( *Xj, *MXj, newDot );
1609 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
1611 for (
int i=0; i<Q.size(); i++) {
1612 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1613 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1617 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1618 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1624 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1625 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1627 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1633 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1634 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1637 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1639 else if (xc <= qcs[i]) {
1641 OPT::Apply( *(this->_Op), *Xj, *MXj);
1648 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1649 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1651 MVT::MvDot( *Xj, *MXj, newDot );
1656 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1662 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1664 MVT::MvScale( *Xj, ONE/diag );
1667 MVT::MvScale( *MXj, ONE/diag );
1675 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1676 Teuchos::RCP<MV> tempMXj;
1677 MVT::MvRandom( *tempXj );
1679 tempMXj = MVT::Clone( X, 1 );
1680 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1686 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1687 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1689 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1692 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
1694 for (
int i=0; i<Q.size(); i++) {
1695 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1699 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1700 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1705 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1706 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1708 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1714 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1715 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1718 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1720 else if (xc <= qcs[i]) {
1722 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1731 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1732 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1734 MVT::MvDot( *tempXj, *tempMXj, newDot );
1738 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1739 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1745 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1747 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1767 template<
class ScalarType,
class MV,
class OP>
1770 using Teuchos::ParameterList;
1771 using Teuchos::parameterList;
1774 RCP<ParameterList> params = parameterList (
"DGKS");
1779 "Maximum number of orthogonalization passes (includes the " 1780 "first). Default is 2, since \"twice is enough\" for Krylov " 1783 "Block reorthogonalization threshold.");
1785 "(Non-block) reorthogonalization threshold.");
1787 "Singular block detection threshold.");
1792 template<
class ScalarType,
class MV,
class OP>
1795 using Teuchos::ParameterList;
1798 RCP<ParameterList> params = getDGKSDefaultParameters<ScalarType, MV, OP>();
1800 params->set (
"maxNumOrthogPasses",
1802 params->set (
"blkTol",
1804 params->set (
"depTol",
1806 params->set (
"singTol",
1814 #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 ...