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 #include "Teuchos_ParameterListAcceptorDefaultBase.hpp" 68 #ifdef BELOS_TEUCHOS_TIME_MONITOR 69 #include "Teuchos_TimeMonitor.hpp" 70 #endif // BELOS_TEUCHOS_TIME_MONITOR 74 template<
class ScalarType,
class MV,
class OP>
77 public Teuchos::ParameterListAcceptorDefaultBase
80 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
81 typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
82 typedef Teuchos::ScalarTraits<ScalarType> SCT;
92 Teuchos::RCP<const OP> Op = Teuchos::null,
93 const int max_ortho_steps = 1,
94 const MagnitudeType blk_tol = 10*MGT::squareroot( MGT::eps() ),
95 const MagnitudeType sing_tol = 10*MGT::eps() )
97 max_ortho_steps_( max_ortho_steps ),
99 sing_tol_( sing_tol ),
102 std::string orthoLabel = label_ +
": Orthogonalization";
103 #ifdef BELOS_TEUCHOS_TIME_MONITOR 104 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
107 std::string updateLabel = label_ +
": Ortho (Update)";
108 #ifdef BELOS_TEUCHOS_TIME_MONITOR 109 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
112 std::string normLabel = label_ +
": Ortho (Norm)";
113 #ifdef BELOS_TEUCHOS_TIME_MONITOR 114 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
117 std::string ipLabel = label_ +
": Ortho (Inner Product)";
118 #ifdef BELOS_TEUCHOS_TIME_MONITOR 119 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
125 const std::string& label =
"Belos",
126 Teuchos::RCP<const OP> Op = Teuchos::null) :
128 max_ortho_steps_ (2),
129 blk_tol_ (10 * MGT::squareroot (MGT::eps())),
130 sing_tol_ (10 * MGT::eps()),
135 std::string orthoLabel = label_ +
": Orthogonalization";
136 #ifdef BELOS_TEUCHOS_TIME_MONITOR 137 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
139 std::string updateLabel = label_ +
": Ortho (Update)";
140 #ifdef BELOS_TEUCHOS_TIME_MONITOR 141 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
143 std::string normLabel = label_ +
": Ortho (Norm)";
144 #ifdef BELOS_TEUCHOS_TIME_MONITOR 145 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
147 std::string ipLabel = label_ +
": Ortho (Inner Product)";
148 #ifdef BELOS_TEUCHOS_TIME_MONITOR 149 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
162 using Teuchos::Exceptions::InvalidParameterName;
163 using Teuchos::ParameterList;
164 using Teuchos::parameterList;
168 RCP<ParameterList> params;
169 if (plist.is_null()) {
170 params = parameterList (*defaultParams);
185 int maxNumOrthogPasses;
186 MagnitudeType blkTol;
187 MagnitudeType singTol;
190 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
191 }
catch (InvalidParameterName&) {
192 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
193 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
204 blkTol = params->get<MagnitudeType> (
"blkTol");
205 }
catch (InvalidParameterName&) {
207 blkTol = params->get<MagnitudeType> (
"depTol");
210 params->remove (
"depTol");
211 }
catch (InvalidParameterName&) {
212 blkTol = defaultParams->get<MagnitudeType> (
"blkTol");
214 params->set (
"blkTol", blkTol);
218 singTol = params->get<MagnitudeType> (
"singTol");
219 }
catch (InvalidParameterName&) {
220 singTol = defaultParams->get<MagnitudeType> (
"singTol");
221 params->set (
"singTol", singTol);
224 max_ortho_steps_ = maxNumOrthogPasses;
228 setMyParamList (params);
231 Teuchos::RCP<const Teuchos::ParameterList>
235 using Teuchos::ParameterList;
236 using Teuchos::parameterList;
239 if (defaultParams_.is_null()) {
240 RCP<ParameterList> params = parameterList (
"IMGS");
244 const int defaultMaxNumOrthogPasses = 2;
245 const MagnitudeType eps = MGT::eps();
246 const MagnitudeType defaultBlkTol =
247 as<MagnitudeType> (10) * MGT::squareroot (eps);
248 const MagnitudeType defaultSingTol = as<MagnitudeType> (10) * eps;
250 params->set (
"maxNumOrthogPasses", defaultMaxNumOrthogPasses,
251 "Maximum number of orthogonalization passes (includes the " 252 "first). Default is 2, since \"twice is enough\" for Krylov " 254 params->set (
"blkTol", defaultBlkTol,
"Block reorthogonalization " 256 params->set (
"singTol", defaultSingTol,
"Singular block detection " 258 defaultParams_ = params;
260 return defaultParams_;
268 Teuchos::RCP<const Teuchos::ParameterList>
272 using Teuchos::ParameterList;
273 using Teuchos::parameterList;
278 RCP<ParameterList> params = parameterList (*defaultParams);
280 const int maxBlkOrtho = 1;
281 const MagnitudeType blkTol = MGT::zero();
282 const MagnitudeType singTol = MGT::zero();
284 params->set (
"maxNumOrthogPasses", maxBlkOrtho);
285 params->set (
"blkTol", blkTol);
286 params->set (
"singTol", singTol);
295 void setBlkTol(
const MagnitudeType blk_tol ) { blk_tol_ = blk_tol; }
298 void setSingTol(
const MagnitudeType sing_tol ) { sing_tol_ = sing_tol; }
339 void project ( MV &X, Teuchos::RCP<MV> MX,
340 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
341 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
347 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
348 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
378 int normalize ( MV &X, Teuchos::RCP<MV> MX,
379 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
384 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
433 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
434 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
435 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
445 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
454 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
460 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
469 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
470 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
479 void setLabel(
const std::string& label);
483 const std::string&
getLabel()
const {
return label_; }
489 int max_ortho_steps_;
491 MagnitudeType blk_tol_;
493 MagnitudeType sing_tol_;
497 #ifdef BELOS_TEUCHOS_TIME_MONITOR 498 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_,
499 timerNorm_, timerScale_, timerInnerProd_;
500 #endif // BELOS_TEUCHOS_TIME_MONITOR 503 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
506 int findBasis(MV &X, Teuchos::RCP<MV> MX,
507 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
508 bool completeBasis,
int howMany = -1 )
const;
511 bool blkOrtho1 ( 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;
516 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
517 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
518 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
533 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
534 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
535 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
536 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
541 template<
class ScalarType,
class MV,
class OP>
544 if (label != label_) {
546 std::string orthoLabel = label_ +
": Orthogonalization";
547 #ifdef BELOS_TEUCHOS_TIME_MONITOR 548 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
551 std::string updateLabel = label_ +
": Ortho (Update)";
552 #ifdef BELOS_TEUCHOS_TIME_MONITOR 553 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
556 std::string normLabel = label_ +
": Ortho (Norm)";
557 #ifdef BELOS_TEUCHOS_TIME_MONITOR 558 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
561 std::string ipLabel = label_ +
": Ortho (Inner Product)";
562 #ifdef BELOS_TEUCHOS_TIME_MONITOR 563 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
570 template<
class ScalarType,
class MV,
class OP>
571 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
573 const ScalarType ONE = SCT::one();
574 int rank = MVT::GetNumberVecs(X);
575 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
577 for (
int i=0; i<rank; i++) {
580 return xTx.normFrobenius();
585 template<
class ScalarType,
class MV,
class OP>
586 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
588 int r1 = MVT::GetNumberVecs(X1);
589 int r2 = MVT::GetNumberVecs(X2);
590 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
592 return xTx.normFrobenius();
597 template<
class ScalarType,
class MV,
class OP>
602 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
603 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
604 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 606 using Teuchos::Array;
608 using Teuchos::is_null;
611 using Teuchos::SerialDenseMatrix;
612 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
613 typedef typename Array< RCP< const MV > >::size_type size_type;
615 #ifdef BELOS_TEUCHOS_TIME_MONITOR 616 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
619 ScalarType ONE = SCT::one();
623 int xc = MVT::GetNumberVecs( X );
624 ptrdiff_t xr = MVT::GetGlobalLength( X );
631 B = rcp (
new serial_dense_matrix_type (xc, xc));
641 for (size_type k = 0; k < nq; ++k)
643 const int numRows = MVT::GetNumberVecs (*Q[k]);
644 const int numCols = xc;
647 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
648 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
650 int err = C[k]->reshape (numRows, numCols);
651 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
652 "IMGS orthogonalization: failed to reshape " 653 "C[" << k <<
"] (the array of block " 654 "coefficients resulting from projecting X " 655 "against Q[1:" << nq <<
"]).");
661 if (MX == Teuchos::null) {
663 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
664 OPT::Apply(*(this->_Op),X,*MX);
669 MX = Teuchos::rcp( &X,
false );
672 int mxc = MVT::GetNumberVecs( *MX );
673 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
676 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
679 for (
int i=0; i<nq; i++) {
680 numbas += MVT::GetNumberVecs( *Q[i] );
684 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
685 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
687 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
688 "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
690 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
691 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
697 bool dep_flg =
false;
700 Teuchos::RCP<MV> tmpX, tmpMX;
701 tmpX = MVT::CloneCopy(X);
703 tmpMX = MVT::CloneCopy(*MX);
710 dep_flg = blkOrtho1( X, MX, C, Q );
714 #ifdef BELOS_TEUCHOS_TIME_MONITOR 715 Teuchos::TimeMonitor normTimer( *timerNorm_ );
717 if ( B == Teuchos::null ) {
718 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
720 std::vector<ScalarType> diag(xc);
721 MVT::MvDot( X, *MX, diag );
722 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
724 MVT::MvScale( X, ONE/(*B)(0,0) );
727 MVT::MvScale( *MX, ONE/(*B)(0,0) );
735 dep_flg = blkOrtho( X, MX, C, Q );
741 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
744 MVT::Assign( *tmpX, X );
746 MVT::Assign( *tmpMX, *MX );
751 rank = findBasis( X, MX, B,
false );
756 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
759 MVT::Assign( *tmpX, X );
761 MVT::Assign( *tmpMX, *MX );
768 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
769 "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
779 template<
class ScalarType,
class MV,
class OP>
781 MV &X, Teuchos::RCP<MV> MX,
782 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
784 #ifdef BELOS_TEUCHOS_TIME_MONITOR 785 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
789 return findBasis(X, MX, B,
true);
795 template<
class ScalarType,
class MV,
class OP>
797 MV &X, Teuchos::RCP<MV> MX,
798 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
799 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
815 #ifdef BELOS_TEUCHOS_TIME_MONITOR 816 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
819 int xc = MVT::GetNumberVecs( X );
820 ptrdiff_t xr = MVT::GetGlobalLength( X );
822 std::vector<int> qcs(nq);
824 if (nq == 0 || xc == 0 || xr == 0) {
827 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
836 if (MX == Teuchos::null) {
838 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
839 OPT::Apply(*(this->_Op),X,*MX);
844 MX = Teuchos::rcp( &X,
false );
846 int mxc = MVT::GetNumberVecs( *MX );
847 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
850 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
851 "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
853 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
854 "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
858 for (
int i=0; i<nq; i++) {
859 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
860 "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
861 qcs[i] = MVT::GetNumberVecs( *Q[i] );
862 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
863 "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
867 if ( C[i] == Teuchos::null ) {
868 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
871 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
872 "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
877 blkOrtho( X, MX, C, Q );
884 template<
class ScalarType,
class MV,
class OP>
886 MV &X, Teuchos::RCP<MV> MX,
887 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
888 bool completeBasis,
int howMany )
const {
905 #ifdef BELOS_TEUCHOS_TIME_MONITOR 906 Teuchos::TimeMonitor normTimer( *timerNorm_ );
909 const ScalarType ONE = SCT::one();
910 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
912 int xc = MVT::GetNumberVecs( X );
913 ptrdiff_t xr = MVT::GetGlobalLength( X );
926 if (MX == Teuchos::null) {
928 MX = MVT::Clone(X,xc);
929 OPT::Apply(*(this->_Op),X,*MX);
936 if ( B == Teuchos::null ) {
937 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
940 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
941 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
944 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
945 "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
946 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
947 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
948 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
949 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
950 TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
951 "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
952 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
953 "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
958 int xstart = xc - howMany;
960 for (
int j = xstart; j < xc; j++) {
969 std::vector<int> index(1);
971 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
972 Teuchos::RCP<MV> MXj;
973 if ((this->_hasOp)) {
975 MXj = MVT::CloneViewNonConst( *MX, index );
983 Teuchos::SerialDenseVector<int,ScalarType> product(numX);
984 Teuchos::SerialDenseVector<int,ScalarType> P2(1);
985 Teuchos::RCP<const MV> prevX, prevMX;
987 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
991 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
992 MVT::MvDot( *Xj, *MXj, oldDot );
994 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
995 "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
998 for (
int ii=0; ii<numX; ii++) {
1001 prevX = MVT::CloneView( X, index );
1003 prevMX = MVT::CloneView( *MX, index );
1006 for (
int i=0; i<max_ortho_steps_; ++i) {
1010 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1011 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1019 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1020 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1022 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1031 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1033 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1038 product[ii] = P2[0];
1040 product[ii] += P2[0];
1047 MVT::MvDot( *Xj, *oldMXj, newDot );
1050 if (completeBasis) {
1054 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1059 std::cout <<
"Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1062 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1063 Teuchos::RCP<MV> tempMXj;
1064 MVT::MvRandom( *tempXj );
1066 tempMXj = MVT::Clone( X, 1 );
1067 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1072 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1075 for (
int ii=0; ii<numX; ii++) {
1078 prevX = MVT::CloneView( X, index );
1080 prevMX = MVT::CloneView( *MX, index );
1083 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1085 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1086 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1091 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1092 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1094 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
1097 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1098 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1100 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
1105 product[ii] = P2[0];
1107 product[ii] += P2[0];
1112 MVT::MvDot( *tempXj, *tempMXj, newDot );
1114 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1116 MVT::Assign( *tempXj, *Xj );
1118 MVT::Assign( *tempMXj, *MXj );
1130 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1139 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1141 MVT::MvScale( *Xj, ONE/diag );
1144 MVT::MvScale( *MXj, ONE/diag );
1158 for (
int i=0; i<numX; i++) {
1159 (*B)(i,j) = product(i);
1170 template<
class ScalarType,
class MV,
class OP>
1172 IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1173 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1174 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1177 int xc = MVT::GetNumberVecs( X );
1178 const ScalarType ONE = SCT::one();
1180 std::vector<int> qcs( nq );
1181 for (
int i=0; i<nq; i++) {
1182 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1186 std::vector<int> index(1);
1187 Teuchos::RCP<const MV> tempQ;
1189 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1191 for (
int i=0; i<nq; i++) {
1194 for (
int ii=0; ii<qcs[i]; ii++) {
1197 tempQ = MVT::CloneView( *Q[i], index );
1198 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1202 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1203 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1209 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1210 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1212 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1219 OPT::Apply( *(this->_Op), X, *MX);
1223 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1224 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1225 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1231 for (
int j = 1; j < max_ortho_steps_; ++j) {
1233 for (
int i=0; i<nq; i++) {
1235 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1238 for (
int ii=0; ii<qcs[i]; ii++) {
1241 tempQ = MVT::CloneView( *Q[i], index );
1242 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1243 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1247 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1248 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1254 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1255 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1257 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1266 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1268 else if (xc <= qcs[i]) {
1270 OPT::Apply( *(this->_Op), X, *MX);
1281 template<
class ScalarType,
class MV,
class OP>
1283 IMGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1284 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1285 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1288 int xc = MVT::GetNumberVecs( X );
1289 bool dep_flg =
false;
1290 const ScalarType ONE = SCT::one();
1292 std::vector<int> qcs( nq );
1293 for (
int i=0; i<nq; i++) {
1294 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1300 std::vector<ScalarType> oldDot( xc );
1301 MVT::MvDot( X, *MX, oldDot );
1303 std::vector<int> index(1);
1304 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1305 Teuchos::RCP<const MV> tempQ;
1308 for (
int i=0; i<nq; i++) {
1311 for (
int ii=0; ii<qcs[i]; ii++) {
1314 tempQ = MVT::CloneView( *Q[i], index );
1315 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1319 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1320 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1326 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1327 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1329 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1336 OPT::Apply( *(this->_Op), X, *MX);
1340 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1341 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1342 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1348 for (
int j = 1; j < max_ortho_steps_; ++j) {
1350 for (
int i=0; i<nq; i++) {
1351 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1354 for (
int ii=0; ii<qcs[i]; ii++) {
1357 tempQ = MVT::CloneView( *Q[i], index );
1358 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1359 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1363 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1364 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1370 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1371 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1373 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1381 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1383 else if (xc <= qcs[i]) {
1385 OPT::Apply( *(this->_Op), X, *MX);
1392 std::vector<ScalarType> newDot(xc);
1393 MVT::MvDot( X, *MX, newDot );
1396 for (
int i=0; i<xc; i++){
1397 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1406 template<
class ScalarType,
class MV,
class OP>
1408 IMGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1409 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1410 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1411 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1413 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1415 const ScalarType ONE = SCT::one();
1416 const ScalarType ZERO = SCT::zero();
1419 int xc = MVT::GetNumberVecs( X );
1420 std::vector<int> indX( 1 );
1421 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1423 std::vector<int> qcs( nq );
1424 for (
int i=0; i<nq; i++) {
1425 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1429 Teuchos::RCP<const MV> lastQ;
1430 Teuchos::RCP<MV> Xj, MXj;
1431 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1434 for (
int j=0; j<xc; j++) {
1436 bool dep_flg =
false;
1440 std::vector<int> index( j );
1441 for (
int ind=0; ind<j; ind++) {
1444 lastQ = MVT::CloneView( X, index );
1447 Q.push_back( lastQ );
1449 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1454 Xj = MVT::CloneViewNonConst( X, indX );
1456 MXj = MVT::CloneViewNonConst( *MX, indX );
1463 MVT::MvDot( *Xj, *MXj, oldDot );
1465 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1466 Teuchos::RCP<const MV> tempQ;
1469 for (
int i=0; i<Q.size(); i++) {
1472 for (
int ii=0; ii<qcs[i]; ii++) {
1475 tempQ = MVT::CloneView( *Q[i], indX );
1477 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1483 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1489 OPT::Apply( *(this->_Op), *Xj, *MXj);
1493 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1494 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1495 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1496 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1502 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1504 for (
int i=0; i<Q.size(); i++) {
1505 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1508 for (
int ii=0; ii<qcs[i]; ii++) {
1511 tempQ = MVT::CloneView( *Q[i], indX );
1513 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1517 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1521 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1528 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1530 else if (xc <= qcs[i]) {
1532 OPT::Apply( *(this->_Op), *Xj, *MXj);
1540 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1546 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1548 MVT::MvScale( *Xj, ONE/diag );
1551 MVT::MvScale( *MXj, ONE/diag );
1559 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1560 Teuchos::RCP<MV> tempMXj;
1561 MVT::MvRandom( *tempXj );
1563 tempMXj = MVT::Clone( X, 1 );
1564 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1569 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1571 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1573 for (
int i=0; i<Q.size(); i++) {
1574 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1577 for (
int ii=0; ii<qcs[i]; ii++) {
1580 tempQ = MVT::CloneView( *Q[i], indX );
1581 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1585 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1593 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1595 else if (xc <= qcs[i]) {
1597 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1604 MVT::MvDot( *tempXj, *tempMXj, newDot );
1607 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1608 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1614 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1616 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1638 #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.
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 ...
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.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
IMGSOrthoManager(const std::string &label="Belos", Teuchos::RCP< const OP > Op=Teuchos::null, const int max_ortho_steps=1, const MagnitudeType blk_tol=10 *MGT::squareroot(MGT::eps()), const MagnitudeType sing_tol=10 *MGT::eps())
Constructor specifying re-orthogonalization tolerance.
void project(MV &X, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
This method calls project(X,Teuchos::null,C,Q); see documentation for that function.
void 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 ...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
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...
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.
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.