51 #ifndef BELOS_IMGS_ORTHOMANAGER_MP_VECTOR_HPP 52 #define BELOS_IMGS_ORTHOMANAGER_MP_VECTOR_HPP 55 #include "BelosIMGSOrthoManager.hpp" 59 template<
class Storage,
class MV,
class OP>
60 class IMGSOrthoManager<
Sacado::MP::Vector<Storage>,MV,OP> :
61 public MatOrthoManager<Sacado::MP::Vector<Storage>,MV,OP>
65 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
MagnitudeType;
66 typedef typename Teuchos::ScalarTraits<MagnitudeType>
MGT;
67 typedef Teuchos::ScalarTraits<ScalarType>
SCT;
68 typedef MultiVecTraits<ScalarType,MV>
MVT;
69 typedef OperatorTraits<ScalarType,MV,OP>
OPT;
77 Teuchos::RCP<const OP> Op = Teuchos::null,
78 const int max_ortho_steps = max_ortho_steps_default_,
82 max_ortho_steps_( max_ortho_steps ),
84 sing_tol_( sing_tol ),
87 #ifdef BELOS_TEUCHOS_TIME_MONITOR 89 ss << label_ +
": IMGS[" << max_ortho_steps_ <<
"]";
91 std::string orthoLabel = ss.str() +
": Orthogonalization";
92 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
94 std::string updateLabel = ss.str() +
": Ortho (Update)";
95 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
97 std::string normLabel = ss.str() +
": Ortho (Norm)";
98 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
100 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
101 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
107 const std::string& label =
"Belos",
108 Teuchos::RCP<const OP> Op = Teuchos::null) :
110 max_ortho_steps_ (max_ortho_steps_default_),
111 blk_tol_ (blk_tol_default_),
112 sing_tol_ (sing_tol_default_),
115 setParameterList (plist);
117 #ifdef BELOS_TEUCHOS_TIME_MONITOR 118 std::stringstream ss;
119 ss << label_ +
": IMGS[" << max_ortho_steps_ <<
"]";
121 std::string orthoLabel = ss.str() +
": Orthogonalization";
122 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
124 std::string updateLabel = ss.str() +
": Ortho (Update)";
125 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
127 std::string normLabel = ss.str() +
": Ortho (Norm)";
128 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
130 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
131 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
144 using Teuchos::Exceptions::InvalidParameterName;
145 using Teuchos::ParameterList;
146 using Teuchos::parameterList;
149 RCP<const ParameterList> defaultParams = getValidParameters();
150 RCP<ParameterList> params;
151 if (plist.is_null()) {
152 params = parameterList (*defaultParams);
167 int maxNumOrthogPasses;
172 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
173 }
catch (InvalidParameterName&) {
174 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
175 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
187 }
catch (InvalidParameterName&) {
192 params->remove (
"depTol");
193 }
catch (InvalidParameterName&) {
196 params->set (
"blkTol", blkTol);
201 }
catch (InvalidParameterName&) {
203 params->set (
"singTol", singTol);
206 max_ortho_steps_ = maxNumOrthogPasses;
210 this->setMyParamList (params);
213 Teuchos::RCP<const Teuchos::ParameterList>
216 if (defaultParams_.is_null()) {
217 defaultParams_ = Belos::getIMGSDefaultParameters<ScalarType, MV, OP>();
220 return defaultParams_;
273 void project ( MV &X, Teuchos::RCP<MV> MX,
274 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
275 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
281 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
282 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
283 project(X,Teuchos::null,C,Q);
312 int normalize ( MV &X, Teuchos::RCP<MV> MX,
313 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
318 int normalize ( MV &X, Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
319 return normalize(X,Teuchos::null,B);
365 projectAndNormalizeWithMxImpl (MV &X,
367 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
368 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
369 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
379 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
381 return orthonormError(X,Teuchos::null);
388 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
389 orthonormError(
const MV &X, Teuchos::RCP<const MV> MX)
const;
394 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
396 return orthogError(X1,Teuchos::null,X2);
403 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
404 orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const;
413 void setLabel(
const std::string& label);
417 const std::string&
getLabel()
const {
return label_; }
451 #ifdef BELOS_TEUCHOS_TIME_MONITOR 452 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
453 #endif // BELOS_TEUCHOS_TIME_MONITOR 459 int findBasis(MV &X, Teuchos::RCP<MV> MX,
460 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
461 bool completeBasis,
int howMany = -1 )
const;
464 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
465 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
466 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
469 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
470 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
471 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
486 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
487 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
488 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
489 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
493 template<
class StorageType,
class MV,
class OP>
494 const int IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::max_ortho_steps_default_ = 1;
496 template<
class StorageType,
class MV,
class OP>
497 const typename IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::MagnitudeType
498 IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::blk_tol_default_
499 = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::MagnitudeType>::squareroot(
502 template<
class StorageType,
class MV,
class OP>
503 const typename IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::MagnitudeType
504 IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::sing_tol_default_
505 = 10*Teuchos::ScalarTraits<typename IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::MagnitudeType>::eps();
507 template<
class StorageType,
class MV,
class OP>
508 const int IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::max_ortho_steps_fast_ = 1;
510 template<
class StorageType,
class MV,
class OP>
511 const typename IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::MagnitudeType
512 IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::blk_tol_fast_
513 = Teuchos::ScalarTraits<typename IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::MagnitudeType>::zero();
515 template<
class StorageType,
class MV,
class OP>
516 const typename IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::MagnitudeType
517 IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::sing_tol_fast_
518 = Teuchos::ScalarTraits<typename IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::MagnitudeType>::zero();
522 template<
class StorageType,
class MV,
class OP>
523 void IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::setLabel(
const std::string& label)
525 if (label != label_) {
527 #ifdef BELOS_TEUCHOS_TIME_MONITOR 528 std::stringstream ss;
529 ss << label_ +
": IMGS[" << max_ortho_steps_ <<
"]";
531 std::string orthoLabel = ss.str() +
": Orthogonalization";
532 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
534 std::string updateLabel = ss.str() +
": Ortho (Update)";
535 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
537 std::string normLabel = ss.str() +
": Ortho (Norm)";
538 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
540 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
541 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
548 template<
class StorageType,
class MV,
class OP>
549 typename Teuchos::ScalarTraits<Sacado::MP::Vector<StorageType> >::magnitudeType
550 IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::orthonormError(
const MV &X, Teuchos::RCP<const MV> MX)
const {
551 const ScalarType ONE = SCT::one();
552 int rank = MVT::GetNumberVecs(X);
553 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
554 MatOrthoManager<ScalarType,MV,OP>::innerProd(X,X,MX,xTx);
555 for (
int i=0; i<rank; i++) {
558 return xTx.normFrobenius();
563 template<
class StorageType,
class MV,
class OP>
564 typename Teuchos::ScalarTraits<Sacado::MP::Vector<StorageType> >::magnitudeType
565 IMGSOrthoManager<Sacado::MP::Vector<StorageType>,MV,OP>::orthogError(
const MV &X1, Teuchos::RCP<const MV> MX1,
const MV &X2)
const {
566 int r1 = MVT::GetNumberVecs(X1);
567 int r2 = MVT::GetNumberVecs(X2);
568 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
569 MatOrthoManager<ScalarType,MV,OP>::innerProd(X2,X1,MX1,xTx);
570 return xTx.normFrobenius();
575 template<
class StorageType,
class MV,
class OP>
577 IMGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::
578 projectAndNormalizeWithMxImpl(MV &X,
580 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
581 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
582 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 584 using Teuchos::Array;
586 using Teuchos::is_null;
589 using Teuchos::SerialDenseMatrix;
590 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
591 typedef typename Array< RCP< const MV > >::size_type size_type;
593 #ifdef BELOS_TEUCHOS_TIME_MONITOR 594 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
597 ScalarType ONE = SCT::one();
598 const MagnitudeType ZERO = MGT::zero();
601 int xc = MVT::GetNumberVecs( X );
602 ptrdiff_t xr = MVT::GetGlobalLength( X );
609 B = rcp (
new serial_dense_matrix_type (xc, xc));
619 for (size_type k = 0; k < nq; ++k)
621 const int numRows = MVT::GetNumberVecs (*Q[k]);
622 const int numCols = xc;
625 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
626 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
628 int err = C[k]->reshape (numRows, numCols);
629 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
630 "IMGS orthogonalization: failed to reshape " 631 "C[" << k <<
"] (the array of block " 632 "coefficients resulting from projecting X " 633 "against Q[1:" << nq <<
"]).");
639 if (MX == Teuchos::null) {
641 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
642 OPT::Apply(*(this->_Op),X,*MX);
647 MX = Teuchos::rcp( &X,
false );
650 int mxc = MVT::GetNumberVecs( *MX );
651 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
654 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::IMGSOrthoManager::projectAndNormalize(): X must be non-empty" );
657 for (
int i=0; i<nq; i++) {
658 numbas += MVT::GetNumberVecs( *Q[i] );
662 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
663 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
665 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
666 "Belos::IMGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
668 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
669 "Belos::IMGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
675 bool dep_flg =
false;
678 Teuchos::RCP<MV> tmpX, tmpMX;
679 tmpX = MVT::CloneCopy(X);
681 tmpMX = MVT::CloneCopy(*MX);
688 dep_flg = blkOrtho1( X, MX, C, Q );
691 if ( B == Teuchos::null ) {
692 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
694 std::vector<ScalarType>
diag(xc);
696 #ifdef BELOS_TEUCHOS_TIME_MONITOR 697 Teuchos::TimeMonitor normTimer( *timerNorm_ );
699 MVT::MvDot( X, *MX,
diag );
701 (*B)(0,0) = SCT::squareroot(SCT::magnitude(
diag[0]));
703 ScalarType scale = ONE;
708 MVT::MvScale( X, scale );
711 MVT::MvScale( *MX, scale );
717 dep_flg = blkOrtho( X, MX, C, Q );
723 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
726 MVT::Assign( *tmpX, X );
728 MVT::Assign( *tmpMX, *MX );
733 rank = findBasis( X, MX, B,
false );
738 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
741 MVT::Assign( *tmpX, X );
743 MVT::Assign( *tmpMX, *MX );
750 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
751 "Belos::IMGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
761 template<
class StorageType,
class MV,
class OP>
762 int IMGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::normalize(
763 MV &X, Teuchos::RCP<MV> MX,
764 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
766 #ifdef BELOS_TEUCHOS_TIME_MONITOR 767 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
771 return findBasis(X, MX, B,
true);
777 template<
class StorageType,
class MV,
class OP>
778 void IMGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::project(
779 MV &X, Teuchos::RCP<MV> MX,
780 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
781 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
797 #ifdef BELOS_TEUCHOS_TIME_MONITOR 798 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
801 int xc = MVT::GetNumberVecs( X );
802 ptrdiff_t xr = MVT::GetGlobalLength( X );
804 std::vector<int> qcs(nq);
806 if (nq == 0 || xc == 0 || xr == 0) {
809 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
818 if (MX == Teuchos::null) {
820 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
821 OPT::Apply(*(this->_Op),X,*MX);
826 MX = Teuchos::rcp( &X,
false );
828 int mxc = MVT::GetNumberVecs( *MX );
829 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
832 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
833 "Belos::IMGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
835 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
836 "Belos::IMGSOrthoManager::project(): Size of X not consistant with MX,Q" );
840 for (
int i=0; i<nq; i++) {
841 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
842 "Belos::IMGSOrthoManager::project(): Q lengths not mutually consistant" );
843 qcs[i] = MVT::GetNumberVecs( *Q[i] );
844 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
845 "Belos::IMGSOrthoManager::project(): Q has less rows than columns" );
849 if ( C[i] == Teuchos::null ) {
850 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
853 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
854 "Belos::IMGSOrthoManager::project(): Size of Q not consistant with size of C" );
859 blkOrtho( X, MX, C, Q );
866 template<
class StorageType,
class MV,
class OP>
867 int IMGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::findBasis(
868 MV &X, Teuchos::RCP<MV> MX,
869 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
870 bool completeBasis,
int howMany )
const {
887 const ScalarType ONE = SCT::one();
888 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
890 int xc = MVT::GetNumberVecs( X );
891 ptrdiff_t xr = MVT::GetGlobalLength( X );
904 if (MX == Teuchos::null) {
906 MX = MVT::Clone(X,xc);
907 OPT::Apply(*(this->_Op),X,*MX);
914 if ( B == Teuchos::null ) {
915 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
918 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
919 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
922 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
923 "Belos::IMGSOrthoManager::findBasis(): X must be non-empty" );
924 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
925 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
926 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
927 "Belos::IMGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
928 TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
929 "Belos::IMGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
930 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
931 "Belos::IMGSOrthoManager::findBasis(): Invalid howMany parameter" );
936 int xstart = xc - howMany;
938 for (
int j = xstart;
j < xc;
j++) {
947 std::vector<int> index(1);
949 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
950 Teuchos::RCP<MV> MXj;
951 if ((this->_hasOp)) {
953 MXj = MVT::CloneViewNonConst( *MX, index );
960 Teuchos::RCP<MV> oldMXj;
962 oldMXj = MVT::CloneCopy( *MXj );
966 Teuchos::SerialDenseVector<int,ScalarType> product(numX);
967 Teuchos::SerialDenseVector<int,ScalarType> P2(1);
968 Teuchos::RCP<const MV> prevX, prevMX;
970 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
975 #ifdef BELOS_TEUCHOS_TIME_MONITOR 976 Teuchos::TimeMonitor normTimer( *timerNorm_ );
978 MVT::MvDot( *Xj, *MXj, oldDot );
981 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
982 "Belos::IMGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
985 for (
int ii=0; ii<numX; ii++) {
988 prevX = MVT::CloneView( X, index );
990 prevMX = MVT::CloneView( *MX, index );
993 for (
int i=0; i<max_ortho_steps_; ++i) {
997 #ifdef BELOS_TEUCHOS_TIME_MONITOR 998 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1000 MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*Xj,MXj,P2);
1006 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1007 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1009 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1017 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1018 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1020 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1025 product[ii] = P2[0];
1027 product[ii] += P2[0];
1035 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1036 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1038 MVT::MvDot( *Xj, *oldMXj, newDot );
1041 newDot[0] = oldDot[0];
1045 if (completeBasis) {
1049 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1054 std::cout <<
"Belos::IMGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1057 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1058 Teuchos::RCP<MV> tempMXj;
1059 MVT::MvRandom( *tempXj );
1061 tempMXj = MVT::Clone( X, 1 );
1062 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1068 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1069 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1071 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_ );
1088 MatOrthoManager<ScalarType,MV,OP>::innerProd(*prevX,*tempXj,tempMXj,P2);
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];
1113 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1114 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1116 MVT::MvDot( *tempXj, *tempMXj, newDot );
1119 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1121 MVT::Assign( *tempXj, *Xj );
1123 MVT::Assign( *tempMXj, *MXj );
1135 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1144 ScalarType
diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1145 ScalarType scale = ONE;
1147 MVT::MvScale( *Xj, scale );
1150 MVT::MvScale( *MXj, scale );
1163 for (
int i=0; i<numX; i++) {
1164 (*B)(i,
j) = product(i);
1175 template<
class StorageType,
class MV,
class OP>
1177 IMGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1178 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1179 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1182 int xc = MVT::GetNumberVecs( X );
1183 const ScalarType ONE = SCT::one();
1185 std::vector<int> qcs( nq );
1186 for (
int i=0; i<nq; i++) {
1187 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1191 std::vector<int> index(1);
1192 Teuchos::RCP<const MV> tempQ;
1194 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1196 for (
int i=0; i<nq; i++) {
1199 for (
int ii=0; ii<qcs[i]; ii++) {
1202 tempQ = MVT::CloneView( *Q[i], index );
1203 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1207 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1208 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1210 MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,X,MX,tempC);
1214 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1215 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1217 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1224 OPT::Apply( *(this->_Op), X, *MX);
1228 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1229 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1230 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1236 for (
int j = 1;
j < max_ortho_steps_; ++
j) {
1238 for (
int i=0; i<nq; i++) {
1240 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],1);
1243 for (
int ii=0; ii<qcs[i]; ii++) {
1246 tempQ = MVT::CloneView( *Q[i], index );
1247 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii, 0 );
1248 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii, 0 );
1252 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1253 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1255 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1259 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1260 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1262 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1271 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1273 else if (xc <= qcs[i]) {
1275 OPT::Apply( *(this->_Op), X, *MX);
1286 template<
class StorageType,
class MV,
class OP>
1288 IMGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1289 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1290 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const 1293 int xc = MVT::GetNumberVecs( X );
1294 bool dep_flg =
false;
1295 const ScalarType ONE = SCT::one();
1297 std::vector<int> qcs( nq );
1298 for (
int i=0; i<nq; i++) {
1299 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1305 std::vector<ScalarType> oldDot( xc );
1307 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1308 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1310 MVT::MvDot( X, *MX, oldDot );
1313 std::vector<int> index(1);
1314 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1315 Teuchos::RCP<const MV> tempQ;
1318 for (
int i=0; i<nq; i++) {
1321 for (
int ii=0; ii<qcs[i]; ii++) {
1324 tempQ = MVT::CloneView( *Q[i], index );
1325 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1329 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1330 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1332 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC);
1336 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1337 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1339 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, X );
1346 OPT::Apply( *(this->_Op), X, *MX);
1350 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1351 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1352 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1358 for (
int j = 1;
j < max_ortho_steps_; ++
j) {
1360 for (
int i=0; i<nq; i++) {
1361 Teuchos::SerialDenseMatrix<int,ScalarType> C2(qcs[i],xc);
1364 for (
int ii=0; ii<qcs[i]; ii++) {
1367 tempQ = MVT::CloneView( *Q[i], index );
1368 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, xc, ii, 0 );
1369 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, xc, ii, 0 );
1373 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1374 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1376 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, X, MX, tempC2 );
1380 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1381 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1383 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, X );
1391 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1393 else if (xc <= qcs[i]) {
1395 OPT::Apply( *(this->_Op), X, *MX);
1402 std::vector<ScalarType> newDot(xc);
1404 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1405 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1407 MVT::MvDot( X, *MX, newDot );
1411 for (
int i=0; i<xc; i++){
1412 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1421 template<
class StorageType,
class MV,
class OP>
1423 IMGSOrthoManager<Sacado::MP::Vector<StorageType>, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1424 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1425 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1426 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const 1428 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1430 const ScalarType ONE = SCT::one();
1431 const ScalarType ZERO = SCT::zero();
1434 int xc = MVT::GetNumberVecs( X );
1435 std::vector<int> indX( 1 );
1436 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1438 std::vector<int> qcs( nq );
1439 for (
int i=0; i<nq; i++) {
1440 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1444 Teuchos::RCP<const MV> lastQ;
1445 Teuchos::RCP<MV> Xj, MXj;
1446 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1449 for (
int j=0;
j<xc;
j++) {
1451 bool dep_flg =
false;
1455 std::vector<int> index(
j );
1456 for (
int ind=0; ind<
j; ind++) {
1459 lastQ = MVT::CloneView( X, index );
1462 Q.push_back( lastQ );
1464 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1469 Xj = MVT::CloneViewNonConst( X, indX );
1471 MXj = MVT::CloneViewNonConst( *MX, indX );
1479 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1480 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1482 MVT::MvDot( *Xj, *MXj, oldDot );
1485 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1486 Teuchos::RCP<const MV> tempQ;
1489 for (
int i=0; i<Q.size(); i++) {
1492 for (
int ii=0; ii<qcs[i]; ii++) {
1495 tempQ = MVT::CloneView( *Q[i], indX );
1497 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], 1, 1, ii,
j );
1500 MatOrthoManager<ScalarType,MV,OP>::innerProd(*tempQ,*Xj,MXj,tempC);
1503 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *Xj );
1509 OPT::Apply( *(this->_Op), *Xj, *MXj);
1513 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1514 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1515 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0,
j );
1516 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1522 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1524 for (
int i=0; i<Q.size(); i++) {
1525 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1528 for (
int ii=0; ii<qcs[i]; ii++) {
1531 tempQ = MVT::CloneView( *Q[i], indX );
1533 Teuchos::SerialDenseMatrix<int,ScalarType> tempC2( Teuchos::View, C2, 1, 1, ii );
1536 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *Xj, MXj, tempC2);
1537 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC2, ONE, *Xj );
1541 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0,
j );
1548 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1550 else if (xc <= qcs[i]) {
1552 OPT::Apply( *(this->_Op), *Xj, *MXj);
1560 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1566 ScalarType
diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1568 MVT::MvScale( *Xj, ONE/
diag );
1571 MVT::MvScale( *MXj, ONE/
diag );
1579 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1580 Teuchos::RCP<MV> tempMXj;
1581 MVT::MvRandom( *tempXj );
1583 tempMXj = MVT::Clone( X, 1 );
1584 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1590 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1591 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1593 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1596 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1598 for (
int i=0; i<Q.size(); i++) {
1599 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1602 for (
int ii=0; ii<qcs[i]; ii++) {
1605 tempQ = MVT::CloneView( *Q[i], indX );
1606 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, product, 1, 1, ii );
1609 MatOrthoManager<ScalarType,MV,OP>::innerProd( *tempQ, *tempXj, tempMXj, tempC );
1610 MVT::MvTimesMatAddMv( -ONE, *tempQ, tempC, ONE, *tempXj );
1618 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1620 else if (xc <= qcs[i]) {
1622 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1630 #ifdef BELOS_TEUCHOS_TIME_MONITOR 1631 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1633 MVT::MvDot( *tempXj, *tempMXj, newDot );
1637 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1638 ScalarType
diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1644 MVT::MvAddMv( ONE/
diag, *tempXj, ZERO, *tempXj, *Xj );
1646 MVT::MvAddMv( ONE/
diag, *tempMXj, ZERO, *tempMXj, *MXj );
1668 #endif // BELOS_IMGS_ORTHOMANAGER_HPP 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 getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
OperatorTraits< ScalarType, MV, OP > OPT
MagnitudeType sing_tol_
Singular block detection threshold.
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
Teuchos::ScalarTraits< ScalarType > SCT
static const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
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.
KOKKOS_INLINE_FUNCTION bool OR(Mask< T > m)
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
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.
KOKKOS_INLINE_FUNCTION MaskedAssign< scalar > mask_assign(bool b, scalar *s)
int max_ortho_steps_
Max number of (re)orthogonalization steps, including the first.
Sacado::MP::Vector< Storage > ScalarType
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.
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
MagnitudeType getSingTol() const
Return parameter for singular block detection.
MultiVecTraits< ScalarType, MV > MVT
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
Teuchos::RCP< Teuchos::ParameterList > defaultParams_
Default parameter list.
~IMGSOrthoManager()
Destructor.
static const int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
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.
std::string label_
Label for timers.
MagnitudeType blk_tol_
Block reorthogonalization tolerance.
Teuchos::ScalarTraits< MagnitudeType > MGT
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...