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,
102 max_blk_ortho_( max_blk_ortho ),
105 sing_tol_( sing_tol ),
108 #ifdef BELOS_TEUCHOS_TIME_MONITOR 109 std::stringstream ss;
110 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
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;
141 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
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");
187 const MagnitudeType blkTol = params->get<MagnitudeType> (
"blkTol");
188 const MagnitudeType depTol = params->get<MagnitudeType> (
"depTol");
189 const MagnitudeType singTol = params->get<MagnitudeType> (
"singTol");
191 max_blk_ortho_ = maxNumOrthogPasses;
196 this->setMyParamList (params);
199 Teuchos::RCP<const Teuchos::ParameterList>
202 if (defaultParams_.is_null()) {
203 defaultParams_ = Belos::getDGKSDefaultParameters<ScalarType, MV, OP>();
206 return defaultParams_;
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);
245 sing_tol_ = 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);
453 const std::string&
getLabel()
const {
return label_; }
485 MagnitudeType blk_tol_;
487 MagnitudeType dep_tol_;
489 MagnitudeType sing_tol_;
493 #ifdef BELOS_TEUCHOS_TIME_MONITOR 494 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
495 #endif // BELOS_TEUCHOS_TIME_MONITOR 498 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
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;
528 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
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>
539 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
541 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
542 Teuchos::ScalarTraits<
typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps() );
544 template<
class ScalarType,
class MV,
class OP>
545 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
547 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one()
548 / Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
549 2*Teuchos::ScalarTraits<
typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one() );
551 template<
class ScalarType,
class MV,
class OP>
552 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
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>
560 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
562 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
564 template<
class ScalarType,
class MV,
class OP>
565 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
567 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
569 template<
class ScalarType,
class MV,
class OP>
570 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
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();
662 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
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();
949 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
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>
1248 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
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>
1371 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
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>
1495 DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
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.
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
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.
~DGKSOrthoManager()
Destructor.
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.
Traits class which defines basic operations on multivectors.
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.
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.
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...
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...
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...
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 ...
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 ...