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();
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
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();
644 B = rcp (
new serial_dense_matrix_type (xc, xc));
654 for (size_type k = 0; k < nq; ++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) {
682 MX = Teuchos::rcp( &X,
false );
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++) {
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;
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);
735 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
748 dep_flg = blkOrtho( X, MX, C, Q );
754 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
764 rank = findBasis( X, MX, B,
false );
769 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
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_);
835 std::vector<int> qcs(nq);
837 if (nq == 0 || xc == 0 || xr == 0) {
849 if (MX == Teuchos::null) {
857 MX = Teuchos::rcp( &X,
false );
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++) {
873 "Belos::ICGSOrthoManager::project(): Q lengths not mutually consistant" );
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 ());
940 if (MX == Teuchos::null) {
950 if ( B == Teuchos::null ) {
951 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
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);
986 Teuchos::RCP<MV> MXj;
997 std::vector<int> prev_idx( numX );
998 Teuchos::RCP<const MV> prevX, prevMX;
1001 for (
int i=0; i<numX; i++) {
1011 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1012 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
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_ );
1050 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1051 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
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;
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_ );
1107 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1108 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1116 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1132 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1140 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1159 for (
int i=0; i<numX; i++) {
1160 (*B)(i,j) = product(i,0);
1171 template<
class ScalarType,
class MV,
class OP>
1174 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1175 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1179 const ScalarType ONE = SCT::one();
1181 std::vector<int> qcs( nq );
1182 for (
int i=0; i<nq; 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_ );
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_ );
1247 else if (xc <= qcs[i]) {
1260 template<
class ScalarType,
class MV,
class OP>
1263 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1264 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 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++) {
1279 std::vector<ScalarType> oldDot( xc );
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_ );
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_ );
1340 else if (xc <= qcs[i]) {
1349 std::vector<ScalarType> newDot(xc);
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>
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();
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++) {
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++) {
1404 Q.push_back( lastQ );
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 );
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 );
1466 else if (xc <= qcs[i]) {
1476 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1482 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1495 Teuchos::RCP<MV> tempXj =
MVT::Clone( X, 1 );
1496 Teuchos::RCP<MV> tempMXj;
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 );
1522 else if (xc <= qcs[i]) {
1535 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1536 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
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.
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
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.
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
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::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...
Class which defines basic traits for the operator type.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void Assign(const MV &A, MV &mv)
mv := A
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
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...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
static void Apply(const OP &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Apply Op to x, putting the result into y.
Teuchos::RCP< const OP > _Op
~ICGSOrthoManager()
Destructor.
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 ...
Exception thrown to signal error in an orthogonalization manager method.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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.
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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.
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 ...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
Class which defines basic traits for the operator type.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Belos header file which uses auto-configuration information to include necessary C++ headers...
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 ...
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 ...
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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const