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,
104 max_blk_ortho_( max_blk_ortho ),
107 sing_tol_( sing_tol ),
110 #ifdef BELOS_TEUCHOS_TIME_MONITOR 111 std::stringstream ss;
112 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
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;
143 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
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");
189 const MagnitudeType blkTol = params->get<MagnitudeType> (
"blkTol");
190 const MagnitudeType depTol = params->get<MagnitudeType> (
"depTol");
191 const MagnitudeType singTol = params->get<MagnitudeType> (
"singTol");
193 max_blk_ortho_ = maxNumOrthogPasses;
198 setMyParamList (params);
201 Teuchos::RCP<const Teuchos::ParameterList>
204 if (defaultParams_.is_null()) {
205 defaultParams_ = Belos::getDGKSDefaultParameters<ScalarType, MV, OP>();
208 return defaultParams_;
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);
247 sing_tol_ = 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);
455 const std::string&
getLabel()
const {
return label_; }
487 MagnitudeType blk_tol_;
489 MagnitudeType dep_tol_;
491 MagnitudeType sing_tol_;
495 #ifdef BELOS_TEUCHOS_TIME_MONITOR 496 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
497 #endif // BELOS_TEUCHOS_TIME_MONITOR 500 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
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;
530 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
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>
541 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
543 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
544 Teuchos::ScalarTraits<
typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps() );
546 template<
class ScalarType,
class MV,
class OP>
547 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
549 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one()
550 / Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
551 2*Teuchos::ScalarTraits<
typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one() );
553 template<
class ScalarType,
class MV,
class OP>
554 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
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>
562 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
564 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
566 template<
class ScalarType,
class MV,
class OP>
567 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
569 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
571 template<
class ScalarType,
class MV,
class OP>
572 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
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();
664 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
667 int xc = MVT::GetNumberVecs( X );
668 ptrdiff_t xr = MVT::GetGlobalLength( X );
675 B = rcp (
new serial_dense_matrix_type (xc, xc));
685 for (size_type k = 0; k < nq; ++k)
687 const int numRows = MVT::GetNumberVecs (*Q[k]);
688 const int numCols = xc;
691 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
692 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
694 int err = C[k]->reshape (numRows, numCols);
695 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
696 "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();
951 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
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>
1250 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
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>
1373 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
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>
1497 DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
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.
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 ...