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"
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_;
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" );
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 );
883 template<
class ScalarType,
class MV,
class OP>
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);
1192 template<
class ScalarType,
class MV,
class OP>
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);
1303 template<
class ScalarType,
class MV,
class OP>
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_)) {
1438 template<
class ScalarType,
class MV,
class OP>
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 );
1691 template<
class ScalarType,
class MV,
class OP>
1694 using Teuchos::ParameterList;
1695 using Teuchos::parameterList;
1698 RCP<ParameterList> params = parameterList (
"IMGS");
1703 "Maximum number of orthogonalization passes (includes the "
1704 "first). Default is 2, since \"twice is enough\" for Krylov "
1707 "Block reorthogonalization threshold.");
1709 "Singular block detection threshold.");
1714 template<
class ScalarType,
class MV,
class OP>
1717 using Teuchos::ParameterList;
1720 RCP<ParameterList> params = getIMGSDefaultParameters<ScalarType, MV, OP>();
1722 params->set (
"maxNumOrthogPasses",
1724 params->set (
"blkTol",
1726 params->set (
"singTol",
Belos header file which uses auto-configuration information to include necessary C++ headers.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
An implementation of the Belos::MatOrthoManager that performs orthogonalization using multiple steps ...
MagnitudeType sing_tol_
Singular block detection threshold.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
int max_ortho_steps_
Max number of (re)orthogonalization steps, including the first.
MultiVecTraits< ScalarType, MV > MVT
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
Teuchos::ScalarTraits< ScalarType > SCT
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
int normalize(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B) const
This method calls normalize(X,Teuchos::null,B); see documentation for that function.
bool blkOrtho(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Routine to compute the block orthogonalization.
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.
~IMGSOrthoManager()
Destructor.
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 ...
MagnitudeType blk_tol_
Block reorthogonalization tolerance.
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 setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
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.
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::RCP< const Teuchos::ParameterList > getValidParameters() const
bool blkOrtho1(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Routine to compute the block orthogonalization.
void project(MV &X, Teuchos::RCP< MV > MX, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
Given a list of (mutually and internally) orthonormal bases Q, this method takes a multivector X and ...
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
static const int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
std::string label_
Label for timers.
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.
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
OperatorTraits< ScalarType, MV, OP > OPT
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 .
Teuchos::ScalarTraits< MagnitudeType > MGT
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.
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...
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
void innerProd(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z) const
Provides the inner product defining the orthogonality concepts, using the provided operator.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Exception thrown to signal error in an orthogonalization manager method.
Teuchos::RCP< Teuchos::ParameterList > getIMGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.
Teuchos::RCP< Teuchos::ParameterList > getIMGSDefaultParameters()
"Default" parameters for robustness and accuracy.