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 72 template<
class ScalarType,
class MV,
class OP>
91 const int max_blk_ortho = 2,
94 const MagnitudeType sing_tol = 10*
MGT::eps() )
96 max_blk_ortho_( max_blk_ortho ),
99 sing_tol_( sing_tol ),
102 #ifdef BELOS_TEUCHOS_TIME_MONITOR 103 std::string orthoLabel = label_ +
": Orthogonalization";
110 const std::string& label =
"Belos",
114 blk_tol_ (
Teuchos::
as<MagnitudeType>(10) *
MGT::squareroot(
MGT::eps())),
115 dep_tol_ (
MGT::one() /
MGT::squareroot (
Teuchos::
as<MagnitudeType>(2))),
116 sing_tol_ (
Teuchos::
as<MagnitudeType>(10) *
MGT::eps()),
121 #ifdef BELOS_TEUCHOS_TIME_MONITOR 122 std::string orthoLabel = label_ +
": Orthogonalization";
138 using Teuchos::parameterList;
142 RCP<ParameterList> params;
145 params = parameterList (*defaultParams);
156 const int maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
157 const MagnitudeType blkTol = params->get<MagnitudeType> (
"blkTol");
158 const MagnitudeType depTol = params->get<MagnitudeType> (
"depTol");
159 const MagnitudeType singTol = params->get<MagnitudeType> (
"singTol");
161 max_blk_ortho_ = maxNumOrthogPasses;
174 using Teuchos::parameterList;
177 if (defaultParams_.
is_null()) {
178 RCP<ParameterList> params = parameterList (
"DGKS");
179 const MagnitudeType eps =
MGT::eps ();
183 const int defaultMaxNumOrthogPasses = 2;
184 const MagnitudeType defaultBlkTol =
186 const MagnitudeType defaultDepTol =
189 const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
191 params->set (
"maxNumOrthogPasses", defaultMaxNumOrthogPasses,
192 "Maximum number of orthogonalization passes (includes the " 193 "first). Default is 2, since \"twice is enough\" for Krylov " 195 params->set (
"blkTol", defaultBlkTol,
"Block reorthogonalization " 197 params->set (
"depTol", defaultDepTol,
198 "(Non-block) reorthogonalization threshold.");
199 params->set (
"singTol", defaultSingTol,
"Singular block detection " 201 defaultParams_ = params;
203 return defaultParams_;
217 RCP<ParameterList> params =
rcp (
new ParameterList (*defaultParams));
219 const int maxBlkOrtho = 1;
220 const MagnitudeType blkTol =
MGT::zero();
221 const MagnitudeType depTol =
MGT::zero();
222 const MagnitudeType singTol =
MGT::zero();
224 params->set (
"maxNumOrthogPasses", maxBlkOrtho);
225 params->set (
"blkTol", blkTol);
226 params->set (
"depTol", depTol);
227 params->set (
"singTol", singTol);
244 params->
set (
"blkTol", blk_tol);
254 params->
set (
"depTol", dep_tol);
264 params->
set (
"singTol", sing_tol);
266 sing_tol_ = sing_tol;
470 void setLabel(
const std::string& label);
474 const std::string&
getLabel()
const {
return label_; }
483 MagnitudeType blk_tol_;
485 MagnitudeType dep_tol_;
487 MagnitudeType sing_tol_;
491 #ifdef BELOS_TEUCHOS_TIME_MONITOR 493 #endif // BELOS_TEUCHOS_TIME_MONITOR 501 bool completeBasis,
int howMany = -1 )
const;
529 template<
class ScalarType,
class MV,
class OP>
532 if (label != label_) {
534 std::string orthoLabel = label_ +
": Orthogonalization";
535 #ifdef BELOS_TEUCHOS_TIME_MONITOR 543 template<
class ScalarType,
class MV,
class OP>
546 const ScalarType ONE = SCT::one();
547 int rank = MVT::GetNumberVecs(X);
550 for (
int i=0; i<rank; i++) {
558 template<
class ScalarType,
class MV,
class OP>
561 int r1 = MVT::GetNumberVecs(X1);
562 int r2 = MVT::GetNumberVecs(X2);
570 template<
class ScalarType,
class MV,
class OP>
586 typedef typename Array< RCP< const MV > >::size_type size_type;
588 #ifdef BELOS_TEUCHOS_TIME_MONITOR 592 ScalarType ONE = SCT::one();
593 ScalarType ZERO = SCT::zero();
596 int xc = MVT::GetNumberVecs( X );
597 ptrdiff_t xr = MVT::GetGlobalLength( X );
604 B =
rcp (
new serial_dense_matrix_type (xc, xc));
614 for (size_type k = 0; k < nq; ++k)
616 const int numRows = MVT::GetNumberVecs (*Q[k]);
617 const int numCols = xc;
620 C[k] =
rcp (
new serial_dense_matrix_type (numRows, numCols));
621 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
623 int err = C[k]->reshape (numRows, numCols);
625 "DGKS orthogonalization: failed to reshape " 626 "C[" << k <<
"] (the array of block " 627 "coefficients resulting from projecting X " 628 "against Q[1:" << nq <<
"]).");
634 if (MX == Teuchos::null) {
636 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
637 OPT::Apply(*(this->_Op),X,*MX);
645 int mxc = MVT::GetNumberVecs( *MX );
646 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
649 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
652 for (
int i=0; i<nq; i++) {
653 numbas += MVT::GetNumberVecs( *Q[i] );
658 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
661 "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
664 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
670 bool dep_flg =
false;
674 tmpX = MVT::CloneCopy(X);
676 tmpMX = MVT::CloneCopy(*MX);
680 dep_flg = blkOrtho( X, MX, C, Q );
685 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
688 MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
690 MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
695 rank = findBasis( X, MX, B,
false );
699 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
702 MVT::MvAddMv( ONE, *tmpX, ZERO, *tmpX, X );
704 MVT::MvAddMv( ONE, *tmpMX, ZERO, *tmpMX, *MX );
711 "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
721 template<
class ScalarType,
class MV,
class OP>
726 #ifdef BELOS_TEUCHOS_TIME_MONITOR 731 return findBasis(X, MX, B,
true);
737 template<
class ScalarType,
class MV,
class OP>
757 #ifdef BELOS_TEUCHOS_TIME_MONITOR 761 int xc = MVT::GetNumberVecs( X );
762 ptrdiff_t xr = MVT::GetGlobalLength( X );
764 std::vector<int> qcs(nq);
766 if (nq == 0 || xc == 0 || xr == 0) {
769 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
778 if (MX == Teuchos::null) {
780 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
781 OPT::Apply(*(this->_Op),X,*MX);
788 int mxc = MVT::GetNumberVecs( *MX );
789 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
793 "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
796 "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
800 for (
int i=0; i<nq; i++) {
802 "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
803 qcs[i] = MVT::GetNumberVecs( *Q[i] );
805 "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
809 if ( C[i] == Teuchos::null ) {
814 "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
819 blkOrtho( X, MX, C, Q );
826 template<
class ScalarType,
class MV,
class OP>
830 bool completeBasis,
int howMany )
const {
847 const ScalarType ONE = SCT::one();
848 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
850 int xc = MVT::GetNumberVecs( X );
851 ptrdiff_t xr = MVT::GetGlobalLength( X );
864 if (MX == Teuchos::null) {
866 MX = MVT::Clone(X,xc);
867 OPT::Apply(*(this->_Op),X,*MX);
874 if ( B == Teuchos::null ) {
878 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
879 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
883 "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
885 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
887 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
889 "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
891 "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
896 int xstart = xc - howMany;
898 for (
int j = xstart; j < xc; j++) {
907 std::vector<int> index(1);
911 if ((this->_hasOp)) {
913 MXj = MVT::CloneViewNonConst( *MX, index );
921 std::vector<int> prev_idx( numX );
925 for (
int i=0; i<numX; i++) {
928 prevX = MVT::CloneView( X, prev_idx );
930 prevMX = MVT::CloneView( *MX, prev_idx );
936 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
941 MVT::MvDot( *Xj, *MXj, oldDot );
944 "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
954 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
961 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
965 MVT::MvDot( *Xj, *MXj, newDot );
968 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(dep_tol_*oldDot[0]) ) {
975 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
976 if ((this->_hasOp)) {
977 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
984 MVT::MvDot( *Xj, *oldMXj, newDot );
991 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
996 std::cout <<
"Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1001 MVT::MvRandom( *tempXj );
1003 tempMXj = MVT::Clone( X, 1 );
1004 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1009 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1011 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
1013 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1015 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1019 MVT::MvDot( *tempXj, *tempMXj, newDot );
1021 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1023 MVT::MvAddMv( ONE, *tempXj, ZERO, *tempXj, *Xj );
1025 MVT::MvAddMv( ONE, *tempMXj, ZERO, *tempMXj, *MXj );
1037 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1045 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1047 if (SCT::magnitude(diag) > ZERO) {
1048 MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
1051 MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
1065 for (
int i=0; i<numX; i++) {
1066 (*B)(i,j) = product(i,0);
1077 template<
class ScalarType,
class MV,
class OP>
1079 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X,
Teuchos::RCP<MV> MX,
1084 int xc = MVT::GetNumberVecs( X );
1085 bool dep_flg =
false;
1086 const ScalarType ONE = SCT::one();
1088 std::vector<int> qcs( nq );
1089 for (
int i=0; i<nq; i++) {
1090 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1096 std::vector<ScalarType> oldDot( xc );
1097 MVT::MvDot( X, *MX, oldDot );
1101 for (
int i=0; i<nq; i++) {
1105 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1110 OPT::Apply( *(this->_Op), X, *MX);
1114 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1115 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1116 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1122 for (
int j = 1; j < max_blk_ortho_; ++j) {
1124 for (
int i=0; i<nq; i++) {
1130 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1136 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1138 else if (xc <= qcs[i]) {
1140 OPT::Apply( *(this->_Op), X, *MX);
1147 std::vector<ScalarType> newDot(xc);
1148 MVT::MvDot( X, *MX, newDot );
1151 for (
int i=0; i<xc; i++){
1152 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1162 template<
class ScalarType,
class MV,
class OP>
1164 DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X,
Teuchos::RCP<MV> MX,
1171 const ScalarType ONE = SCT::one();
1172 const ScalarType ZERO = SCT::zero();
1175 int xc = MVT::GetNumberVecs( X );
1176 std::vector<int> indX( 1 );
1177 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1179 std::vector<int> qcs( nq );
1180 for (
int i=0; i<nq; i++) {
1181 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1190 for (
int j=0; j<xc; j++) {
1192 bool dep_flg =
false;
1196 std::vector<int> index( j );
1197 for (
int ind=0; ind<j; ind++) {
1200 lastQ = MVT::CloneView( X, index );
1203 Q.push_back( lastQ );
1205 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1210 Xj = MVT::CloneViewNonConst( X, indX );
1212 MXj = MVT::CloneViewNonConst( *MX, indX );
1219 MVT::MvDot( *Xj, *MXj, oldDot );
1223 for (
int i=0; i<Q.size(); i++) {
1231 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1236 OPT::Apply( *(this->_Op), *Xj, *MXj);
1240 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1241 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1242 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1248 MVT::MvDot( *Xj, *MXj, newDot );
1253 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
1255 for (
int i=0; i<Q.size(); i++) {
1262 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1268 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1270 else if (xc <= qcs[i]) {
1272 OPT::Apply( *(this->_Op), *Xj, *MXj);
1278 MVT::MvDot( *Xj, *MXj, newDot );
1283 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1289 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1291 MVT::MvAddMv( ONE/diag, *Xj, ZERO, *Xj, *Xj );
1294 MVT::MvAddMv( ONE/diag, *MXj, ZERO, *MXj, *MXj );
1304 MVT::MvRandom( *tempXj );
1306 tempMXj = MVT::Clone( X, 1 );
1307 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1312 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1314 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
1316 for (
int i=0; i<Q.size(); i++) {
1321 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1327 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1329 else if (xc <= qcs[i]) {
1331 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1339 MVT::MvDot( *tempXj, *tempMXj, newDot );
1342 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1343 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1349 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1351 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1373 #endif // BELOS_DGKS_ORTHOMANAGER_HPP void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() 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 ...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
bool is_null(const boost::shared_ptr< T > &p)
static magnitudeType eps()
~DGKSOrthoManager()
Destructor.
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
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.
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)
void validateParametersAndSetDefaults(ParameterList const &validParamList, int const depth=1000)
MagnitudeType getDepTol() const
Return parameter for re-orthogonalization threshhold.
RCP< ParameterList > getNonconstParameterList()
DGKSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_blk_ortho=2, const MagnitudeType blk_tol=10 *MGT::squareroot(MGT::eps()), const MagnitudeType dep_tol=MGT::one()/MGT::squareroot(2 *MGT::one()), const MagnitudeType sing_tol=10 *MGT::eps())
Constructor specifying re-orthogonalization tolerance.
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 ...
TypeTo as(const TypeFrom &t)
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...
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 ...