47 #ifndef BELOS_DGKS_ORTHOMANAGER_HPP 48 #define BELOS_DGKS_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_blk_ortho = 2,
92 const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
93 const MagnitudeType dep_tol = MGT::one()/MGT::squareroot( 2*MGT::one() ),
102 #ifdef BELOS_TEUCHOS_TIME_MONITOR 103 std::string orthoLabel =
label_ +
": Orthogonalization";
104 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter( orthoLabel );
110 const std::string& label =
"Belos",
111 Teuchos::RCP<const OP> Op = Teuchos::null)
121 #ifdef BELOS_TEUCHOS_TIME_MONITOR 122 std::string orthoLabel =
label_ +
": Orthogonalization";
123 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter( orthoLabel );
137 using Teuchos::ParameterList;
138 using Teuchos::parameterList;
142 RCP<ParameterList> params;
143 if (plist.is_null()) {
145 params = parameterList (*defaultParams);
148 params->validateParametersAndSetDefaults (*defaultParams);
156 const int maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
166 setMyParamList (params);
169 Teuchos::RCP<const Teuchos::ParameterList>
173 using Teuchos::ParameterList;
174 using Teuchos::parameterList;
178 RCP<ParameterList> params = parameterList (
"DGKS");
183 const int defaultMaxNumOrthogPasses = 2;
185 as<MagnitudeType> (10) * MGT::squareroot (eps);
187 MGT::one() / MGT::squareroot (as<MagnitudeType> (2));
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 " 208 Teuchos::RCP<const Teuchos::ParameterList>
211 using Teuchos::ParameterList;
217 RCP<ParameterList> params = rcp (
new ParameterList (*defaultParams));
219 const int maxBlkOrtho = 1;
224 params->set (
"maxNumOrthogPasses", maxBlkOrtho);
225 params->set (
"blkTol", blkTol);
226 params->set (
"depTol", depTol);
227 params->set (
"singTol", singTol);
238 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
239 if (! params.is_null()) {
244 params->set (
"blkTol", blk_tol);
252 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
253 if (! params.is_null()) {
254 params->set (
"depTol", dep_tol);
262 Teuchos::RCP<Teuchos::ParameterList> params = getNonconstParameterList();
263 if (! params.is_null()) {
264 params->set (
"singTol", sing_tol);
311 void project ( MV &X, Teuchos::RCP<MV> MX,
312 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
313 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
319 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
320 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
350 int normalize ( MV &X, Teuchos::RCP<MV> MX,
351 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
356 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
420 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
421 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
422 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
434 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
445 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
451 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
460 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
461 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
470 void setLabel(
const std::string& label);
491 #ifdef BELOS_TEUCHOS_TIME_MONITOR 492 Teuchos::RCP<Teuchos::Time> timerOrtho_;
493 #endif // BELOS_TEUCHOS_TIME_MONITOR 499 int findBasis(MV &X, Teuchos::RCP<MV> MX,
500 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
501 bool completeBasis,
int howMany = -1 )
const;
504 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
505 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
506 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
522 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
523 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
524 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
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 536 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
543 template<
class ScalarType,
class MV,
class OP>
544 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
546 const ScalarType ONE = SCT::one();
547 int rank = MVT::GetNumberVecs(X);
548 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
550 for (
int i=0; i<rank; i++) {
553 return xTx.normFrobenius();
558 template<
class ScalarType,
class MV,
class OP>
559 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
561 int r1 = MVT::GetNumberVecs(X1);
562 int r2 = MVT::GetNumberVecs(X2);
563 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
565 return xTx.normFrobenius();
570 template<
class ScalarType,
class MV,
class OP>
575 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
576 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
577 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 579 using Teuchos::Array;
581 using Teuchos::is_null;
584 using Teuchos::SerialDenseMatrix;
585 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
586 typedef typename Array< RCP< const MV > >::size_type size_type;
588 #ifdef BELOS_TEUCHOS_TIME_MONITOR 589 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
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);
624 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
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);
642 MX = Teuchos::rcp( &X,
false );
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] );
657 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
658 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
660 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
661 "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
663 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
664 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
670 bool dep_flg =
false;
673 Teuchos::RCP<MV> tmpX, tmpMX;
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 );
710 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
711 "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
721 template<
class ScalarType,
class MV,
class OP>
723 MV &X, Teuchos::RCP<MV> MX,
724 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
726 #ifdef BELOS_TEUCHOS_TIME_MONITOR 727 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
731 return findBasis(X, MX, B,
true);
737 template<
class ScalarType,
class MV,
class OP>
739 MV &X, Teuchos::RCP<MV> MX,
740 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
741 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
757 #ifdef BELOS_TEUCHOS_TIME_MONITOR 758 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
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);
786 MX = Teuchos::rcp( &X,
false );
788 int mxc = MVT::GetNumberVecs( *MX );
789 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
792 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
793 "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
795 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
796 "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
800 for (
int i=0; i<nq; i++) {
801 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
802 "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
803 qcs[i] = MVT::GetNumberVecs( *Q[i] );
804 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
805 "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
809 if ( C[i] == Teuchos::null ) {
810 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
813 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
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>
828 MV &X, Teuchos::RCP<MV> MX,
829 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
830 bool completeBasis,
int howMany )
const {
847 const ScalarType ONE = SCT::one();
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 ) {
875 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
878 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
879 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
882 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
883 "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
884 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
885 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
886 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
887 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
888 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
889 "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
890 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
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);
909 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
910 Teuchos::RCP<MV> MXj;
911 if ((this->_hasOp)) {
913 MXj = MVT::CloneViewNonConst( *MX, index );
921 std::vector<int> prev_idx( numX );
922 Teuchos::RCP<const MV> prevX, prevMX;
925 for (
int i=0; i<numX; i++) {
928 prevX = MVT::CloneView( X, prev_idx );
930 prevMX = MVT::CloneView( *MX, prev_idx );
935 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
936 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
940 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
941 MVT::MvDot( *Xj, *MXj, oldDot );
943 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO,
OrthoError,
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]) ) {
971 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
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;
999 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1000 Teuchos::RCP<MV> tempMXj;
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>
1080 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1081 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 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 );
1099 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
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++) {
1125 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
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>
1165 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1166 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1167 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1169 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
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] );
1185 Teuchos::RCP<const MV> lastQ;
1186 Teuchos::RCP<MV> Xj, MXj;
1187 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
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 );
1221 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1223 for (
int i=0; i<Q.size(); i++) {
1226 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
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++) {
1256 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1257 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
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 );
1302 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1303 Teuchos::RCP<MV> tempMXj;
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++) {
1317 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
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.
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.
OperatorTraits< ScalarType, MV, OP > OPT
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
MagnitudeType sing_tol_
Singular block detection threshold.
Teuchos::ScalarTraits< ScalarType > SCT
~DGKSOrthoManager()
Destructor.
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
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.
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.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Traits class which defines basic operations on multivectors.
MagnitudeType dep_tol_
(Non-block) reorthogonalization threshold.
Teuchos::ScalarTraits< MagnitudeType > MGT
std::string label_
Label for timer(s).
Exception thrown to signal error in an orthogonalization manager method.
MagnitudeType blk_tol_
Block reorthogonalization threshold.
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.
MagnitudeType getDepTol() const
Return parameter for re-orthogonalization threshhold.
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.
int blkOrthoSing(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 > > QQ) const
Project X against QQ and normalize X, one vector at a time.
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...
int max_blk_ortho_
Max number of (re)orthogonalization steps, including the first.
Class which defines basic traits for the operator type.
int findBasis(MV &X, Teuchos::RCP< MV > MX, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > C, bool completeBasis, int howMany=-1) const
Routine to find an orthonormal basis for X.
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...
bool blkOrtho(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Routine to compute the block orthogonalization.
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 ...
MultiVecTraits< ScalarType, MV > MVT
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 ...