47 #ifndef BELOS_DGKS_ORTHOMANAGER_HPP
48 #define BELOS_DGKS_ORTHOMANAGER_HPP
65 #ifdef BELOS_TEUCHOS_TIME_MONITOR
67 #endif // BELOS_TEUCHOS_TIME_MONITOR
72 template<
class ScalarType,
class MV,
class OP>
76 template<
class ScalarType,
class MV,
class OP>
79 template<
class ScalarType,
class MV,
class OP>
102 max_blk_ortho_( max_blk_ortho ),
105 sing_tol_( sing_tol ),
108 #ifdef BELOS_TEUCHOS_TIME_MONITOR
109 std::stringstream ss;
110 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
112 std::string orthoLabel = ss.str() +
": Orthogonalization";
115 std::string updateLabel = ss.str() +
": Ortho (Update)";
118 std::string normLabel = ss.str() +
": Ortho (Norm)";
121 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
128 const std::string& label =
"Belos",
139 #ifdef BELOS_TEUCHOS_TIME_MONITOR
140 std::stringstream ss;
141 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
143 std::string orthoLabel = ss.str() +
": Orthogonalization";
146 std::string updateLabel = ss.str() +
": Ortho (Update)";
149 std::string normLabel = ss.str() +
": Ortho (Norm)";
152 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
168 using Teuchos::parameterList;
172 RCP<ParameterList> params;
175 params = parameterList (*defaultParams);
186 const int maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
187 const MagnitudeType blkTol = params->get<MagnitudeType> (
"blkTol");
188 const MagnitudeType depTol = params->get<MagnitudeType> (
"depTol");
189 const MagnitudeType singTol = params->get<MagnitudeType> (
"singTol");
191 max_blk_ortho_ = maxNumOrthogPasses;
202 if (defaultParams_.
is_null()) {
203 defaultParams_ = Belos::getDGKSDefaultParameters<ScalarType, MV, OP>();
206 return defaultParams_;
223 params->
set (
"blkTol", blk_tol);
233 params->
set (
"depTol", dep_tol);
243 params->
set (
"singTol", sing_tol);
245 sing_tol_ = sing_tol;
449 void setLabel(
const std::string& label);
453 const std::string&
getLabel()
const {
return label_; }
485 MagnitudeType blk_tol_;
487 MagnitudeType dep_tol_;
489 MagnitudeType sing_tol_;
493 #ifdef BELOS_TEUCHOS_TIME_MONITOR
495 #endif // BELOS_TEUCHOS_TIME_MONITOR
503 bool completeBasis,
int howMany = -1 )
const;
535 template<
class ScalarType,
class MV,
class OP>
538 template<
class ScalarType,
class MV,
class OP>
539 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
544 template<
class ScalarType,
class MV,
class OP>
545 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
551 template<
class ScalarType,
class MV,
class OP>
552 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
556 template<
class ScalarType,
class MV,
class OP>
559 template<
class ScalarType,
class MV,
class OP>
560 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
564 template<
class ScalarType,
class MV,
class OP>
565 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
569 template<
class ScalarType,
class MV,
class OP>
570 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
576 template<
class ScalarType,
class MV,
class OP>
579 if (label != label_) {
581 #ifdef BELOS_TEUCHOS_TIME_MONITOR
582 std::stringstream ss;
583 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
585 std::string orthoLabel = ss.str() +
": Orthogonalization";
588 std::string updateLabel = ss.str() +
": Ortho (Update)";
591 std::string normLabel = ss.str() +
": Ortho (Norm)";
594 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
602 template<
class ScalarType,
class MV,
class OP>
605 const ScalarType ONE = SCT::one();
606 int rank = MVT::GetNumberVecs(X);
609 #ifdef BELOS_TEUCHOS_TIME_MONITOR
614 for (
int i=0; i<rank; i++) {
622 template<
class ScalarType,
class MV,
class OP>
625 int r1 = MVT::GetNumberVecs(X1);
626 int r2 = MVT::GetNumberVecs(X2);
629 #ifdef BELOS_TEUCHOS_TIME_MONITOR
639 template<
class ScalarType,
class MV,
class OP>
655 typedef typename Array< RCP< const MV > >::size_type size_type;
657 #ifdef BELOS_TEUCHOS_TIME_MONITOR
661 ScalarType ONE = SCT::one();
662 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
665 int xc = MVT::GetNumberVecs( X );
666 ptrdiff_t xr = MVT::GetGlobalLength( X );
673 B =
rcp (
new serial_dense_matrix_type (xc, xc));
683 for (size_type k = 0; k < nq; ++k)
685 const int numRows = MVT::GetNumberVecs (*Q[k]);
686 const int numCols = xc;
689 C[k] =
rcp (
new serial_dense_matrix_type (numRows, numCols));
690 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
692 int err = C[k]->reshape (numRows, numCols);
694 "DGKS orthogonalization: failed to reshape "
695 "C[" << k <<
"] (the array of block "
696 "coefficients resulting from projecting X "
697 "against Q[1:" << nq <<
"]).");
703 if (MX == Teuchos::null) {
705 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
706 OPT::Apply(*(this->_Op),X,*MX);
714 int mxc = MVT::GetNumberVecs( *MX );
715 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
718 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
721 for (
int i=0; i<nq; i++) {
722 numbas += MVT::GetNumberVecs( *Q[i] );
727 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
730 "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
732 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
733 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
739 bool dep_flg =
false;
745 dep_flg = blkOrtho1( X, MX, C, Q );
748 if ( B == Teuchos::null ) {
751 std::vector<ScalarType> diag(1);
753 #ifdef BELOS_TEUCHOS_TIME_MONITOR
756 MVT::MvDot( X, *MX, diag );
758 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
760 if (SCT::magnitude((*B)(0,0)) > ZERO) {
762 MVT::MvScale( X, ONE/(*B)(0,0) );
765 MVT::MvScale( *MX, ONE/(*B)(0,0) );
773 tmpX = MVT::CloneCopy(X);
775 tmpMX = MVT::CloneCopy(*MX);
779 dep_flg = blkOrtho( X, MX, C, Q );
784 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
787 MVT::Assign( *tmpX, X );
789 MVT::Assign( *tmpMX, *MX );
794 rank = findBasis( X, MX, B,
false );
798 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
801 MVT::Assign( *tmpX, X );
803 MVT::Assign( *tmpMX, *MX );
810 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
811 "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
821 template<
class ScalarType,
class MV,
class OP>
826 #ifdef BELOS_TEUCHOS_TIME_MONITOR
831 return findBasis(X, MX, B,
true);
838 template<
class ScalarType,
class MV,
class OP>
858 #ifdef BELOS_TEUCHOS_TIME_MONITOR
862 int xc = MVT::GetNumberVecs( X );
863 ptrdiff_t xr = MVT::GetGlobalLength( X );
865 std::vector<int> qcs(nq);
867 if (nq == 0 || xc == 0 || xr == 0) {
870 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
879 if (MX == Teuchos::null) {
881 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
882 OPT::Apply(*(this->_Op),X,*MX);
889 int mxc = MVT::GetNumberVecs( *MX );
890 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
894 "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
897 "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
900 for (
int i=0; i<nq; i++) {
902 "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
903 qcs[i] = MVT::GetNumberVecs( *Q[i] );
905 "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
908 if ( C[i] == Teuchos::null ) {
913 "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
918 blkOrtho( X, MX, C, Q );
925 template<
class ScalarType,
class MV,
class OP>
929 bool completeBasis,
int howMany )
const {
946 const ScalarType ONE = SCT::one();
947 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
949 int xc = MVT::GetNumberVecs( X );
950 ptrdiff_t xr = MVT::GetGlobalLength( X );
963 if (MX == Teuchos::null) {
965 MX = MVT::Clone(X,xc);
966 OPT::Apply(*(this->_Op),X,*MX);
973 if ( B == Teuchos::null ) {
977 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
978 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
982 "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
984 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
986 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
988 "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
990 "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
995 int xstart = xc - howMany;
997 for (
int j = xstart; j < xc; j++) {
1003 bool addVec =
false;
1006 std::vector<int> index(1);
1010 if ((this->_hasOp)) {
1012 MXj = MVT::CloneViewNonConst( *MX, index );
1020 std::vector<int> prev_idx( numX );
1025 for (
int i=0; i<numX; i++) {
1028 prevX = MVT::CloneView( X, prev_idx );
1030 prevMX = MVT::CloneView( *MX, prev_idx );
1033 oldMXj = MVT::CloneCopy( *MXj );
1038 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1043 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1046 MVT::MvDot( *Xj, *MXj, oldDot );
1050 "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1057 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1065 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1068 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1076 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1079 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1084 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1087 MVT::MvDot( *Xj, *MXj, newDot );
1091 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1096 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1104 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1107 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1109 if ((this->_hasOp)) {
1110 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1113 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1121 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1124 MVT::MvDot( *Xj, *oldMXj, newDot );
1127 newDot[0] = oldDot[0];
1131 if (completeBasis) {
1135 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1140 std::cout <<
"Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1145 MVT::MvRandom( *tempXj );
1147 tempMXj = MVT::Clone( X, 1 );
1148 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1154 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1157 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1160 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
1162 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1168 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1171 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1174 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1177 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1182 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1185 MVT::MvDot( *tempXj, *tempMXj, newDot );
1188 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1190 MVT::Assign( *tempXj, *Xj );
1192 MVT::Assign( *tempMXj, *MXj );
1204 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1212 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1214 if (SCT::magnitude(diag) > ZERO) {
1215 MVT::MvScale( *Xj, ONE/diag );
1218 MVT::MvScale( *MXj, ONE/diag );
1232 for (
int i=0; i<numX; i++) {
1233 (*B)(i,j) = product(i,0);
1244 template<
class ScalarType,
class MV,
class OP>
1246 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X,
Teuchos::RCP<MV> MX,
1251 int xc = MVT::GetNumberVecs( X );
1252 const ScalarType ONE = SCT::one();
1254 std::vector<int> qcs( nq );
1255 for (
int i=0; i<nq; i++) {
1256 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1260 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1262 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1265 MVT::MvDot( X, *MX, oldDot );
1270 for (
int i=0; i<nq; i++) {
1273 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1280 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1283 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1289 OPT::Apply( *(this->_Op), X, *MX);
1293 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1294 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1296 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1299 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1306 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1309 MVT::MvDot( X, *MX, newDot );
1323 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1326 for (
int i=0; i<nq; i++) {
1331 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1339 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1342 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1348 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1352 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1354 else if (xc <= qcs[i]) {
1356 OPT::Apply( *(this->_Op), X, *MX);
1367 template<
class ScalarType,
class MV,
class OP>
1369 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X,
Teuchos::RCP<MV> MX,
1374 int xc = MVT::GetNumberVecs( X );
1375 bool dep_flg =
false;
1376 const ScalarType ONE = SCT::one();
1378 std::vector<int> qcs( nq );
1379 for (
int i=0; i<nq; i++) {
1380 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1386 std::vector<ScalarType> oldDot( xc );
1388 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1391 MVT::MvDot( X, *MX, oldDot );
1396 for (
int i=0; i<nq; i++) {
1399 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1406 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1409 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1415 OPT::Apply( *(this->_Op), X, *MX);
1419 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1420 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1422 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1425 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1432 for (
int j = 1; j < max_blk_ortho_; ++j) {
1434 for (
int i=0; i<nq; i++) {
1439 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1447 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1450 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1456 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1460 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1462 else if (xc <= qcs[i]) {
1464 OPT::Apply( *(this->_Op), X, *MX);
1471 std::vector<ScalarType> newDot(xc);
1473 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1476 MVT::MvDot( X, *MX, newDot );
1480 for (
int i=0; i<xc; i++){
1481 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1491 template<
class ScalarType,
class MV,
class OP>
1493 DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X,
Teuchos::RCP<MV> MX,
1500 const ScalarType ONE = SCT::one();
1501 const ScalarType ZERO = SCT::zero();
1504 int xc = MVT::GetNumberVecs( X );
1505 std::vector<int> indX( 1 );
1506 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1508 std::vector<int> qcs( nq );
1509 for (
int i=0; i<nq; i++) {
1510 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1519 for (
int j=0; j<xc; j++) {
1521 bool dep_flg =
false;
1525 std::vector<int> index( j );
1526 for (
int ind=0; ind<j; ind++) {
1529 lastQ = MVT::CloneView( X, index );
1532 Q.push_back( lastQ );
1534 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1539 Xj = MVT::CloneViewNonConst( X, indX );
1541 MXj = MVT::CloneViewNonConst( *MX, indX );
1549 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1552 MVT::MvDot( *Xj, *MXj, oldDot );
1557 for (
int i=0; i<Q.size(); i++) {
1564 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1571 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1574 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1580 OPT::Apply( *(this->_Op), *Xj, *MXj);
1584 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1585 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1587 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1590 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1598 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1601 MVT::MvDot( *Xj, *MXj, newDot );
1607 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
1609 for (
int i=0; i<Q.size(); i++) {
1615 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1622 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1625 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1631 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1635 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1637 else if (xc <= qcs[i]) {
1639 OPT::Apply( *(this->_Op), *Xj, *MXj);
1646 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1649 MVT::MvDot( *Xj, *MXj, newDot );
1654 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1660 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1662 MVT::MvScale( *Xj, ONE/diag );
1665 MVT::MvScale( *MXj, ONE/diag );
1675 MVT::MvRandom( *tempXj );
1677 tempMXj = MVT::Clone( X, 1 );
1678 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1684 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1687 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1690 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
1692 for (
int i=0; i<Q.size(); i++) {
1697 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1703 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1706 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1712 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1716 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1718 else if (xc <= qcs[i]) {
1720 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1729 #ifdef BELOS_TEUCHOS_TIME_MONITOR
1732 MVT::MvDot( *tempXj, *tempMXj, newDot );
1736 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1737 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1743 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1745 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1765 template<
class ScalarType,
class MV,
class OP>
1769 using Teuchos::parameterList;
1772 RCP<ParameterList> params = parameterList (
"DGKS");
1777 "Maximum number of orthogonalization passes (includes the "
1778 "first). Default is 2, since \"twice is enough\" for Krylov "
1781 "Block reorthogonalization threshold.");
1783 "(Non-block) reorthogonalization threshold.");
1785 "Singular block detection threshold.");
1790 template<
class ScalarType,
class MV,
class OP>
1796 RCP<ParameterList> params = getDGKSDefaultParameters<ScalarType, MV, OP>();
1798 params->set (
"maxNumOrthogPasses",
1800 params->set (
"blkTol",
1802 params->set (
"depTol",
1804 params->set (
"singTol",
1812 #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).
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
bool is_null(const boost::shared_ptr< T > &p)
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 ...
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.
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
~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.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
bool is_null(const std::shared_ptr< T > &p)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Declaration of basic traits for the multivector type.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
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.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
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...
MagnitudeType getDepTol() const
Return parameter for re-orthogonalization threshhold.
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_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void setMyParamList(const RCP< ParameterList > ¶mList)
static const int max_blk_ortho_fast_
Max number of (re)orthogonalization steps, including the first (fast).
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
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.
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...
RCP< ParameterList > getNonconstParameterList()
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
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 ...
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
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.
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.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
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.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...