47 #ifndef BELOS_DGKS_ORTHOMANAGER_HPP 48 #define BELOS_DGKS_ORTHOMANAGER_HPP 65 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp" 66 #ifdef BELOS_TEUCHOS_TIME_MONITOR 68 #endif // BELOS_TEUCHOS_TIME_MONITOR 73 template<
class ScalarType,
class MV,
class OP>
77 template<
class ScalarType,
class MV,
class OP>
80 template<
class ScalarType,
class MV,
class OP>
104 max_blk_ortho_( max_blk_ortho ),
107 sing_tol_( sing_tol ),
110 #ifdef BELOS_TEUCHOS_TIME_MONITOR 111 std::stringstream ss;
112 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
114 std::string orthoLabel = ss.str() +
": Orthogonalization";
117 std::string updateLabel = ss.str() +
": Ortho (Update)";
120 std::string normLabel = ss.str() +
": Ortho (Norm)";
123 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
130 const std::string& label =
"Belos",
141 #ifdef BELOS_TEUCHOS_TIME_MONITOR 142 std::stringstream ss;
143 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
145 std::string orthoLabel = ss.str() +
": Orthogonalization";
148 std::string updateLabel = ss.str() +
": Ortho (Update)";
151 std::string normLabel = ss.str() +
": Ortho (Norm)";
154 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
170 using Teuchos::parameterList;
174 RCP<ParameterList> params;
177 params = parameterList (*defaultParams);
188 const int maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
189 const MagnitudeType blkTol = params->get<MagnitudeType> (
"blkTol");
190 const MagnitudeType depTol = params->get<MagnitudeType> (
"depTol");
191 const MagnitudeType singTol = params->get<MagnitudeType> (
"singTol");
193 max_blk_ortho_ = maxNumOrthogPasses;
204 if (defaultParams_.is_null()) {
205 defaultParams_ = Belos::getDGKSDefaultParameters<ScalarType, MV, OP>();
208 return defaultParams_;
225 params->
set (
"blkTol", blk_tol);
235 params->
set (
"depTol", dep_tol);
245 params->
set (
"singTol", sing_tol);
247 sing_tol_ = sing_tol;
451 void setLabel(
const std::string& label);
455 const std::string&
getLabel()
const {
return label_; }
487 MagnitudeType blk_tol_;
489 MagnitudeType dep_tol_;
491 MagnitudeType sing_tol_;
495 #ifdef BELOS_TEUCHOS_TIME_MONITOR 497 #endif // BELOS_TEUCHOS_TIME_MONITOR 505 bool completeBasis,
int howMany = -1 )
const;
537 template<
class ScalarType,
class MV,
class OP>
540 template<
class ScalarType,
class MV,
class OP>
541 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
546 template<
class ScalarType,
class MV,
class OP>
547 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
553 template<
class ScalarType,
class MV,
class OP>
554 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
558 template<
class ScalarType,
class MV,
class OP>
561 template<
class ScalarType,
class MV,
class OP>
562 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
566 template<
class ScalarType,
class MV,
class OP>
567 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
571 template<
class ScalarType,
class MV,
class OP>
572 const typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType
578 template<
class ScalarType,
class MV,
class OP>
581 if (label != label_) {
583 #ifdef BELOS_TEUCHOS_TIME_MONITOR 584 std::stringstream ss;
585 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
587 std::string orthoLabel = ss.str() +
": Orthogonalization";
590 std::string updateLabel = ss.str() +
": Ortho (Update)";
593 std::string normLabel = ss.str() +
": Ortho (Norm)";
596 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
604 template<
class ScalarType,
class MV,
class OP>
607 const ScalarType ONE = SCT::one();
608 int rank = MVT::GetNumberVecs(X);
611 #ifdef BELOS_TEUCHOS_TIME_MONITOR 616 for (
int i=0; i<rank; i++) {
624 template<
class ScalarType,
class MV,
class OP>
627 int r1 = MVT::GetNumberVecs(X1);
628 int r2 = MVT::GetNumberVecs(X2);
631 #ifdef BELOS_TEUCHOS_TIME_MONITOR 641 template<
class ScalarType,
class MV,
class OP>
657 typedef typename Array< RCP< const MV > >::size_type size_type;
659 #ifdef BELOS_TEUCHOS_TIME_MONITOR 663 ScalarType ONE = SCT::one();
664 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
667 int xc = MVT::GetNumberVecs( X );
668 ptrdiff_t xr = MVT::GetGlobalLength( X );
675 B =
rcp (
new serial_dense_matrix_type (xc, xc));
685 for (size_type k = 0; k < nq; ++k)
687 const int numRows = MVT::GetNumberVecs (*Q[k]);
688 const int numCols = xc;
691 C[k] =
rcp (
new serial_dense_matrix_type (numRows, numCols));
692 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
694 int err = C[k]->reshape (numRows, numCols);
696 "DGKS orthogonalization: failed to reshape " 697 "C[" << k <<
"] (the array of block " 698 "coefficients resulting from projecting X " 699 "against Q[1:" << nq <<
"]).");
705 if (MX == Teuchos::null) {
707 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
708 OPT::Apply(*(this->_Op),X,*MX);
716 int mxc = MVT::GetNumberVecs( *MX );
717 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
720 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
723 for (
int i=0; i<nq; i++) {
724 numbas += MVT::GetNumberVecs( *Q[i] );
729 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
732 "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
735 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
741 bool dep_flg =
false;
747 dep_flg = blkOrtho1( X, MX, C, Q );
750 if ( B == Teuchos::null ) {
753 std::vector<ScalarType> diag(1);
755 #ifdef BELOS_TEUCHOS_TIME_MONITOR 758 MVT::MvDot( X, *MX, diag );
760 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
762 if (SCT::magnitude((*B)(0,0)) > ZERO) {
764 MVT::MvScale( X, ONE/(*B)(0,0) );
767 MVT::MvScale( *MX, ONE/(*B)(0,0) );
775 tmpX = MVT::CloneCopy(X);
777 tmpMX = MVT::CloneCopy(*MX);
781 dep_flg = blkOrtho( X, MX, C, Q );
786 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
789 MVT::Assign( *tmpX, X );
791 MVT::Assign( *tmpMX, *MX );
796 rank = findBasis( X, MX, B,
false );
800 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
803 MVT::Assign( *tmpX, X );
805 MVT::Assign( *tmpMX, *MX );
813 "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
823 template<
class ScalarType,
class MV,
class OP>
828 #ifdef BELOS_TEUCHOS_TIME_MONITOR 833 return findBasis(X, MX, B,
true);
840 template<
class ScalarType,
class MV,
class OP>
860 #ifdef BELOS_TEUCHOS_TIME_MONITOR 864 int xc = MVT::GetNumberVecs( X );
865 ptrdiff_t xr = MVT::GetGlobalLength( X );
867 std::vector<int> qcs(nq);
869 if (nq == 0 || xc == 0 || xr == 0) {
872 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
881 if (MX == Teuchos::null) {
883 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
884 OPT::Apply(*(this->_Op),X,*MX);
891 int mxc = MVT::GetNumberVecs( *MX );
892 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
896 "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
899 "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
903 for (
int i=0; i<nq; i++) {
905 "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
906 qcs[i] = MVT::GetNumberVecs( *Q[i] );
908 "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
912 if ( C[i] == Teuchos::null ) {
917 "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
922 blkOrtho( X, MX, C, Q );
929 template<
class ScalarType,
class MV,
class OP>
933 bool completeBasis,
int howMany )
const {
950 const ScalarType ONE = SCT::one();
951 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
953 int xc = MVT::GetNumberVecs( X );
954 ptrdiff_t xr = MVT::GetGlobalLength( X );
967 if (MX == Teuchos::null) {
969 MX = MVT::Clone(X,xc);
970 OPT::Apply(*(this->_Op),X,*MX);
977 if ( B == Teuchos::null ) {
981 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
982 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
986 "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
988 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
990 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
992 "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
994 "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
999 int xstart = xc - howMany;
1001 for (
int j = xstart; j < xc; j++) {
1007 bool addVec =
false;
1010 std::vector<int> index(1);
1014 if ((this->_hasOp)) {
1016 MXj = MVT::CloneViewNonConst( *MX, index );
1024 std::vector<int> prev_idx( numX );
1029 for (
int i=0; i<numX; i++) {
1032 prevX = MVT::CloneView( X, prev_idx );
1034 prevMX = MVT::CloneView( *MX, prev_idx );
1037 oldMXj = MVT::CloneCopy( *MXj );
1042 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1047 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1050 MVT::MvDot( *Xj, *MXj, oldDot );
1054 "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1061 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1069 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1072 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1080 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1083 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1088 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1091 MVT::MvDot( *Xj, *MXj, newDot );
1095 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1100 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1108 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1111 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1113 if ((this->_hasOp)) {
1114 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1117 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1125 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1128 MVT::MvDot( *Xj, *oldMXj, newDot );
1131 newDot[0] = oldDot[0];
1135 if (completeBasis) {
1139 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1144 std::cout <<
"Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1149 MVT::MvRandom( *tempXj );
1151 tempMXj = MVT::Clone( X, 1 );
1152 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1158 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1161 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1164 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
1166 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1172 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1175 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1178 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1181 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1186 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1189 MVT::MvDot( *tempXj, *tempMXj, newDot );
1192 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1194 MVT::Assign( *tempXj, *Xj );
1196 MVT::Assign( *tempMXj, *MXj );
1208 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1216 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1218 if (SCT::magnitude(diag) > ZERO) {
1219 MVT::MvScale( *Xj, ONE/diag );
1222 MVT::MvScale( *MXj, ONE/diag );
1236 for (
int i=0; i<numX; i++) {
1237 (*B)(i,j) = product(i,0);
1248 template<
class ScalarType,
class MV,
class OP>
1250 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X,
Teuchos::RCP<MV> MX,
1255 int xc = MVT::GetNumberVecs( X );
1256 const ScalarType ONE = SCT::one();
1258 std::vector<int> qcs( nq );
1259 for (
int i=0; i<nq; i++) {
1260 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1264 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1266 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1269 MVT::MvDot( X, *MX, oldDot );
1274 for (
int i=0; i<nq; i++) {
1277 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1284 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1287 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1293 OPT::Apply( *(this->_Op), X, *MX);
1297 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1298 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1300 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1303 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1310 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1313 MVT::MvDot( X, *MX, newDot );
1327 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1330 for (
int i=0; i<nq; i++) {
1335 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1343 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1346 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1352 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1356 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1358 else if (xc <= qcs[i]) {
1360 OPT::Apply( *(this->_Op), X, *MX);
1371 template<
class ScalarType,
class MV,
class OP>
1373 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X,
Teuchos::RCP<MV> MX,
1378 int xc = MVT::GetNumberVecs( X );
1379 bool dep_flg =
false;
1380 const ScalarType ONE = SCT::one();
1382 std::vector<int> qcs( nq );
1383 for (
int i=0; i<nq; i++) {
1384 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1390 std::vector<ScalarType> oldDot( xc );
1392 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1395 MVT::MvDot( X, *MX, oldDot );
1400 for (
int i=0; i<nq; i++) {
1403 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1410 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1413 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1419 OPT::Apply( *(this->_Op), X, *MX);
1423 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1424 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1426 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1429 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1436 for (
int j = 1; j < max_blk_ortho_; ++j) {
1438 for (
int i=0; i<nq; i++) {
1443 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1451 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1454 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1460 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1464 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1466 else if (xc <= qcs[i]) {
1468 OPT::Apply( *(this->_Op), X, *MX);
1475 std::vector<ScalarType> newDot(xc);
1477 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1480 MVT::MvDot( X, *MX, newDot );
1484 for (
int i=0; i<xc; i++){
1485 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1495 template<
class ScalarType,
class MV,
class OP>
1497 DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X,
Teuchos::RCP<MV> MX,
1504 const ScalarType ONE = SCT::one();
1505 const ScalarType ZERO = SCT::zero();
1508 int xc = MVT::GetNumberVecs( X );
1509 std::vector<int> indX( 1 );
1510 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1512 std::vector<int> qcs( nq );
1513 for (
int i=0; i<nq; i++) {
1514 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1523 for (
int j=0; j<xc; j++) {
1525 bool dep_flg =
false;
1529 std::vector<int> index( j );
1530 for (
int ind=0; ind<j; ind++) {
1533 lastQ = MVT::CloneView( X, index );
1536 Q.push_back( lastQ );
1538 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1543 Xj = MVT::CloneViewNonConst( X, indX );
1545 MXj = MVT::CloneViewNonConst( *MX, indX );
1553 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1556 MVT::MvDot( *Xj, *MXj, oldDot );
1561 for (
int i=0; i<Q.size(); i++) {
1568 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1575 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1578 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1584 OPT::Apply( *(this->_Op), *Xj, *MXj);
1588 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1589 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1591 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1594 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1602 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1605 MVT::MvDot( *Xj, *MXj, newDot );
1611 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
1613 for (
int i=0; i<Q.size(); i++) {
1619 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1626 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1629 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1635 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1639 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1641 else if (xc <= qcs[i]) {
1643 OPT::Apply( *(this->_Op), *Xj, *MXj);
1650 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1653 MVT::MvDot( *Xj, *MXj, newDot );
1658 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1664 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1666 MVT::MvScale( *Xj, ONE/diag );
1669 MVT::MvScale( *MXj, ONE/diag );
1679 MVT::MvRandom( *tempXj );
1681 tempMXj = MVT::Clone( X, 1 );
1682 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1688 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1691 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1694 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
1696 for (
int i=0; i<Q.size(); i++) {
1701 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1707 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1710 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1716 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1720 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1722 else if (xc <= qcs[i]) {
1724 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1733 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1736 MVT::MvDot( *tempXj, *tempMXj, newDot );
1740 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1741 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1747 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1749 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1769 template<
class ScalarType,
class MV,
class OP>
1773 using Teuchos::parameterList;
1776 RCP<ParameterList> params = parameterList (
"DGKS");
1781 "Maximum number of orthogonalization passes (includes the " 1782 "first). Default is 2, since \"twice is enough\" for Krylov " 1785 "Block reorthogonalization threshold.");
1787 "(Non-block) reorthogonalization threshold.");
1789 "Singular block detection threshold.");
1794 template<
class ScalarType,
class MV,
class OP>
1800 RCP<ParameterList> params = getDGKSDefaultParameters<ScalarType, MV, OP>();
1802 params->set (
"maxNumOrthogPasses",
1804 params->set (
"blkTol",
1806 params->set (
"depTol",
1808 params->set (
"singTol",
1816 #endif // BELOS_DGKS_ORTHOMANAGER_HPP void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
static const MagnitudeType dep_tol_fast_
(Non-block) reorthogonalization threshold (fast).
int normalize(MV &X, Teuchos::RCP< MV > MX, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method takes a multivector X and attempts to compute an orthonormal basis for ...
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
bool is_null(const boost::shared_ptr< T > &p)
Teuchos::RCP< Teuchos::ParameterList > getDGKSDefaultParameters()
"Default" parameters for robustness and accuracy.
DGKSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_blk_ortho=max_blk_ortho_default_, const MagnitudeType blk_tol=blk_tol_default_, const MagnitudeType dep_tol=dep_tol_default_, const MagnitudeType sing_tol=sing_tol_default_)
Constructor specifying re-orthogonalization tolerance.
~DGKSOrthoManager()
Destructor.
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
static const MagnitudeType dep_tol_default_
(Non-block) reorthogonalization threshold (default).
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
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::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
DGKSOrthoManager(const Teuchos::RCP< Teuchos::ParameterList > &plist, const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null)
Constructor that takes a list of parameters.
Class which defines basic traits for the operator type.
static const int max_blk_ortho_default_
Max number of (re)orthogonalization steps, including the first (default).
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
Traits class which defines basic operations on multivectors.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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 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)
MagnitudeType getDepTol() const
Return parameter for re-orthogonalization threshhold.
RCP< ParameterList > getNonconstParameterList()
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
int normalize(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method calls normalize(X,Teuchos::null,B); see documentation for that function.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
virtual int projectAndNormalizeWithMxImpl(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for ...
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
Class which defines basic traits for the operator type.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Provides the inner product defining the orthogonality concepts, using the provided operator...
Teuchos::RCP< Teuchos::ParameterList > getDGKSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
Belos header file which uses auto-configuration information to include necessary C++ headers...
MagnitudeType getSingTol() const
Return parameter for singular block detection.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogError(const MV &X1, const MV &X2) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
void project(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a list of (mutually and internally) orthonormal bases Q, this method takes a multivector X and ...