47 #ifndef BELOS_ICGS_ORTHOMANAGER_HPP 48 #define BELOS_ICGS_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 72 template<
class ScalarType,
class MV,
class OP>
75 public Teuchos::ParameterListAcceptorDefaultBase
78 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
79 typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
80 typedef Teuchos::ScalarTraits<ScalarType> SCT;
90 Teuchos::RCP<const OP> Op = Teuchos::null,
91 const int max_ortho_steps = 2,
92 const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
93 const MagnitudeType sing_tol = 10*MGT::eps() )
95 max_ortho_steps_( max_ortho_steps ),
97 sing_tol_( sing_tol ),
100 #ifdef BELOS_TEUCHOS_TIME_MONITOR 101 std::string orthoLabel = label_ +
": Orthogonalization";
102 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
104 std::string updateLabel = label_ +
": Ortho (Update)";
105 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
107 std::string normLabel = label_ +
": Ortho (Norm)";
108 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
110 std::string ipLabel = label_ +
": Ortho (Inner Product)";
111 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
117 const std::string& label =
"Belos",
118 Teuchos::RCP<const OP> Op = Teuchos::null) :
120 max_ortho_steps_ (2),
121 blk_tol_ (10 * MGT::squareroot (MGT::eps())),
122 sing_tol_ (10 * MGT::eps()),
127 #ifdef BELOS_TEUCHOS_TIME_MONITOR 128 std::string orthoLabel = label_ +
": Orthogonalization";
129 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
131 std::string updateLabel = label_ +
": Ortho (Update)";
132 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
134 std::string normLabel = label_ +
": Ortho (Norm)";
135 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
137 std::string ipLabel = label_ +
": Ortho (Inner Product)";
138 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
152 using Teuchos::Exceptions::InvalidParameterName;
153 using Teuchos::ParameterList;
154 using Teuchos::parameterList;
158 RCP<ParameterList> params;
159 if (plist.is_null()) {
160 params = parameterList (*defaultParams);
175 int maxNumOrthogPasses;
176 MagnitudeType blkTol;
177 MagnitudeType singTol;
180 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
181 }
catch (InvalidParameterName&) {
182 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
183 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
194 blkTol = params->get<MagnitudeType> (
"blkTol");
195 }
catch (InvalidParameterName&) {
197 blkTol = params->get<MagnitudeType> (
"depTol");
200 params->remove (
"depTol");
201 }
catch (InvalidParameterName&) {
202 blkTol = defaultParams->get<MagnitudeType> (
"blkTol");
204 params->set (
"blkTol", blkTol);
208 singTol = params->get<MagnitudeType> (
"singTol");
209 }
catch (InvalidParameterName&) {
210 singTol = defaultParams->get<MagnitudeType> (
"singTol");
211 params->set (
"singTol", singTol);
214 max_ortho_steps_ = maxNumOrthogPasses;
218 setMyParamList (params);
221 Teuchos::RCP<const Teuchos::ParameterList>
225 using Teuchos::ParameterList;
226 using Teuchos::parameterList;
229 if (defaultParams_.is_null()) {
230 RCP<ParameterList> params = parameterList (
"ICGS");
234 const int defaultMaxNumOrthogPasses = 2;
235 const MagnitudeType eps = MGT::eps();
236 const MagnitudeType defaultBlkTol =
237 as<MagnitudeType> (10) * MGT::squareroot (eps);
238 const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
240 params->set (
"maxNumOrthogPasses", defaultMaxNumOrthogPasses,
241 "Maximum number of orthogonalization passes (includes the " 242 "first). Default is 2, since \"twice is enough\" for Krylov " 244 params->set (
"blkTol", defaultBlkTol,
"Block reorthogonalization " 246 params->set (
"singTol", defaultSingTol,
"Singular block detection " 248 defaultParams_ = params;
250 return defaultParams_;
258 Teuchos::RCP<const Teuchos::ParameterList>
262 using Teuchos::ParameterList;
263 using Teuchos::parameterList;
268 RCP<ParameterList> params = parameterList (*defaultParams);
270 const int maxBlkOrtho = 1;
271 const MagnitudeType blkTol = MGT::zero();
272 const MagnitudeType singTol = MGT::zero();
274 params->set (
"maxNumOrthogPasses", maxBlkOrtho);
275 params->set (
"blkTol", blkTol);
276 params->set (
"singTol", singTol);
287 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
288 if (! params.is_null()) {
293 params->set (
"blkTol", blk_tol);
301 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
302 if (! params.is_null()) {
307 params->set (
"singTol", sing_tol);
309 sing_tol_ = sing_tol;
351 void project ( MV &X, Teuchos::RCP<MV> MX,
352 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
353 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
359 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
360 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
390 int normalize ( MV &X, Teuchos::RCP<MV> MX,
391 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
396 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
446 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
447 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
448 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
458 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
467 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
473 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
482 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
483 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
492 void setLabel(
const std::string& label);
496 const std::string&
getLabel()
const {
return label_; }
502 int max_ortho_steps_;
504 MagnitudeType blk_tol_;
506 MagnitudeType sing_tol_;
510 #ifdef BELOS_TEUCHOS_TIME_MONITOR 511 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_,
512 timerNorm_, timerScale_, timerInnerProd_;
513 #endif // BELOS_TEUCHOS_TIME_MONITOR 516 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
519 int findBasis(MV &X, Teuchos::RCP<MV> MX,
520 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
521 bool completeBasis,
int howMany = -1 )
const;
524 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
525 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
526 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
529 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
530 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
531 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
546 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
547 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
548 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
549 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
554 template<
class ScalarType,
class MV,
class OP>
557 if (label != label_) {
559 std::string orthoLabel = label_ +
": Orthogonalization";
560 #ifdef BELOS_TEUCHOS_TIME_MONITOR 561 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
564 std::string updateLabel = label_ +
": Ortho (Update)";
565 #ifdef BELOS_TEUCHOS_TIME_MONITOR 566 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
569 std::string normLabel = label_ +
": Ortho (Norm)";
570 #ifdef BELOS_TEUCHOS_TIME_MONITOR 571 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
574 std::string ipLabel = label_ +
": Ortho (Inner Product)";
575 #ifdef BELOS_TEUCHOS_TIME_MONITOR 576 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
583 template<
class ScalarType,
class MV,
class OP>
584 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
586 const ScalarType ONE = SCT::one();
587 int rank = MVT::GetNumberVecs(X);
588 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
590 for (
int i=0; i<rank; i++) {
593 return xTx.normFrobenius();
598 template<
class ScalarType,
class MV,
class OP>
599 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
601 int r1 = MVT::GetNumberVecs(X1);
602 int r2 = MVT::GetNumberVecs(X2);
603 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
605 return xTx.normFrobenius();
610 template<
class ScalarType,
class MV,
class OP>
615 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
616 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
617 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 619 using Teuchos::Array;
621 using Teuchos::is_null;
624 using Teuchos::SerialDenseMatrix;
625 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
626 typedef typename Array< RCP< const MV > >::size_type size_type;
628 #ifdef BELOS_TEUCHOS_TIME_MONITOR 629 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
632 ScalarType ONE = SCT::one();
636 int xc = MVT::GetNumberVecs( X );
637 ptrdiff_t xr = MVT::GetGlobalLength( X );
644 B = rcp (
new serial_dense_matrix_type (xc, xc));
654 for (size_type k = 0; k < nq; ++k)
656 const int numRows = MVT::GetNumberVecs (*Q[k]);
657 const int numCols = xc;
660 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
661 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
663 int err = C[k]->reshape (numRows, numCols);
664 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
665 "IMGS orthogonalization: failed to reshape " 666 "C[" << k <<
"] (the array of block " 667 "coefficients resulting from projecting X " 668 "against Q[1:" << nq <<
"]).");
674 if (MX == Teuchos::null) {
676 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
677 OPT::Apply(*(this->_Op),X,*MX);
682 MX = Teuchos::rcp( &X,
false );
685 int mxc = MVT::GetNumberVecs( *MX );
686 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
689 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::ICGSOrthoManager::projectAndNormalize(): X must be non-empty" );
692 for (
int i=0; i<nq; i++) {
693 numbas += MVT::GetNumberVecs( *Q[i] );
697 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
698 "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
700 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
701 "Belos::ICGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
703 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
704 "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
710 bool dep_flg =
false;
713 Teuchos::RCP<MV> tmpX, tmpMX;
714 tmpX = MVT::CloneCopy(X);
716 tmpMX = MVT::CloneCopy(*MX);
723 dep_flg = blkOrtho1( X, MX, C, Q );
727 #ifdef BELOS_TEUCHOS_TIME_MONITOR 728 Teuchos::TimeMonitor normTimer( *timerNorm_ );
730 if ( B == Teuchos::null ) {
731 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
733 std::vector<ScalarType> diag(xc);
734 MVT::MvDot( X, *MX, diag );
735 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
737 MVT::MvScale( X, ONE/(*B)(0,0) );
740 MVT::MvScale( *MX, ONE/(*B)(0,0) );
748 dep_flg = blkOrtho( X, MX, C, Q );
754 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
757 MVT::Assign( *tmpX, X );
759 MVT::Assign( *tmpMX, *MX );
764 rank = findBasis( X, MX, B,
false );
769 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
772 MVT::Assign( *tmpX, X );
774 MVT::Assign( *tmpMX, *MX );
781 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
782 "Belos::ICGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
792 template<
class ScalarType,
class MV,
class OP>
794 MV &X, Teuchos::RCP<MV> MX,
795 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
797 #ifdef BELOS_TEUCHOS_TIME_MONITOR 798 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
802 return findBasis(X, MX, B,
true);
808 template<
class ScalarType,
class MV,
class OP>
810 MV &X, Teuchos::RCP<MV> MX,
811 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
812 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
828 #ifdef BELOS_TEUCHOS_TIME_MONITOR 829 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
832 int xc = MVT::GetNumberVecs( X );
833 ptrdiff_t xr = MVT::GetGlobalLength( X );
835 std::vector<int> qcs(nq);
837 if (nq == 0 || xc == 0 || xr == 0) {
840 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
849 if (MX == Teuchos::null) {
851 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
852 OPT::Apply(*(this->_Op),X,*MX);
857 MX = Teuchos::rcp( &X,
false );
859 int mxc = MVT::GetNumberVecs( *MX );
860 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
863 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
864 "Belos::ICGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
866 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
867 "Belos::ICGSOrthoManager::project(): Size of X not consistant with MX,Q" );
871 for (
int i=0; i<nq; i++) {
872 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
873 "Belos::ICGSOrthoManager::project(): Q lengths not mutually consistant" );
874 qcs[i] = MVT::GetNumberVecs( *Q[i] );
875 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
876 "Belos::ICGSOrthoManager::project(): Q has less rows than columns" );
880 if ( C[i] == Teuchos::null ) {
881 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
884 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
885 "Belos::ICGSOrthoManager::project(): Size of Q not consistant with size of C" );
890 blkOrtho( X, MX, C, Q );
897 template<
class ScalarType,
class MV,
class OP>
902 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
919 #ifdef BELOS_TEUCHOS_TIME_MONITOR 920 Teuchos::TimeMonitor normTimer (*timerNorm_);
921 #endif // BELOS_TEUCHOS_TIME_MONITOR 923 const ScalarType ONE = SCT::one ();
924 const MagnitudeType ZERO = SCT::magnitude (SCT::zero ());
926 const int xc = MVT::GetNumberVecs (X);
927 const ptrdiff_t xr = MVT::GetGlobalLength (X);
940 if (MX == Teuchos::null) {
942 MX = MVT::Clone(X,xc);
943 OPT::Apply(*(this->_Op),X,*MX);
950 if ( B == Teuchos::null ) {
951 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
954 const int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
955 const ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
958 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
959 "Belos::ICGSOrthoManager::findBasis(): X must be non-empty" );
960 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
961 "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
962 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
963 "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
964 TEUCHOS_TEST_FOR_EXCEPTION( as<ptrdiff_t> (xc) > xr, std::invalid_argument,
965 "Belos::ICGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
966 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
967 "Belos::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
972 int xstart = xc - howMany;
974 for (
int j = xstart; j < xc; j++) {
983 std::vector<int> index(1);
985 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
986 Teuchos::RCP<MV> MXj;
989 MXj = MVT::CloneViewNonConst( *MX, index );
997 std::vector<int> prev_idx( numX );
998 Teuchos::RCP<const MV> prevX, prevMX;
1001 for (
int i=0; i<numX; i++) {
1004 prevX = MVT::CloneView( X, prev_idx );
1006 prevMX = MVT::CloneView( *MX, prev_idx );
1011 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1012 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1016 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
1017 MVT::MvDot( *Xj, *MXj, oldDot );
1019 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
1020 "Belos::ICGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1024 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1026 for (
int i=0; i<max_ortho_steps_; ++i) {
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1031 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1039 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1040 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1042 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1050 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1051 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1053 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1066 MVT::MvDot( *Xj, *oldMXj, newDot );
1069 if (completeBasis) {
1073 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1078 std::cout <<
"Belos::ICGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1081 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1082 Teuchos::RCP<MV> tempMXj;
1083 MVT::MvRandom( *tempXj );
1085 tempMXj = MVT::Clone( X, 1 );
1086 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1091 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1093 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1095 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1096 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1101 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1102 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1104 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1107 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1108 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1110 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1114 MVT::MvDot( *tempXj, *tempMXj, newDot );
1116 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1118 MVT::Assign( *tempXj, *Xj );
1120 MVT::Assign( *tempMXj, *MXj );
1132 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1140 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1142 MVT::MvScale( *Xj, ONE/diag );
1145 MVT::MvScale( *MXj, ONE/diag );
1159 for (
int i=0; i<numX; i++) {
1160 (*B)(i,j) = product(i,0);
1171 template<
class ScalarType,
class MV,
class OP>
1173 ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1174 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1175 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1178 int xc = MVT::GetNumberVecs( X );
1179 const ScalarType ONE = SCT::one();
1181 std::vector<int> qcs( nq );
1182 for (
int i=0; i<nq; i++) {
1183 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1188 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1190 for (
int i=0; i<nq; i++) {
1193 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1194 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1200 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1201 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1203 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1209 OPT::Apply( *(this->_Op), X, *MX);
1213 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1214 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1215 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1221 for (
int j = 1; j < max_ortho_steps_; ++j) {
1223 for (
int i=0; i<nq; i++) {
1224 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
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, *Q[i], C2, ONE, X );
1245 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1247 else if (xc <= qcs[i]) {
1249 OPT::Apply( *(this->_Op), X, *MX);
1260 template<
class ScalarType,
class MV,
class OP>
1262 ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1263 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1264 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1267 int xc = MVT::GetNumberVecs( X );
1268 bool dep_flg =
false;
1269 const ScalarType ONE = SCT::one();
1271 std::vector<int> qcs( nq );
1272 for (
int i=0; i<nq; i++) {
1273 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1279 std::vector<ScalarType> oldDot( xc );
1280 MVT::MvDot( X, *MX, oldDot );
1282 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1284 for (
int i=0; i<nq; i++) {
1287 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1288 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1294 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1295 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1297 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1302 OPT::Apply( *(this->_Op), X, *MX);
1306 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1307 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1308 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1314 for (
int j = 1; j < max_ortho_steps_; ++j) {
1316 for (
int i=0; i<nq; i++) {
1317 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1321 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1322 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1328 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1329 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1331 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1338 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1340 else if (xc <= qcs[i]) {
1342 OPT::Apply( *(this->_Op), X, *MX);
1349 std::vector<ScalarType> newDot(xc);
1350 MVT::MvDot( X, *MX, newDot );
1353 for (
int i=0; i<xc; i++){
1354 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1363 template<
class ScalarType,
class MV,
class OP>
1365 ICGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1366 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1367 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1368 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1370 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1372 const ScalarType ONE = SCT::one();
1373 const ScalarType ZERO = SCT::zero();
1376 int xc = MVT::GetNumberVecs( X );
1377 std::vector<int> indX( 1 );
1378 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1380 std::vector<int> qcs( nq );
1381 for (
int i=0; i<nq; i++) {
1382 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1386 Teuchos::RCP<const MV> lastQ;
1387 Teuchos::RCP<MV> Xj, MXj;
1388 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1391 for (
int j=0; j<xc; j++) {
1393 bool dep_flg =
false;
1397 std::vector<int> index( j );
1398 for (
int ind=0; ind<j; ind++) {
1401 lastQ = MVT::CloneView( X, index );
1404 Q.push_back( lastQ );
1406 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1411 Xj = MVT::CloneViewNonConst( X, indX );
1413 MXj = MVT::CloneViewNonConst( *MX, indX );
1420 MVT::MvDot( *Xj, *MXj, oldDot );
1422 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1424 for (
int i=0; i<Q.size(); i++) {
1427 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1432 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1437 OPT::Apply( *(this->_Op), *Xj, *MXj);
1441 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1442 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1443 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1449 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1451 for (
int i=0; i<Q.size(); i++) {
1452 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1453 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1458 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1464 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1466 else if (xc <= qcs[i]) {
1468 OPT::Apply( *(this->_Op), *Xj, *MXj);
1476 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1482 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1484 MVT::MvScale( *Xj, ONE/diag );
1487 MVT::MvScale( *MXj, ONE/diag );
1495 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1496 Teuchos::RCP<MV> tempMXj;
1497 MVT::MvRandom( *tempXj );
1499 tempMXj = MVT::Clone( X, 1 );
1500 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1505 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1507 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1509 for (
int i=0; i<Q.size(); i++) {
1510 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1514 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1520 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1522 else if (xc <= qcs[i]) {
1524 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1532 MVT::MvDot( *tempXj, *tempMXj, newDot );
1535 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1536 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1542 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1544 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1566 #endif // BELOS_ICGS_ORTHOMANAGER_HPP void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
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::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...
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
ICGSOrthoManager(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.
Declaration of basic traits for the multivector type.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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 ...
Class which defines basic traits for the operator type.
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
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.
~ICGSOrthoManager()
Destructor.
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 setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
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::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector, measured as the Frobenius norm of ...
Belos header file which uses auto-configuration information to include necessary C++ headers...
ICGSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_ortho_steps=2, const MagnitudeType blk_tol=10 *MGT::squareroot(MGT::eps()), const MagnitudeType sing_tol=10 *MGT::eps())
Constructor specifying re-orthogonalization tolerance.
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 ...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...