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_;
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;
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." );
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" );
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" );
866 if ( C[i] == Teuchos::null ) {
867 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
870 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
871 "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
876 blkOrtho( X, MX, C, Q );
885 MV &X, Teuchos::RCP<MV> MX,
886 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
887 bool completeBasis,
int howMany )
const {
904 const ScalarType ONE = SCT::one();
907 int xc = MVT::GetNumberVecs( X );
908 ptrdiff_t xr = MVT::GetGlobalLength( X );
921 if (MX == Teuchos::null) {
923 MX = MVT::Clone(X,xc);
924 OPT::Apply(*(this->_Op),X,*MX);
931 if ( B == Teuchos::null ) {
932 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
935 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
936 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
939 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
940 "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
941 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
942 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
943 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
944 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
945 TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
946 "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
947 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
948 "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
953 int xstart = xc - howMany;
955 for (
int j = xstart; j < xc; j++) {
964 std::vector<int> index(1);
966 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
967 Teuchos::RCP<MV> MXj;
968 if ((this->_hasOp)) {
970 MXj = MVT::CloneViewNonConst( *MX, index );
977 Teuchos::RCP<MV> oldMXj;
979 oldMXj = MVT::CloneCopy( *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 );
992#ifdef BELOS_TEUCHOS_TIME_MONITOR
993 Teuchos::TimeMonitor normTimer( *timerNorm_ );
995 MVT::MvDot( *Xj, *MXj, oldDot );
998 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO,
OrthoError,
999 "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1002 for (
int ii=0; ii<numX; ii++) {
1005 prevX = MVT::CloneView( X, index );
1007 prevMX = MVT::CloneView( *MX, index );
1010 for (
int i=0; i<max_ortho_steps_; ++i) {
1014#ifdef BELOS_TEUCHOS_TIME_MONITOR
1015 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1023#ifdef BELOS_TEUCHOS_TIME_MONITOR
1024 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1026 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1034#ifdef BELOS_TEUCHOS_TIME_MONITOR
1035 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1037 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1042 product[ii] = P2[0];
1044 product[ii] += P2[0];
1052#ifdef BELOS_TEUCHOS_TIME_MONITOR
1053 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1055 MVT::MvDot( *Xj, *oldMXj, newDot );
1058 newDot[0] = oldDot[0];
1062 if (completeBasis) {
1066 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1071 std::cout <<
"Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1074 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1075 Teuchos::RCP<MV> tempMXj;
1076 MVT::MvRandom( *tempXj );
1078 tempMXj = MVT::Clone( X, 1 );
1079 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1085#ifdef BELOS_TEUCHOS_TIME_MONITOR
1086 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1088 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1092 for (
int ii=0; ii<numX; ii++) {
1095 prevX = MVT::CloneView( X, index );
1097 prevMX = MVT::CloneView( *MX, index );
1100 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1102#ifdef BELOS_TEUCHOS_TIME_MONITOR
1103 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1108#ifdef BELOS_TEUCHOS_TIME_MONITOR
1109 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1111 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *tempXj );
1114#ifdef BELOS_TEUCHOS_TIME_MONITOR
1115 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1117 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *tempMXj );
1122 product[ii] = P2[0];
1124 product[ii] += P2[0];
1130#ifdef BELOS_TEUCHOS_TIME_MONITOR
1131 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1133 MVT::MvDot( *tempXj, *tempMXj, newDot );
1136 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1138 MVT::Assign( *tempXj, *Xj );
1140 MVT::Assign( *tempMXj, *MXj );
1152 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1161 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1162 if (SCT::magnitude(diag) > ZERO) {
1163 MVT::MvScale( *Xj, ONE/diag );
1166 MVT::MvScale( *MXj, ONE/diag );
1180 for (
int i=0; i<numX; i++) {
1181 (*B)(i,j) = product(i);
1195 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1196 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
1199 int xc = MVT::GetNumberVecs( X );
1200 const ScalarType ONE = SCT::one();
1202 std::vector<int> qcs( nq );
1203 for (
int i=0; i<nq; i++) {
1204 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1208 std::vector<int> index(1);
1209 Teuchos::RCP<const MV> tempQ;
1211 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1213 for (
int i=0; i<nq; i++) {
1216 for (
int ii=0; ii<qcs[i]; ii++) {
1219 tempQ = MVT::CloneView( *Q[i], index );
1220 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1224#ifdef BELOS_TEUCHOS_TIME_MONITOR
1225 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1231#ifdef BELOS_TEUCHOS_TIME_MONITOR
1232 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1234 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1241 OPT::Apply( *(this->_Op), X, *MX);
1245 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1246 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1247 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1253 for (
int j = 1; j < max_ortho_steps_; ++j) {
1255 for (
int i=0; i<nq; i++) {
1257 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1260 for (
int ii=0; ii<qcs[i]; ii++) {
1263 tempQ = MVT::CloneView( *Q[i], index );
1264 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1265 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1269#ifdef BELOS_TEUCHOS_TIME_MONITOR
1270 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1276#ifdef BELOS_TEUCHOS_TIME_MONITOR
1277 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1279 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1288 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1290 else if (xc <= qcs[i]) {
1292 OPT::Apply( *(this->_Op), X, *MX);
1306 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1307 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
1310 int xc = MVT::GetNumberVecs( X );
1311 bool dep_flg =
false;
1312 const ScalarType ONE = SCT::one();
1314 std::vector<int> qcs( nq );
1315 for (
int i=0; i<nq; i++) {
1316 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1322 std::vector<ScalarType> oldDot( xc );
1324#ifdef BELOS_TEUCHOS_TIME_MONITOR
1325 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1327 MVT::MvDot( X, *MX, oldDot );
1330 std::vector<int> index(1);
1331 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1332 Teuchos::RCP<const MV> tempQ;
1335 for (
int i=0; i<nq; i++) {
1338 for (
int ii=0; ii<qcs[i]; ii++) {
1341 tempQ = MVT::CloneView( *Q[i], index );
1342 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1346#ifdef BELOS_TEUCHOS_TIME_MONITOR
1347 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1353#ifdef BELOS_TEUCHOS_TIME_MONITOR
1354 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1356 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1363 OPT::Apply( *(this->_Op), X, *MX);
1367 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1368 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1369 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1375 for (
int j = 1; j < max_ortho_steps_; ++j) {
1377 for (
int i=0; i<nq; i++) {
1378 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1381 for (
int ii=0; ii<qcs[i]; ii++) {
1384 tempQ = MVT::CloneView( *Q[i], index );
1385 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1386 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1390#ifdef BELOS_TEUCHOS_TIME_MONITOR
1391 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1397#ifdef BELOS_TEUCHOS_TIME_MONITOR
1398 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1400 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1408 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1410 else if (xc <= qcs[i]) {
1412 OPT::Apply( *(this->_Op), X, *MX);
1419 std::vector<ScalarType> newDot(xc);
1421#ifdef BELOS_TEUCHOS_TIME_MONITOR
1422 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1424 MVT::MvDot( X, *MX, newDot );
1428 for (
int i=0; i<xc; i++){
1429 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1441 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1442 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1443 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const
1445 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1447 const ScalarType ONE = SCT::one();
1448 const ScalarType ZERO = SCT::zero();
1451 int xc = MVT::GetNumberVecs( X );
1452 std::vector<int> indX( 1 );
1453 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1455 std::vector<int> qcs( nq );
1456 for (
int i=0; i<nq; i++) {
1457 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1461 Teuchos::RCP<const MV> lastQ;
1462 Teuchos::RCP<MV> Xj, MXj;
1463 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1466 for (
int j=0; j<xc; j++) {
1468 bool dep_flg =
false;
1472 std::vector<int> index( j );
1473 for (
int ind=0; ind<j; ind++) {
1476 lastQ = MVT::CloneView( X, index );
1479 Q.push_back( lastQ );
1481 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1486 Xj = MVT::CloneViewNonConst( X, indX );
1488 MXj = MVT::CloneViewNonConst( *MX, indX );
1496#ifdef BELOS_TEUCHOS_TIME_MONITOR
1497 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1499 MVT::MvDot( *Xj, *MXj, oldDot );
1502 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1503 Teuchos::RCP<const MV> tempQ;
1506 for (
int i=0; i<Q.size(); i++) {
1509 for (
int ii=0; ii<qcs[i]; ii++) {
1512 tempQ = MVT::CloneView( *Q[i], indX );
1514 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, j );
1520 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1526 OPT::Apply( *(this->_Op), *Xj, *MXj);
1530 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1531 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1532 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1533 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1539 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1541 for (
int i=0; i<Q.size(); i++) {
1542 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1545 for (
int ii=0; ii<qcs[i]; ii++) {
1548 tempQ = MVT::CloneView( *Q[i], indX );
1550 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1554 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1558 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1565 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1567 else if (xc <= qcs[i]) {
1569 OPT::Apply( *(this->_Op), *Xj, *MXj);
1578#ifdef BELOS_TEUCHOS_TIME_MONITOR
1579 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1581 MVT::MvDot( *Xj, *MXj, newDot );
1585 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1591 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1593 MVT::MvScale( *Xj, ONE/diag );
1596 MVT::MvScale( *MXj, ONE/diag );
1604 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1605 Teuchos::RCP<MV> tempMXj;
1606 MVT::MvRandom( *tempXj );
1608 tempMXj = MVT::Clone( X, 1 );
1609 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1615#ifdef BELOS_TEUCHOS_TIME_MONITOR
1616 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1618 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1621 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1623 for (
int i=0; i<Q.size(); i++) {
1624 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1627 for (
int ii=0; ii<qcs[i]; ii++) {
1630 tempQ = MVT::CloneView( *Q[i], indX );
1631 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1635 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1643 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1645 else if (xc <= qcs[i]) {
1647 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1655#ifdef BELOS_TEUCHOS_TIME_MONITOR
1656 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1658 MVT::MvDot( *tempXj, *tempMXj, newDot );
1662 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1663 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1669 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1671 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );