47 #ifndef BELOS_IMGS_ORTHOMANAGER_HPP 48 #define BELOS_IMGS_ORTHOMANAGER_HPP 63 #include "Teuchos_SerialDenseMatrix.hpp" 64 #include "Teuchos_SerialDenseVector.hpp" 66 #include "Teuchos_as.hpp" 67 #ifdef BELOS_TEUCHOS_TIME_MONITOR 68 #include "Teuchos_TimeMonitor.hpp" 69 #endif // BELOS_TEUCHOS_TIME_MONITOR 74 template<
class ScalarType,
class MV,
class OP>
78 template<
class ScalarType,
class MV,
class OP>
81 template<
class ScalarType,
class MV,
class OP>
86 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
87 typedef typename Teuchos::ScalarTraits<MagnitudeType>
MGT;
88 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
98 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) :
138 #ifdef BELOS_TEUCHOS_TIME_MONITOR 139 std::stringstream ss;
142 std::string orthoLabel = ss.str() +
": Orthogonalization";
143 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
145 std::string updateLabel = ss.str() +
": Ortho (Update)";
146 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
148 std::string normLabel = ss.str() +
": Ortho (Norm)";
149 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
151 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
152 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
165 using Teuchos::Exceptions::InvalidParameterName;
166 using Teuchos::ParameterList;
167 using Teuchos::parameterList;
171 RCP<ParameterList> params;
172 if (plist.is_null()) {
173 params = parameterList (*defaultParams);
188 int maxNumOrthogPasses;
193 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
194 }
catch (InvalidParameterName&) {
195 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
196 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
208 }
catch (InvalidParameterName&) {
213 params->remove (
"depTol");
214 }
catch (InvalidParameterName&) {
217 params->set (
"blkTol", blkTol);
222 }
catch (InvalidParameterName&) {
224 params->set (
"singTol", singTol);
231 this->setMyParamList (params);
234 Teuchos::RCP<const Teuchos::ParameterList>
238 defaultParams_ = Belos::getIMGSDefaultParameters<ScalarType, MV, OP>();
294 void project ( MV &X, Teuchos::RCP<MV> MX,
295 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
296 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
302 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
303 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
333 int normalize ( MV &X, Teuchos::RCP<MV> MX,
334 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
339 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
388 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
389 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
390 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
400 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
409 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
415 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
424 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
425 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
434 void setLabel(
const std::string& label);
472 #ifdef BELOS_TEUCHOS_TIME_MONITOR 473 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
474 #endif // BELOS_TEUCHOS_TIME_MONITOR 480 int findBasis(MV &X, Teuchos::RCP<MV> MX,
481 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
482 bool completeBasis,
int howMany = -1 )
const;
485 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
486 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
487 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
490 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
491 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
492 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
508 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
509 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
510 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
514 template<
class ScalarType,
class MV,
class OP>
517 template<
class ScalarType,
class MV,
class OP>
520 = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
523 template<
class ScalarType,
class MV,
class OP>
526 = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
528 template<
class ScalarType,
class MV,
class OP>
531 template<
class ScalarType,
class MV,
class OP>
534 = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
536 template<
class ScalarType,
class MV,
class OP>
539 = Teuchos::ScalarTraits<typename IMGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
543 template<
class ScalarType,
class MV,
class OP>
546 if (label != label_) {
548 #ifdef BELOS_TEUCHOS_TIME_MONITOR 549 std::stringstream ss;
550 ss << label_ +
": IMGS[" << max_ortho_steps_ <<
"]";
552 std::string orthoLabel = ss.str() +
": Orthogonalization";
553 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
555 std::string updateLabel = ss.str() +
": Ortho (Update)";
556 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
558 std::string normLabel = ss.str() +
": Ortho (Norm)";
559 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
561 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
562 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
569 template<
class ScalarType,
class MV,
class OP>
570 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
572 const ScalarType ONE = SCT::one();
573 int rank = MVT::GetNumberVecs(X);
574 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
576 for (
int i=0; i<rank; i++) {
579 return xTx.normFrobenius();
584 template<
class ScalarType,
class MV,
class OP>
585 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
587 int r1 = MVT::GetNumberVecs(X1);
588 int r2 = MVT::GetNumberVecs(X2);
589 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
591 return xTx.normFrobenius();
596 template<
class ScalarType,
class MV,
class OP>
601 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
602 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
603 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 605 using Teuchos::Array;
607 using Teuchos::is_null;
610 using Teuchos::SerialDenseMatrix;
611 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
612 typedef typename Array< RCP< const MV > >::size_type size_type;
614 #ifdef BELOS_TEUCHOS_TIME_MONITOR 615 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
618 ScalarType ONE = SCT::one();
622 int xc = MVT::GetNumberVecs( X );
623 ptrdiff_t xr = MVT::GetGlobalLength( X );
630 B = rcp (
new serial_dense_matrix_type (xc, xc));
640 for (size_type k = 0; k < nq; ++k)
642 const int numRows = MVT::GetNumberVecs (*Q[k]);
643 const int numCols = xc;
646 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
647 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
649 int err = C[k]->reshape (numRows, numCols);
650 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
651 "IMGS orthogonalization: failed to reshape " 652 "C[" << k <<
"] (the array of block " 653 "coefficients resulting from projecting X " 654 "against Q[1:" << nq <<
"]).");
660 if (MX == Teuchos::null) {
662 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
663 OPT::Apply(*(this->_Op),X,*MX);
668 MX = Teuchos::rcp( &X,
false );
671 int mxc = MVT::GetNumberVecs( *MX );
672 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
675 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
678 for (
int i=0; i<nq; i++) {
679 numbas += MVT::GetNumberVecs( *Q[i] );
683 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
684 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
686 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
687 "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
689 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
690 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
696 bool dep_flg =
false;
699 Teuchos::RCP<MV> tmpX, tmpMX;
700 tmpX = MVT::CloneCopy(X);
702 tmpMX = MVT::CloneCopy(*MX);
709 dep_flg = blkOrtho1( X, MX, C, Q );
712 if ( B == Teuchos::null ) {
713 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
715 std::vector<ScalarType> diag(xc);
717 #ifdef BELOS_TEUCHOS_TIME_MONITOR 718 Teuchos::TimeMonitor normTimer( *timerNorm_ );
720 MVT::MvDot( X, *MX, diag );
722 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
724 if (SCT::magnitude((*B)(0,0)) > ZERO) {
726 MVT::MvScale( X, ONE/(*B)(0,0) );
729 MVT::MvScale( *MX, ONE/(*B)(0,0) );
736 dep_flg = blkOrtho( X, MX, C, Q );
742 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
745 MVT::Assign( *tmpX, X );
747 MVT::Assign( *tmpMX, *MX );
752 rank = findBasis( X, MX, B,
false );
757 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
760 MVT::Assign( *tmpX, X );
762 MVT::Assign( *tmpMX, *MX );
769 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
770 "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
780 template<
class ScalarType,
class MV,
class OP>
782 MV &X, Teuchos::RCP<MV> MX,
783 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
785 #ifdef BELOS_TEUCHOS_TIME_MONITOR 786 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
790 return findBasis(X, MX, B,
true);
796 template<
class ScalarType,
class MV,
class OP>
798 MV &X, Teuchos::RCP<MV> MX,
799 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
800 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
816 #ifdef BELOS_TEUCHOS_TIME_MONITOR 817 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
820 int xc = MVT::GetNumberVecs( X );
821 ptrdiff_t xr = MVT::GetGlobalLength( X );
823 std::vector<int> qcs(nq);
825 if (nq == 0 || xc == 0 || xr == 0) {
828 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
837 if (MX == Teuchos::null) {
839 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
840 OPT::Apply(*(this->_Op),X,*MX);
845 MX = Teuchos::rcp( &X,
false );
847 int mxc = MVT::GetNumberVecs( *MX );
848 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
851 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
852 "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
854 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
855 "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
859 for (
int i=0; i<nq; i++) {
860 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
861 "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
862 qcs[i] = MVT::GetNumberVecs( *Q[i] );
863 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
864 "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
868 if ( C[i] == Teuchos::null ) {
869 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
872 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
873 "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
878 blkOrtho( X, MX, C, Q );
885 template<
class ScalarType,
class MV,
class OP>
887 MV &X, Teuchos::RCP<MV> MX,
888 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
889 bool completeBasis,
int howMany )
const {
906 const ScalarType ONE = SCT::one();
909 int xc = MVT::GetNumberVecs( X );
910 ptrdiff_t xr = MVT::GetGlobalLength( X );
923 if (MX == Teuchos::null) {
925 MX = MVT::Clone(X,xc);
926 OPT::Apply(*(this->_Op),X,*MX);
933 if ( B == Teuchos::null ) {
934 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
937 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
938 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
941 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
942 "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
943 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
944 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
945 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
946 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
947 TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
948 "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
949 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
950 "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
955 int xstart = xc - howMany;
957 for (
int j = xstart; j < xc; j++) {
966 std::vector<int> index(1);
968 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
969 Teuchos::RCP<MV> MXj;
970 if ((this->_hasOp)) {
972 MXj = MVT::CloneViewNonConst( *MX, index );
979 Teuchos::RCP<MV> oldMXj;
981 oldMXj = MVT::CloneCopy( *MXj );
985 Teuchos::SerialDenseVector<int,ScalarType> product(numX);
986 Teuchos::SerialDenseVector<int,ScalarType> P2(1);
987 Teuchos::RCP<const MV> prevX, prevMX;
989 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
994 #ifdef BELOS_TEUCHOS_TIME_MONITOR 995 Teuchos::TimeMonitor normTimer( *timerNorm_ );
997 MVT::MvDot( *Xj, *MXj, oldDot );
1000 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO,
OrthoError,
1001 "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1004 for (
int ii=0; ii<numX; ii++) {
1007 prevX = MVT::CloneView( X, index );
1009 prevMX = MVT::CloneView( *MX, index );
1012 for (
int i=0; i<max_ortho_steps_; ++i) {
1016 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1017 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1025 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1026 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1028 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1036 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1037 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1039 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1044 product[ii] = P2[0];
1046 product[ii] += P2[0];
1054 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1055 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1057 MVT::MvDot( *Xj, *oldMXj, newDot );
1060 newDot[0] = oldDot[0];
1064 if (completeBasis) {
1068 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1073 std::cout <<
"Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1076 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1077 Teuchos::RCP<MV> tempMXj;
1078 MVT::MvRandom( *tempXj );
1080 tempMXj = MVT::Clone( X, 1 );
1081 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1087 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1088 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1090 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1094 for (
int ii=0; ii<numX; ii++) {
1097 prevX = MVT::CloneView( X, index );
1099 prevMX = MVT::CloneView( *MX, index );
1102 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1104 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1105 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1110 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1111 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1113 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
1116 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1117 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1119 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
1124 product[ii] = P2[0];
1126 product[ii] += P2[0];
1132 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1133 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1135 MVT::MvDot( *tempXj, *tempMXj, newDot );
1138 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1140 MVT::Assign( *tempXj, *Xj );
1142 MVT::Assign( *tempMXj, *MXj );
1154 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1163 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1164 if (SCT::magnitude(diag) > ZERO) {
1165 MVT::MvScale( *Xj, ONE/diag );
1168 MVT::MvScale( *MXj, ONE/diag );
1182 for (
int i=0; i<numX; i++) {
1183 (*B)(i,j) = product(i);
1194 template<
class ScalarType,
class MV,
class OP>
1197 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1198 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1201 int xc = MVT::GetNumberVecs( X );
1202 const ScalarType ONE = SCT::one();
1204 std::vector<int> qcs( nq );
1205 for (
int i=0; i<nq; i++) {
1206 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1210 std::vector<int> index(1);
1211 Teuchos::RCP<const MV> tempQ;
1213 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1215 for (
int i=0; i<nq; i++) {
1218 for (
int ii=0; ii<qcs[i]; ii++) {
1221 tempQ = MVT::CloneView( *Q[i], index );
1222 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1226 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1227 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1233 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1234 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1236 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1243 OPT::Apply( *(this->_Op), X, *MX);
1247 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1248 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1249 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1255 for (
int j = 1; j < max_ortho_steps_; ++j) {
1257 for (
int i=0; i<nq; i++) {
1259 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1262 for (
int ii=0; ii<qcs[i]; ii++) {
1265 tempQ = MVT::CloneView( *Q[i], index );
1266 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1267 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1271 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1272 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1278 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1279 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1281 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1290 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1292 else if (xc <= qcs[i]) {
1294 OPT::Apply( *(this->_Op), X, *MX);
1305 template<
class ScalarType,
class MV,
class OP>
1308 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1309 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1312 int xc = MVT::GetNumberVecs( X );
1313 bool dep_flg =
false;
1314 const ScalarType ONE = SCT::one();
1316 std::vector<int> qcs( nq );
1317 for (
int i=0; i<nq; i++) {
1318 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1324 std::vector<ScalarType> oldDot( xc );
1326 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1327 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1329 MVT::MvDot( X, *MX, oldDot );
1332 std::vector<int> index(1);
1333 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1334 Teuchos::RCP<const MV> tempQ;
1337 for (
int i=0; i<nq; i++) {
1340 for (
int ii=0; ii<qcs[i]; ii++) {
1343 tempQ = MVT::CloneView( *Q[i], index );
1344 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1348 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1349 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1355 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1356 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1358 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1365 OPT::Apply( *(this->_Op), X, *MX);
1369 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1370 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1371 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1377 for (
int j = 1; j < max_ortho_steps_; ++j) {
1379 for (
int i=0; i<nq; i++) {
1380 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1383 for (
int ii=0; ii<qcs[i]; ii++) {
1386 tempQ = MVT::CloneView( *Q[i], index );
1387 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1388 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1392 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1393 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1399 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1400 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1402 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1410 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1412 else if (xc <= qcs[i]) {
1414 OPT::Apply( *(this->_Op), X, *MX);
1421 std::vector<ScalarType> newDot(xc);
1423 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1424 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1426 MVT::MvDot( X, *MX, newDot );
1430 for (
int i=0; i<xc; i++){
1431 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1440 template<
class ScalarType,
class MV,
class OP>
1443 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1444 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1445 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1447 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1449 const ScalarType ONE = SCT::one();
1450 const ScalarType ZERO = SCT::zero();
1453 int xc = MVT::GetNumberVecs( X );
1454 std::vector<int> indX( 1 );
1455 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1457 std::vector<int> qcs( nq );
1458 for (
int i=0; i<nq; i++) {
1459 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1463 Teuchos::RCP<const MV> lastQ;
1464 Teuchos::RCP<MV> Xj, MXj;
1465 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1468 for (
int j=0; j<xc; j++) {
1470 bool dep_flg =
false;
1474 std::vector<int> index( j );
1475 for (
int ind=0; ind<j; ind++) {
1478 lastQ = MVT::CloneView( X, index );
1481 Q.push_back( lastQ );
1483 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1488 Xj = MVT::CloneViewNonConst( X, indX );
1490 MXj = MVT::CloneViewNonConst( *MX, indX );
1498 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1499 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1501 MVT::MvDot( *Xj, *MXj, oldDot );
1504 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1505 Teuchos::RCP<const MV> tempQ;
1508 for (
int i=0; i<Q.size(); i++) {
1511 for (
int ii=0; ii<qcs[i]; ii++) {
1514 tempQ = MVT::CloneView( *Q[i], indX );
1516 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1522 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1528 OPT::Apply( *(this->_Op), *Xj, *MXj);
1532 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1533 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1534 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1535 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1541 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1543 for (
int i=0; i<Q.size(); i++) {
1544 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1547 for (
int ii=0; ii<qcs[i]; ii++) {
1550 tempQ = MVT::CloneView( *Q[i], indX );
1552 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1556 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1560 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1567 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1569 else if (xc <= qcs[i]) {
1571 OPT::Apply( *(this->_Op), *Xj, *MXj);
1580 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1581 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1583 MVT::MvDot( *Xj, *MXj, newDot );
1587 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1593 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1595 MVT::MvScale( *Xj, ONE/diag );
1598 MVT::MvScale( *MXj, ONE/diag );
1606 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1607 Teuchos::RCP<MV> tempMXj;
1608 MVT::MvRandom( *tempXj );
1610 tempMXj = MVT::Clone( X, 1 );
1611 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1617 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1618 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1620 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1623 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1625 for (
int i=0; i<Q.size(); i++) {
1626 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1629 for (
int ii=0; ii<qcs[i]; ii++) {
1632 tempQ = MVT::CloneView( *Q[i], indX );
1633 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1637 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1645 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1647 else if (xc <= qcs[i]) {
1649 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1657 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1658 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1660 MVT::MvDot( *tempXj, *tempMXj, newDot );
1664 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1665 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1671 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1673 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1693 template<
class ScalarType,
class MV,
class OP>
1696 using Teuchos::ParameterList;
1697 using Teuchos::parameterList;
1700 RCP<ParameterList> params = parameterList (
"IMGS");
1705 "Maximum number of orthogonalization passes (includes the " 1706 "first). Default is 2, since \"twice is enough\" for Krylov " 1709 "Block reorthogonalization threshold.");
1711 "Singular block detection threshold.");
1716 template<
class ScalarType,
class MV,
class OP>
1719 using Teuchos::ParameterList;
1722 RCP<ParameterList> params = getIMGSDefaultParameters<ScalarType, MV, OP>();
1724 params->set (
"maxNumOrthogPasses",
1726 params->set (
"blkTol",
1728 params->set (
"singTol",
1736 #endif // BELOS_IMGS_ORTHOMANAGER_HPP void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
MagnitudeType blk_tol_
Block reorthogonalization tolerance.
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormError(const MV &X) const
This method computes the error in orthonormality of a multivector, measured as the Frobenius norm of ...
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.
Declaration of basic traits for the multivector type.
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...
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 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 ...
Class which defines basic traits for the operator type.
Traits class which defines basic operations on multivectors.
IMGSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_ortho_steps=max_ortho_steps_default_, const MagnitudeType blk_tol=blk_tol_default_, const MagnitudeType sing_tol=sing_tol_default_)
Constructor specifying re-orthogonalization tolerance.
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.
std::string label_
Label for timers.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::ScalarTraits< MagnitudeType > MGT
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.
Exception thrown to signal error in an orthogonalization manager method.
MagnitudeType sing_tol_
Singular block detection 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.
int max_ortho_steps_
Max number of (re)orthogonalization steps, including the first.
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.
MultiVecTraits< ScalarType, MV > MVT
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
Teuchos::RCP< Teuchos::ParameterList > getIMGSDefaultParameters()
"Default" parameters for robustness and accuracy.
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
static const int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
IMGSOrthoManager(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.
~IMGSOrthoManager()
Destructor.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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 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...
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
Belos header file which uses auto-configuration information to include necessary C++ headers...
Teuchos::RCP< Teuchos::ParameterList > getIMGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
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.
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.