47#ifndef BELOS_DGKS_ORTHOMANAGER_HPP
48#define BELOS_DGKS_ORTHOMANAGER_HPP
64#include "Teuchos_as.hpp"
65#ifdef BELOS_TEUCHOS_TIME_MONITOR
66#include "Teuchos_TimeMonitor.hpp"
72 template<
class ScalarType,
class MV,
class OP>
76 template<
class ScalarType,
class MV,
class OP>
79 template<
class ScalarType,
class MV,
class OP>
84 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
85 typedef typename Teuchos::ScalarTraits<MagnitudeType>
MGT;
86 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
96 Teuchos::RCP<const OP> Op = Teuchos::null,
108#ifdef BELOS_TEUCHOS_TIME_MONITOR
109 std::stringstream ss;
112 std::string orthoLabel = ss.str() +
": Orthogonalization";
113 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
115 std::string updateLabel = ss.str() +
": Ortho (Update)";
116 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
118 std::string normLabel = ss.str() +
": Ortho (Norm)";
119 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
121 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
122 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
128 const std::string& label =
"Belos",
129 Teuchos::RCP<const OP> Op = Teuchos::null)
139#ifdef BELOS_TEUCHOS_TIME_MONITOR
140 std::stringstream ss;
143 std::string orthoLabel = ss.str() +
": Orthogonalization";
144 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
146 std::string updateLabel = ss.str() +
": Ortho (Update)";
147 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
149 std::string normLabel = ss.str() +
": Ortho (Norm)";
150 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
152 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
153 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
167 using Teuchos::ParameterList;
168 using Teuchos::parameterList;
172 RCP<ParameterList> params;
173 if (plist.is_null()) {
175 params = parameterList (*defaultParams);
178 params->validateParametersAndSetDefaults (*defaultParams);
186 const int maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
196 this->setMyParamList (params);
199 Teuchos::RCP<const Teuchos::ParameterList>
203 defaultParams_ = Belos::getDGKSDefaultParameters<ScalarType, MV, OP>();
217 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
218 if (! params.is_null()) {
223 params->set (
"blkTol", blk_tol);
231 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
232 if (! params.is_null()) {
233 params->set (
"depTol", dep_tol);
241 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
242 if (! params.is_null()) {
243 params->set (
"singTol", sing_tol);
290 void project ( MV &X, Teuchos::RCP<MV> MX,
291 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
292 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
298 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
299 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
329 int normalize ( MV &X, Teuchos::RCP<MV> MX,
330 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
335 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
399 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
400 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
401 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
413 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
424 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
430 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
439 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
440 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
449 void setLabel(
const std::string& label);
493#ifdef BELOS_TEUCHOS_TIME_MONITOR
494 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
501 int findBasis(MV &X, Teuchos::RCP<MV> MX,
502 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
503 bool completeBasis,
int howMany = -1 )
const;
506 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
507 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
508 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
511 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
512 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
513 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
529 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
530 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
531 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
535 template<
class ScalarType,
class MV,
class OP>
538 template<
class ScalarType,
class MV,
class OP>
541 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
544 template<
class ScalarType,
class MV,
class OP>
547 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::one()
548 / Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
551 template<
class ScalarType,
class MV,
class OP>
554 = 10*Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
556 template<
class ScalarType,
class MV,
class OP>
559 template<
class ScalarType,
class MV,
class OP>
562 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
564 template<
class ScalarType,
class MV,
class OP>
567 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
569 template<
class ScalarType,
class MV,
class OP>
572 = Teuchos::ScalarTraits<typename DGKSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
576 template<
class ScalarType,
class MV,
class OP>
579 if (label != label_) {
581#ifdef BELOS_TEUCHOS_TIME_MONITOR
582 std::stringstream ss;
583 ss << label_ +
": DGKS[" << max_blk_ortho_ <<
"]";
585 std::string orthoLabel = ss.str() +
": Orthogonalization";
586 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
588 std::string updateLabel = ss.str() +
": Ortho (Update)";
589 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
591 std::string normLabel = ss.str() +
": Ortho (Norm)";
592 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
594 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
595 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
602 template<
class ScalarType,
class MV,
class OP>
603 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
605 const ScalarType ONE = SCT::one();
606 int rank = MVT::GetNumberVecs(X);
607 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
609#ifdef BELOS_TEUCHOS_TIME_MONITOR
610 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
614 for (
int i=0; i<rank; i++) {
617 return xTx.normFrobenius();
622 template<
class ScalarType,
class MV,
class OP>
623 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
625 int r1 = MVT::GetNumberVecs(X1);
626 int r2 = MVT::GetNumberVecs(X2);
627 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
629#ifdef BELOS_TEUCHOS_TIME_MONITOR
630 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
634 return xTx.normFrobenius();
639 template<
class ScalarType,
class MV,
class OP>
644 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
645 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
646 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
648 using Teuchos::Array;
650 using Teuchos::is_null;
653 using Teuchos::SerialDenseMatrix;
654 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
655 typedef typename Array< RCP< const MV > >::size_type size_type;
657#ifdef BELOS_TEUCHOS_TIME_MONITOR
658 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
661 ScalarType ONE = SCT::one();
665 int xc = MVT::GetNumberVecs( X );
666 ptrdiff_t xr = MVT::GetGlobalLength( X );
673 B = rcp (
new serial_dense_matrix_type (xc, xc));
683 for (size_type k = 0; k < nq; ++k)
685 const int numRows = MVT::GetNumberVecs (*Q[k]);
686 const int numCols = xc;
689 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
690 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
692 int err = C[k]->reshape (numRows, numCols);
693 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
694 "DGKS orthogonalization: failed to reshape "
695 "C[" << k <<
"] (the array of block "
696 "coefficients resulting from projecting X "
697 "against Q[1:" << nq <<
"]).");
703 if (MX == Teuchos::null) {
705 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
706 OPT::Apply(*(this->_Op),X,*MX);
711 MX = Teuchos::rcp( &X,
false );
714 int mxc = MVT::GetNumberVecs( *MX );
715 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
718 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::DGKSOrthoManager::projectAndNormalize(): X must be non-empty" );
721 for (
int i=0; i<nq; i++) {
722 numbas += MVT::GetNumberVecs( *Q[i] );
726 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
727 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
729 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
730 "Belos::DGKSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
732 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
733 "Belos::DGKSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
739 bool dep_flg =
false;
745 dep_flg = blkOrtho1( X, MX, C, Q );
748 if ( B == Teuchos::null ) {
749 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
751 std::vector<ScalarType> diag(1);
753#ifdef BELOS_TEUCHOS_TIME_MONITOR
754 Teuchos::TimeMonitor normTimer( *timerNorm_ );
756 MVT::MvDot( X, *MX, diag );
758 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
760 if (SCT::magnitude((*B)(0,0)) > ZERO) {
762 MVT::MvScale( X, ONE/(*B)(0,0) );
765 MVT::MvScale( *MX, ONE/(*B)(0,0) );
772 Teuchos::RCP<MV> tmpX, tmpMX;
773 tmpX = MVT::CloneCopy(X);
775 tmpMX = MVT::CloneCopy(*MX);
779 dep_flg = blkOrtho( X, MX, C, Q );
784 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
787 MVT::Assign( *tmpX, X );
789 MVT::Assign( *tmpMX, *MX );
794 rank = findBasis( X, MX, B,
false );
798 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
801 MVT::Assign( *tmpX, X );
803 MVT::Assign( *tmpMX, *MX );
810 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
811 "Belos::DGKSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
821 template<
class ScalarType,
class MV,
class OP>
823 MV &X, Teuchos::RCP<MV> MX,
824 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
826#ifdef BELOS_TEUCHOS_TIME_MONITOR
827 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
831 return findBasis(X, MX, B,
true);
838 template<
class ScalarType,
class MV,
class OP>
840 MV &X, Teuchos::RCP<MV> MX,
841 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
842 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
858#ifdef BELOS_TEUCHOS_TIME_MONITOR
859 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
862 int xc = MVT::GetNumberVecs( X );
863 ptrdiff_t xr = MVT::GetGlobalLength( X );
865 std::vector<int> qcs(nq);
867 if (nq == 0 || xc == 0 || xr == 0) {
870 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
879 if (MX == Teuchos::null) {
881 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
882 OPT::Apply(*(this->_Op),X,*MX);
887 MX = Teuchos::rcp( &X,
false );
889 int mxc = MVT::GetNumberVecs( *MX );
890 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
893 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
894 "Belos::DGKSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
896 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
897 "Belos::DGKSOrthoManager::project(): Size of X not consistant with MX,Q" );
900 for (
int i=0; i<nq; i++) {
901 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
902 "Belos::DGKSOrthoManager::project(): Q lengths not mutually consistant" );
903 qcs[i] = MVT::GetNumberVecs( *Q[i] );
904 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
905 "Belos::DGKSOrthoManager::project(): Q has less rows than columns" );
908 if ( C[i] == Teuchos::null ) {
909 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
912 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
913 "Belos::DGKSOrthoManager::project(): Size of Q not consistant with size of C" );
918 blkOrtho( X, MX, C, Q );
925 template<
class ScalarType,
class MV,
class OP>
927 MV &X, Teuchos::RCP<MV> MX,
928 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
929 bool completeBasis,
int howMany )
const {
946 const ScalarType ONE = SCT::one();
949 int xc = MVT::GetNumberVecs( X );
950 ptrdiff_t xr = MVT::GetGlobalLength( X );
963 if (MX == Teuchos::null) {
965 MX = MVT::Clone(X,xc);
966 OPT::Apply(*(this->_Op),X,*MX);
973 if ( B == Teuchos::null ) {
974 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
977 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
978 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
981 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
982 "Belos::DGKSOrthoManager::findBasis(): X must be non-empty" );
983 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
984 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of B" );
985 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
986 "Belos::DGKSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
987 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<ptrdiff_t
>(xc) > xr, std::invalid_argument,
988 "Belos::DGKSOrthoManager::findBasis(): Size of X not feasible for normalization" );
989 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
990 "Belos::DGKSOrthoManager::findBasis(): Invalid howMany parameter" );
995 int xstart = xc - howMany;
997 for (
int j = xstart; j < xc; j++) {
1003 bool addVec =
false;
1006 std::vector<int> index(1);
1008 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1009 Teuchos::RCP<MV> MXj;
1010 if ((this->_hasOp)) {
1012 MXj = MVT::CloneViewNonConst( *MX, index );
1020 std::vector<int> prev_idx( numX );
1021 Teuchos::RCP<const MV> prevX, prevMX;
1022 Teuchos::RCP<MV> oldMXj;
1025 for (
int i=0; i<numX; i++) {
1028 prevX = MVT::CloneView( X, prev_idx );
1030 prevMX = MVT::CloneView( *MX, prev_idx );
1033 oldMXj = MVT::CloneCopy( *MXj );
1037 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1038 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1043#ifdef BELOS_TEUCHOS_TIME_MONITOR
1044 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1046 MVT::MvDot( *Xj, *MXj, oldDot );
1049 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO,
OrthoError,
1050 "Belos::DGKSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1057#ifdef BELOS_TEUCHOS_TIME_MONITOR
1058 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1065#ifdef BELOS_TEUCHOS_TIME_MONITOR
1066 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1068 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1076#ifdef BELOS_TEUCHOS_TIME_MONITOR
1077 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1079 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1084#ifdef BELOS_TEUCHOS_TIME_MONITOR
1085 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1087 MVT::MvDot( *Xj, *MXj, newDot );
1091 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1094 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1096#ifdef BELOS_TEUCHOS_TIME_MONITOR
1097 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1104#ifdef BELOS_TEUCHOS_TIME_MONITOR
1105 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1107 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1109 if ((this->_hasOp)) {
1110#ifdef BELOS_TEUCHOS_TIME_MONITOR
1111 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1113 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1121#ifdef BELOS_TEUCHOS_TIME_MONITOR
1122 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1124 MVT::MvDot( *Xj, *oldMXj, newDot );
1127 newDot[0] = oldDot[0];
1131 if (completeBasis) {
1135 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1140 std::cout <<
"Belos::DGKSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1143 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1144 Teuchos::RCP<MV> tempMXj;
1145 MVT::MvRandom( *tempXj );
1147 tempMXj = MVT::Clone( X, 1 );
1148 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1154#ifdef BELOS_TEUCHOS_TIME_MONITOR
1155 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1157 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1160 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++){
1162#ifdef BELOS_TEUCHOS_TIME_MONITOR
1163 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1168#ifdef BELOS_TEUCHOS_TIME_MONITOR
1169 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1171 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1174#ifdef BELOS_TEUCHOS_TIME_MONITOR
1175 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1177 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1182#ifdef BELOS_TEUCHOS_TIME_MONITOR
1183 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1185 MVT::MvDot( *tempXj, *tempMXj, newDot );
1188 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1190 MVT::Assign( *tempXj, *Xj );
1192 MVT::Assign( *tempMXj, *MXj );
1204 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1212 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1214 if (SCT::magnitude(diag) > ZERO) {
1215 MVT::MvScale( *Xj, ONE/diag );
1218 MVT::MvScale( *MXj, ONE/diag );
1232 for (
int i=0; i<numX; i++) {
1233 (*B)(i,j) = product(i,0);
1244 template<
class ScalarType,
class MV,
class OP>
1247 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1248 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
1251 int xc = MVT::GetNumberVecs( X );
1252 const ScalarType ONE = SCT::one();
1254 std::vector<int> qcs( nq );
1255 for (
int i=0; i<nq; i++) {
1256 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1260 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1262#ifdef BELOS_TEUCHOS_TIME_MONITOR
1263 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1265 MVT::MvDot( X, *MX, oldDot );
1268 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1270 for (
int i=0; i<nq; i++) {
1273#ifdef BELOS_TEUCHOS_TIME_MONITOR
1274 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1280#ifdef BELOS_TEUCHOS_TIME_MONITOR
1281 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1283 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1289 OPT::Apply( *(this->_Op), X, *MX);
1293 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1294 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1296#ifdef BELOS_TEUCHOS_TIME_MONITOR
1297 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1299 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1306#ifdef BELOS_TEUCHOS_TIME_MONITOR
1307 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1309 MVT::MvDot( X, *MX, newDot );
1323 if ( MGT::squareroot(SCT::magnitude(newDot[0])) < dep_tol_*MGT::squareroot(SCT::magnitude(oldDot[0])) ) {
1326 for (
int i=0; i<nq; i++) {
1327 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1331#ifdef BELOS_TEUCHOS_TIME_MONITOR
1332 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1339#ifdef BELOS_TEUCHOS_TIME_MONITOR
1340 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1342 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1348#ifdef BELOS_TEUCHOS_TIME_MONITOR
1349 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1352 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1354 else if (xc <= qcs[i]) {
1356 OPT::Apply( *(this->_Op), X, *MX);
1367 template<
class ScalarType,
class MV,
class OP>
1370 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1371 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
1374 int xc = MVT::GetNumberVecs( X );
1375 bool dep_flg =
false;
1376 const ScalarType ONE = SCT::one();
1378 std::vector<int> qcs( nq );
1379 for (
int i=0; i<nq; i++) {
1380 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1386 std::vector<ScalarType> oldDot( xc );
1388#ifdef BELOS_TEUCHOS_TIME_MONITOR
1389 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1391 MVT::MvDot( X, *MX, oldDot );
1394 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1396 for (
int i=0; i<nq; i++) {
1399#ifdef BELOS_TEUCHOS_TIME_MONITOR
1400 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1406#ifdef BELOS_TEUCHOS_TIME_MONITOR
1407 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1409 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1415 OPT::Apply( *(this->_Op), X, *MX);
1419 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1420 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1422#ifdef BELOS_TEUCHOS_TIME_MONITOR
1423 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1425 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1432 for (
int j = 1; j < max_blk_ortho_; ++j) {
1434 for (
int i=0; i<nq; i++) {
1435 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
1439#ifdef BELOS_TEUCHOS_TIME_MONITOR
1440 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1447#ifdef BELOS_TEUCHOS_TIME_MONITOR
1448 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1450 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1456#ifdef BELOS_TEUCHOS_TIME_MONITOR
1457 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1460 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1462 else if (xc <= qcs[i]) {
1464 OPT::Apply( *(this->_Op), X, *MX);
1471 std::vector<ScalarType> newDot(xc);
1473#ifdef BELOS_TEUCHOS_TIME_MONITOR
1474 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1476 MVT::MvDot( X, *MX, newDot );
1480 for (
int i=0; i<xc; i++){
1481 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1491 template<
class ScalarType,
class MV,
class OP>
1494 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1495 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1496 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const
1498 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1500 const ScalarType ONE = SCT::one();
1501 const ScalarType ZERO = SCT::zero();
1504 int xc = MVT::GetNumberVecs( X );
1505 std::vector<int> indX( 1 );
1506 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1508 std::vector<int> qcs( nq );
1509 for (
int i=0; i<nq; i++) {
1510 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1514 Teuchos::RCP<const MV> lastQ;
1515 Teuchos::RCP<MV> Xj, MXj;
1516 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1519 for (
int j=0; j<xc; j++) {
1521 bool dep_flg =
false;
1525 std::vector<int> index( j );
1526 for (
int ind=0; ind<j; ind++) {
1529 lastQ = MVT::CloneView( X, index );
1532 Q.push_back( lastQ );
1534 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1539 Xj = MVT::CloneViewNonConst( X, indX );
1541 MXj = MVT::CloneViewNonConst( *MX, indX );
1549#ifdef BELOS_TEUCHOS_TIME_MONITOR
1550 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1552 MVT::MvDot( *Xj, *MXj, oldDot );
1555 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1557 for (
int i=0; i<Q.size(); i++) {
1560 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1564#ifdef BELOS_TEUCHOS_TIME_MONITOR
1565 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1571#ifdef BELOS_TEUCHOS_TIME_MONITOR
1572 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1574 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1580 OPT::Apply( *(this->_Op), *Xj, *MXj);
1584 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1585 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1587#ifdef BELOS_TEUCHOS_TIME_MONITOR
1588 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1590 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1598#ifdef BELOS_TEUCHOS_TIME_MONITOR
1599 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1601 MVT::MvDot( *Xj, *MXj, newDot );
1607 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*dep_tol_) ) {
1609 for (
int i=0; i<Q.size(); i++) {
1610 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1611 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1615#ifdef BELOS_TEUCHOS_TIME_MONITOR
1616 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1622#ifdef BELOS_TEUCHOS_TIME_MONITOR
1623 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1625 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1631#ifdef BELOS_TEUCHOS_TIME_MONITOR
1632 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1635 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1637 else if (xc <= qcs[i]) {
1639 OPT::Apply( *(this->_Op), *Xj, *MXj);
1646#ifdef BELOS_TEUCHOS_TIME_MONITOR
1647 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1649 MVT::MvDot( *Xj, *MXj, newDot );
1654 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1660 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1662 MVT::MvScale( *Xj, ONE/diag );
1665 MVT::MvScale( *MXj, ONE/diag );
1673 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1674 Teuchos::RCP<MV> tempMXj;
1675 MVT::MvRandom( *tempXj );
1677 tempMXj = MVT::Clone( X, 1 );
1678 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1684#ifdef BELOS_TEUCHOS_TIME_MONITOR
1685 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1687 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1690 for (
int num_orth=0; num_orth<max_blk_ortho_; num_orth++) {
1692 for (
int i=0; i<Q.size(); i++) {
1693 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1697#ifdef BELOS_TEUCHOS_TIME_MONITOR
1698 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1703#ifdef BELOS_TEUCHOS_TIME_MONITOR
1704 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1706 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1712#ifdef BELOS_TEUCHOS_TIME_MONITOR
1713 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1716 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1718 else if (xc <= qcs[i]) {
1720 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1729#ifdef BELOS_TEUCHOS_TIME_MONITOR
1730 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1732 MVT::MvDot( *tempXj, *tempMXj, newDot );
1736 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1737 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1743 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1745 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1765 template<
class ScalarType,
class MV,
class OP>
1768 using Teuchos::ParameterList;
1769 using Teuchos::parameterList;
1772 RCP<ParameterList> params = parameterList (
"DGKS");
1777 "Maximum number of orthogonalization passes (includes the "
1778 "first). Default is 2, since \"twice is enough\" for Krylov "
1781 "Block reorthogonalization threshold.");
1783 "(Non-block) reorthogonalization threshold.");
1785 "Singular block detection threshold.");
1790 template<
class ScalarType,
class MV,
class OP>
1793 using Teuchos::ParameterList;
1796 RCP<ParameterList> params = getDGKSDefaultParameters<ScalarType, MV, OP>();
1798 params->set (
"maxNumOrthogPasses",
1800 params->set (
"blkTol",
1802 params->set (
"depTol",
1804 params->set (
"singTol",
Belos header file which uses auto-configuration information to include necessary C++ headers.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
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.
MagnitudeType blk_tol_
Block reorthogonalization threshold.
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 ...
static const int max_blk_ortho_default_
Max number of (re)orthogonalization steps, including the first (default).
std::string label_
Label for timer(s).
~DGKSOrthoManager()
Destructor.
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
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.
static const MagnitudeType dep_tol_default_
(Non-block) reorthogonalization threshold (default).
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.
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
static const int max_blk_ortho_fast_
Max number of (re)orthogonalization steps, including the first (fast).
OperatorTraits< ScalarType, MV, OP > OPT
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
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.
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector.
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
MultiVecTraits< ScalarType, MV > MVT
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
int max_blk_ortho_
Max number of (re)orthogonalization steps, including the first.
MagnitudeType dep_tol_
(Non-block) reorthogonalization threshold.
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 .
MagnitudeType getDepTol() const
Return parameter for re-orthogonalization threshhold.
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
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...
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.
bool blkOrtho1(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.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
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.
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.
Teuchos::ScalarTraits< MagnitudeType > MGT
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 ,...
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
MagnitudeType sing_tol_
Singular block detection threshold.
Teuchos::ScalarTraits< ScalarType > SCT
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
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.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Exception thrown to signal error in an orthogonalization manager method.
Teuchos::RCP< Teuchos::ParameterList > getDGKSDefaultParameters()
"Default" parameters for robustness and accuracy.
Teuchos::RCP< Teuchos::ParameterList > getDGKSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.