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() ),
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";
104 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter( orthoLabel );
110 const std::string& label =
"Belos",
111 Teuchos::RCP<const OP> Op = Teuchos::null)
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";
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");
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;
166 setMyParamList (params);
169 Teuchos::RCP<const Teuchos::ParameterList>
173 using Teuchos::ParameterList;
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 =
185 as<MagnitudeType> (10) * MGT::squareroot (eps);
186 const MagnitudeType defaultDepTol =
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 " 201 defaultParams_ = params;
203 return defaultParams_;
208 Teuchos::RCP<const Teuchos::ParameterList>
211 using Teuchos::ParameterList;
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);
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);
266 sing_tol_ = 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);
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 492 Teuchos::RCP<Teuchos::Time> timerOrtho_;
493 #endif // BELOS_TEUCHOS_TIME_MONITOR 496 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
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;
521 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
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();
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 ) {
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>
1079 DGKSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
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>
1164 DGKSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
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.
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
~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.
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.
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.
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...
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 ...