47 #ifndef BELOS_IMGS_ORTHOMANAGER_HPP 48 #define BELOS_IMGS_ORTHOMANAGER_HPP 63 #include "Teuchos_SerialDenseMatrix.hpp" 64 #include "Teuchos_SerialDenseVector.hpp" 66 #include "Teuchos_as.hpp" 67 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp" 68 #ifdef BELOS_TEUCHOS_TIME_MONITOR 69 #include "Teuchos_TimeMonitor.hpp" 70 #endif // BELOS_TEUCHOS_TIME_MONITOR 75 template<
class ScalarType,
class MV,
class OP>
79 template<
class ScalarType,
class MV,
class OP>
82 template<
class ScalarType,
class MV,
class OP>
85 public Teuchos::ParameterListAcceptorDefaultBase
88 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
89 typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
90 typedef Teuchos::ScalarTraits<ScalarType> SCT;
100 Teuchos::RCP<const OP> Op = Teuchos::null,
105 max_ortho_steps_( max_ortho_steps ),
107 sing_tol_( sing_tol ),
110 #ifdef BELOS_TEUCHOS_TIME_MONITOR 111 std::stringstream ss;
112 ss << label_ +
": IMGS[" << max_ortho_steps_ <<
"]";
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) :
140 #ifdef BELOS_TEUCHOS_TIME_MONITOR 141 std::stringstream ss;
142 ss << label_ +
": IMGS[" << max_ortho_steps_ <<
"]";
144 std::string orthoLabel = ss.str() +
": Orthogonalization";
145 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
147 std::string updateLabel = ss.str() +
": Ortho (Update)";
148 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
150 std::string normLabel = ss.str() +
": Ortho (Norm)";
151 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
153 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
154 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
167 using Teuchos::Exceptions::InvalidParameterName;
168 using Teuchos::ParameterList;
169 using Teuchos::parameterList;
173 RCP<ParameterList> params;
174 if (plist.is_null()) {
175 params = parameterList (*defaultParams);
190 int maxNumOrthogPasses;
191 MagnitudeType blkTol;
192 MagnitudeType singTol;
195 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
196 }
catch (InvalidParameterName&) {
197 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
198 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
209 blkTol = params->get<MagnitudeType> (
"blkTol");
210 }
catch (InvalidParameterName&) {
212 blkTol = params->get<MagnitudeType> (
"depTol");
215 params->remove (
"depTol");
216 }
catch (InvalidParameterName&) {
217 blkTol = defaultParams->get<MagnitudeType> (
"blkTol");
219 params->set (
"blkTol", blkTol);
223 singTol = params->get<MagnitudeType> (
"singTol");
224 }
catch (InvalidParameterName&) {
225 singTol = defaultParams->get<MagnitudeType> (
"singTol");
226 params->set (
"singTol", singTol);
229 max_ortho_steps_ = maxNumOrthogPasses;
233 setMyParamList (params);
236 Teuchos::RCP<const Teuchos::ParameterList>
239 if (defaultParams_.is_null()) {
240 defaultParams_ = Belos::getIMGSDefaultParameters<ScalarType, MV, OP>();
243 return defaultParams_;
252 void setBlkTol(
const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
255 void setSingTol(
const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
296 void project ( MV &X, Teuchos::RCP<MV> MX,
297 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
298 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
304 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
305 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
335 int normalize ( MV &X, Teuchos::RCP<MV> MX,
336 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
341 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
390 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
391 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
392 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
402 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
411 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
417 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
426 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
427 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
436 void setLabel(
const std::string& label);
440 const std::string&
getLabel()
const {
return label_; }
466 int max_ortho_steps_;
468 MagnitudeType blk_tol_;
470 MagnitudeType sing_tol_;
474 #ifdef BELOS_TEUCHOS_TIME_MONITOR 475 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
476 #endif // BELOS_TEUCHOS_TIME_MONITOR 479 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
482 int findBasis(MV &X, Teuchos::RCP<MV> MX,
483 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
484 bool completeBasis,
int howMany = -1 )
const;
487 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
488 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
489 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
492 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
493 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
494 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
509 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
510 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
511 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
512 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
516 template<
class ScalarType,
class MV,
class OP>
519 template<
class ScalarType,
class MV,
class OP>
520 const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
522 = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
523 Teuchos::ScalarTraits<
typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps() );
525 template<
class ScalarType,
class MV,
class OP>
526 const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
528 = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
530 template<
class ScalarType,
class MV,
class OP>
533 template<
class ScalarType,
class MV,
class OP>
534 const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
536 = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
538 template<
class ScalarType,
class MV,
class OP>
539 const typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
541 = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
545 template<
class ScalarType,
class MV,
class OP>
548 if (label != label_) {
550 #ifdef BELOS_TEUCHOS_TIME_MONITOR 551 std::stringstream ss;
552 ss << label_ +
": IMGS[" << max_ortho_steps_ <<
"]";
554 std::string orthoLabel = ss.str() +
": Orthogonalization";
555 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
557 std::string updateLabel = ss.str() +
": Ortho (Update)";
558 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
560 std::string normLabel = ss.str() +
": Ortho (Norm)";
561 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
563 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
564 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
571 template<
class ScalarType,
class MV,
class OP>
572 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
574 const ScalarType ONE = SCT::one();
575 int rank = MVT::GetNumberVecs(X);
576 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
578 for (
int i=0; i<rank; i++) {
581 return xTx.normFrobenius();
586 template<
class ScalarType,
class MV,
class OP>
587 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
589 int r1 = MVT::GetNumberVecs(X1);
590 int r2 = MVT::GetNumberVecs(X2);
591 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
593 return xTx.normFrobenius();
598 template<
class ScalarType,
class MV,
class OP>
603 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
604 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
605 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 607 using Teuchos::Array;
609 using Teuchos::is_null;
612 using Teuchos::SerialDenseMatrix;
613 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
614 typedef typename Array< RCP< const MV > >::size_type size_type;
616 #ifdef BELOS_TEUCHOS_TIME_MONITOR 617 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
620 ScalarType ONE = SCT::one();
621 const MagnitudeType ZERO = MGT::zero();
624 int xc = MVT::GetNumberVecs( X );
625 ptrdiff_t xr = MVT::GetGlobalLength( X );
632 B = rcp (
new serial_dense_matrix_type (xc, xc));
642 for (size_type k = 0; k < nq; ++k)
644 const int numRows = MVT::GetNumberVecs (*Q[k]);
645 const int numCols = xc;
648 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
649 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
651 int err = C[k]->reshape (numRows, numCols);
652 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
653 "IMGS orthogonalization: failed to reshape " 654 "C[" << k <<
"] (the array of block " 655 "coefficients resulting from projecting X " 656 "against Q[1:" << nq <<
"]).");
662 if (MX == Teuchos::null) {
664 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
665 OPT::Apply(*(this->_Op),X,*MX);
670 MX = Teuchos::rcp( &X,
false );
673 int mxc = MVT::GetNumberVecs( *MX );
674 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
677 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
680 for (
int i=0; i<nq; i++) {
681 numbas += MVT::GetNumberVecs( *Q[i] );
685 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
686 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
688 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
689 "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
691 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
692 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
698 bool dep_flg =
false;
701 Teuchos::RCP<MV> tmpX, tmpMX;
702 tmpX = MVT::CloneCopy(X);
704 tmpMX = MVT::CloneCopy(*MX);
711 dep_flg = blkOrtho1( X, MX, C, Q );
714 if ( B == Teuchos::null ) {
715 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
717 std::vector<ScalarType> diag(xc);
719 #ifdef BELOS_TEUCHOS_TIME_MONITOR 720 Teuchos::TimeMonitor normTimer( *timerNorm_ );
722 MVT::MvDot( X, *MX, diag );
724 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
726 if (SCT::magnitude((*B)(0,0)) > ZERO) {
728 MVT::MvScale( X, ONE/(*B)(0,0) );
731 MVT::MvScale( *MX, ONE/(*B)(0,0) );
738 dep_flg = blkOrtho( X, MX, C, Q );
744 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
747 MVT::Assign( *tmpX, X );
749 MVT::Assign( *tmpMX, *MX );
754 rank = findBasis( X, MX, B,
false );
759 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
762 MVT::Assign( *tmpX, X );
764 MVT::Assign( *tmpMX, *MX );
771 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
772 "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
782 template<
class ScalarType,
class MV,
class OP>
784 MV &X, Teuchos::RCP<MV> MX,
785 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
787 #ifdef BELOS_TEUCHOS_TIME_MONITOR 788 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
792 return findBasis(X, MX, B,
true);
798 template<
class ScalarType,
class MV,
class OP>
800 MV &X, Teuchos::RCP<MV> MX,
801 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
802 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
818 #ifdef BELOS_TEUCHOS_TIME_MONITOR 819 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
822 int xc = MVT::GetNumberVecs( X );
823 ptrdiff_t xr = MVT::GetGlobalLength( X );
825 std::vector<int> qcs(nq);
827 if (nq == 0 || xc == 0 || xr == 0) {
830 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
839 if (MX == Teuchos::null) {
841 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
842 OPT::Apply(*(this->_Op),X,*MX);
847 MX = Teuchos::rcp( &X,
false );
849 int mxc = MVT::GetNumberVecs( *MX );
850 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
853 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
854 "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
856 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
857 "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
861 for (
int i=0; i<nq; i++) {
862 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
863 "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
864 qcs[i] = MVT::GetNumberVecs( *Q[i] );
865 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
866 "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
870 if ( C[i] == Teuchos::null ) {
871 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
874 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
875 "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
880 blkOrtho( X, MX, C, Q );
887 template<
class ScalarType,
class MV,
class OP>
889 MV &X, Teuchos::RCP<MV> MX,
890 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
891 bool completeBasis,
int howMany )
const {
908 const ScalarType ONE = SCT::one();
909 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
911 int xc = MVT::GetNumberVecs( X );
912 ptrdiff_t xr = MVT::GetGlobalLength( X );
925 if (MX == Teuchos::null) {
927 MX = MVT::Clone(X,xc);
928 OPT::Apply(*(this->_Op),X,*MX);
935 if ( B == Teuchos::null ) {
936 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
939 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
940 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
943 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
944 "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
945 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
946 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
947 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
948 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
949 TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
950 "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
951 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
952 "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
957 int xstart = xc - howMany;
959 for (
int j = xstart; j < xc; j++) {
968 std::vector<int> index(1);
970 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
971 Teuchos::RCP<MV> MXj;
972 if ((this->_hasOp)) {
974 MXj = MVT::CloneViewNonConst( *MX, index );
981 Teuchos::RCP<MV> oldMXj;
983 oldMXj = MVT::CloneCopy( *MXj );
987 Teuchos::SerialDenseVector<int,ScalarType> product(numX);
988 Teuchos::SerialDenseVector<int,ScalarType> P2(1);
989 Teuchos::RCP<const MV> prevX, prevMX;
991 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
996 #ifdef BELOS_TEUCHOS_TIME_MONITOR 997 Teuchos::TimeMonitor normTimer( *timerNorm_ );
999 MVT::MvDot( *Xj, *MXj, oldDot );
1002 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
1003 "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1006 for (
int ii=0; ii<numX; ii++) {
1009 prevX = MVT::CloneView( X, index );
1011 prevMX = MVT::CloneView( *MX, index );
1014 for (
int i=0; i<max_ortho_steps_; ++i) {
1018 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1019 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1027 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1028 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1030 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1038 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1039 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1041 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1046 product[ii] = P2[0];
1048 product[ii] += P2[0];
1056 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1057 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1059 MVT::MvDot( *Xj, *oldMXj, newDot );
1062 newDot[0] = oldDot[0];
1066 if (completeBasis) {
1070 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1075 std::cout <<
"Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1078 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1079 Teuchos::RCP<MV> tempMXj;
1080 MVT::MvRandom( *tempXj );
1082 tempMXj = MVT::Clone( X, 1 );
1083 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1089 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1090 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1092 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1096 for (
int ii=0; ii<numX; ii++) {
1099 prevX = MVT::CloneView( X, index );
1101 prevMX = MVT::CloneView( *MX, index );
1104 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1106 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1107 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1112 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1113 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1115 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
1118 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1119 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1121 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
1126 product[ii] = P2[0];
1128 product[ii] += P2[0];
1134 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1135 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1137 MVT::MvDot( *tempXj, *tempMXj, newDot );
1140 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1142 MVT::Assign( *tempXj, *Xj );
1144 MVT::Assign( *tempMXj, *MXj );
1156 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1165 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1167 MVT::MvScale( *Xj, ONE/diag );
1170 MVT::MvScale( *MXj, ONE/diag );
1184 for (
int i=0; i<numX; i++) {
1185 (*B)(i,j) = product(i);
1196 template<
class ScalarType,
class MV,
class OP>
1198 IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1199 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1200 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1203 int xc = MVT::GetNumberVecs( X );
1204 const ScalarType ONE = SCT::one();
1206 std::vector<int> qcs( nq );
1207 for (
int i=0; i<nq; i++) {
1208 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1212 std::vector<int> index(1);
1213 Teuchos::RCP<const MV> tempQ;
1215 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1217 for (
int i=0; i<nq; i++) {
1220 for (
int ii=0; ii<qcs[i]; ii++) {
1223 tempQ = MVT::CloneView( *Q[i], index );
1224 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1228 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1229 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1235 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1236 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1238 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1245 OPT::Apply( *(this->_Op), X, *MX);
1249 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1250 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1251 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1257 for (
int j = 1; j < max_ortho_steps_; ++j) {
1259 for (
int i=0; i<nq; i++) {
1261 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1264 for (
int ii=0; ii<qcs[i]; ii++) {
1267 tempQ = MVT::CloneView( *Q[i], index );
1268 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1269 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1273 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1274 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1280 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1281 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1283 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1292 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1294 else if (xc <= qcs[i]) {
1296 OPT::Apply( *(this->_Op), X, *MX);
1307 template<
class ScalarType,
class MV,
class OP>
1309 IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1310 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1311 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1314 int xc = MVT::GetNumberVecs( X );
1315 bool dep_flg =
false;
1316 const ScalarType ONE = SCT::one();
1318 std::vector<int> qcs( nq );
1319 for (
int i=0; i<nq; i++) {
1320 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1326 std::vector<ScalarType> oldDot( xc );
1328 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1329 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1331 MVT::MvDot( X, *MX, oldDot );
1334 std::vector<int> index(1);
1335 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1336 Teuchos::RCP<const MV> tempQ;
1339 for (
int i=0; i<nq; i++) {
1342 for (
int ii=0; ii<qcs[i]; ii++) {
1345 tempQ = MVT::CloneView( *Q[i], index );
1346 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1350 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1351 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1357 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1358 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1360 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1367 OPT::Apply( *(this->_Op), X, *MX);
1371 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1372 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1373 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1379 for (
int j = 1; j < max_ortho_steps_; ++j) {
1381 for (
int i=0; i<nq; i++) {
1382 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1385 for (
int ii=0; ii<qcs[i]; ii++) {
1388 tempQ = MVT::CloneView( *Q[i], index );
1389 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1390 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1394 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1395 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1401 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1402 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1404 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1412 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1414 else if (xc <= qcs[i]) {
1416 OPT::Apply( *(this->_Op), X, *MX);
1423 std::vector<ScalarType> newDot(xc);
1425 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1426 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1428 MVT::MvDot( X, *MX, newDot );
1432 for (
int i=0; i<xc; i++){
1433 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1442 template<
class ScalarType,
class MV,
class OP>
1444 IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1445 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1446 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1447 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1449 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1451 const ScalarType ONE = SCT::one();
1452 const ScalarType ZERO = SCT::zero();
1455 int xc = MVT::GetNumberVecs( X );
1456 std::vector<int> indX( 1 );
1457 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1459 std::vector<int> qcs( nq );
1460 for (
int i=0; i<nq; i++) {
1461 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1465 Teuchos::RCP<const MV> lastQ;
1466 Teuchos::RCP<MV> Xj, MXj;
1467 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1470 for (
int j=0; j<xc; j++) {
1472 bool dep_flg =
false;
1476 std::vector<int> index( j );
1477 for (
int ind=0; ind<j; ind++) {
1480 lastQ = MVT::CloneView( X, index );
1483 Q.push_back( lastQ );
1485 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1490 Xj = MVT::CloneViewNonConst( X, indX );
1492 MXj = MVT::CloneViewNonConst( *MX, indX );
1500 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1501 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1503 MVT::MvDot( *Xj, *MXj, oldDot );
1506 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1507 Teuchos::RCP<const MV> tempQ;
1510 for (
int i=0; i<Q.size(); i++) {
1513 for (
int ii=0; ii<qcs[i]; ii++) {
1516 tempQ = MVT::CloneView( *Q[i], indX );
1518 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1524 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1530 OPT::Apply( *(this->_Op), *Xj, *MXj);
1534 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1535 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1536 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1537 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1543 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1545 for (
int i=0; i<Q.size(); i++) {
1546 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1549 for (
int ii=0; ii<qcs[i]; ii++) {
1552 tempQ = MVT::CloneView( *Q[i], indX );
1554 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1558 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1562 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1569 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1571 else if (xc <= qcs[i]) {
1573 OPT::Apply( *(this->_Op), *Xj, *MXj);
1581 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1587 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1589 MVT::MvScale( *Xj, ONE/diag );
1592 MVT::MvScale( *MXj, ONE/diag );
1600 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1601 Teuchos::RCP<MV> tempMXj;
1602 MVT::MvRandom( *tempXj );
1604 tempMXj = MVT::Clone( X, 1 );
1605 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1611 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1612 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1614 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1617 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1619 for (
int i=0; i<Q.size(); i++) {
1620 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1623 for (
int ii=0; ii<qcs[i]; ii++) {
1626 tempQ = MVT::CloneView( *Q[i], indX );
1627 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1631 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1639 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1641 else if (xc <= qcs[i]) {
1643 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1651 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1652 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1654 MVT::MvDot( *tempXj, *tempMXj, newDot );
1658 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1659 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1665 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1667 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1687 template<
class ScalarType,
class MV,
class OP>
1690 using Teuchos::ParameterList;
1691 using Teuchos::parameterList;
1694 RCP<ParameterList> params = parameterList (
"IMGS");
1699 "Maximum number of orthogonalization passes (includes the " 1700 "first). Default is 2, since \"twice is enough\" for Krylov " 1703 "Block reorthogonalization threshold.");
1705 "Singular block detection threshold.");
1710 template<
class ScalarType,
class MV,
class OP>
1713 using Teuchos::ParameterList;
1716 RCP<ParameterList> params = getIMGSDefaultParameters<ScalarType, MV, OP>();
1718 params->set (
"maxNumOrthogPasses",
1720 params->set (
"blkTol",
1722 params->set (
"singTol",
1730 #endif // BELOS_IMGS_ORTHOMANAGER_HPP void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector, measured as the Frobenius norm of ...
Declaration of basic traits for the multivector type.
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...
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 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 ...
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
IMGSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_ortho_steps=max_ortho_steps_default_, const MagnitudeType blk_tol=blk_tol_default_, const MagnitudeType sing_tol=sing_tol_default_)
Constructor specifying re-orthogonalization tolerance.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
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.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
Teuchos::RCP< Teuchos::ParameterList > getIMGSDefaultParameters()
"Default" parameters for robustness and accuracy.
static const int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
IMGSOrthoManager(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.
~IMGSOrthoManager()
Destructor.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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 ...
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...
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< Teuchos::ParameterList > getIMGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
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.
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.