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();
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
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();
631 B = rcp (
new serial_dense_matrix_type (xc, xc));
641 for (size_type k = 0; k < nq; ++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) {
669 MX = Teuchos::rcp( &X,
false );
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++) {
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;
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);
722 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
735 dep_flg = blkOrtho( X, MX, C, Q );
741 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
751 rank = findBasis( X, MX, B,
false );
756 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
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_);
822 std::vector<int> qcs(nq);
824 if (nq == 0 || xc == 0 || xr == 0) {
836 if (MX == Teuchos::null) {
844 MX = Teuchos::rcp( &X,
false );
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++) {
860 "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
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());
926 if (MX == Teuchos::null) {
936 if ( B == Teuchos::null ) {
937 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
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);
972 Teuchos::RCP<MV> MXj;
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 );
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++) {
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_ );
1030 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1031 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1038 product[ii] = P2[0];
1040 product[ii] += P2[0];
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;
1075 for (
int ii=0; ii<numX; ii++) {
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_ );
1097 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1098 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1105 product[ii] = P2[0];
1107 product[ii] += P2[0];
1114 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1130 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1139 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1158 for (
int i=0; i<numX; i++) {
1159 (*B)(i,j) = product(i);
1170 template<
class ScalarType,
class MV,
class OP>
1173 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1174 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1178 const ScalarType ONE = SCT::one();
1180 std::vector<int> qcs( nq );
1181 for (
int i=0; i<nq; 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++) {
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_ );
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++) {
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_ );
1268 else if (xc <= qcs[i]) {
1281 template<
class ScalarType,
class MV,
class OP>
1284 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1285 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 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++) {
1300 std::vector<ScalarType> oldDot( xc );
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++) {
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_ );
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++) {
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_ );
1383 else if (xc <= qcs[i]) {
1392 std::vector<ScalarType> newDot(xc);
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>
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();
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++) {
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++) {
1447 Q.push_back( lastQ );
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++) {
1477 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1495 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
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++) {
1513 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1521 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1530 else if (xc <= qcs[i]) {
1540 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1546 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1559 Teuchos::RCP<MV> tempXj =
MVT::Clone( X, 1 );
1560 Teuchos::RCP<MV> tempMXj;
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++) {
1581 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1595 else if (xc <= qcs[i]) {
1607 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1608 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1616 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1638 #endif // BELOS_IMGS_ORTHOMANAGER_HPP void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
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...
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
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.
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 ...
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Creates a new const MV that shares the selected contents of mv (shallow copy).
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
static void Assign(const MV &A, MV &mv)
mv := A
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Replace mv with .
Traits class which defines basic operations on multivectors.
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...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.
static void Apply(const OP &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Apply Op to x, putting the result into y.
Teuchos::RCP< const OP > _Op
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 ...
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
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.
Exception thrown to signal error in an orthogonalization manager method.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
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 ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
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 ...
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.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
Class which defines basic traits for the operator type.
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 void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
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 ...
~IMGSOrthoManager()
Destructor.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
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< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...