47#ifndef BELOS_ICGS_ORTHOMANAGER_HPP
48#define BELOS_ICGS_ORTHOMANAGER_HPP
64#include "Teuchos_as.hpp"
65#ifdef BELOS_TEUCHOS_TIME_MONITOR
66#include "Teuchos_TimeMonitor.hpp"
72 template<
class ScalarType,
class MV,
class OP>
76 template<
class ScalarType,
class MV,
class OP>
79 template<
class ScalarType,
class MV,
class OP>
84 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
85 typedef typename Teuchos::ScalarTraits<MagnitudeType> MGT;
86 typedef Teuchos::ScalarTraits<ScalarType> SCT;
96 Teuchos::RCP<const OP> Op = Teuchos::null,
101 max_ortho_steps_( max_ortho_steps ),
103 sing_tol_( sing_tol ),
106#ifdef BELOS_TEUCHOS_TIME_MONITOR
107 std::stringstream ss;
108 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
110 std::string orthoLabel = ss.str() +
": Orthogonalization";
111 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
113 std::string updateLabel = ss.str() +
": Ortho (Update)";
114 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
116 std::string normLabel = ss.str() +
": Ortho (Norm)";
117 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
119 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
120 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
126 const std::string& label =
"Belos",
127 Teuchos::RCP<const OP> Op = Teuchos::null) :
136#ifdef BELOS_TEUCHOS_TIME_MONITOR
137 std::stringstream ss;
138 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
140 std::string orthoLabel = ss.str() +
": Orthogonalization";
141 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
143 std::string updateLabel = ss.str() +
": Ortho (Update)";
144 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
146 std::string normLabel = ss.str() +
": Ortho (Norm)";
147 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
149 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
150 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
164 using Teuchos::Exceptions::InvalidParameterName;
165 using Teuchos::ParameterList;
166 using Teuchos::parameterList;
170 RCP<ParameterList> params;
171 if (plist.is_null()) {
172 params = parameterList (*defaultParams);
187 int maxNumOrthogPasses;
188 MagnitudeType blkTol;
189 MagnitudeType singTol;
192 maxNumOrthogPasses = params->get<
int> (
"maxNumOrthogPasses");
193 }
catch (InvalidParameterName&) {
194 maxNumOrthogPasses = defaultParams->get<
int> (
"maxNumOrthogPasses");
195 params->set (
"maxNumOrthogPasses", maxNumOrthogPasses);
206 blkTol = params->get<MagnitudeType> (
"blkTol");
207 }
catch (InvalidParameterName&) {
209 blkTol = params->get<MagnitudeType> (
"depTol");
212 params->remove (
"depTol");
213 }
catch (InvalidParameterName&) {
214 blkTol = defaultParams->get<MagnitudeType> (
"blkTol");
216 params->set (
"blkTol", blkTol);
220 singTol = params->get<MagnitudeType> (
"singTol");
221 }
catch (InvalidParameterName&) {
222 singTol = defaultParams->get<MagnitudeType> (
"singTol");
223 params->set (
"singTol", singTol);
226 max_ortho_steps_ = maxNumOrthogPasses;
230 this->setMyParamList (params);
233 Teuchos::RCP<const Teuchos::ParameterList>
236 if (defaultParams_.is_null()) {
237 defaultParams_ = Belos::getICGSDefaultParameters<ScalarType, MV, OP>();
240 return defaultParams_;
249 Teuchos::RCP<const Teuchos::ParameterList>
253 using Teuchos::ParameterList;
254 using Teuchos::parameterList;
259 RCP<ParameterList> params = parameterList (*defaultParams);
274 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
275 if (! params.is_null()) {
280 params->set (
"blkTol", blk_tol);
288 Teuchos::RCP<Teuchos::ParameterList> params = this->getNonconstParameterList();
289 if (! params.is_null()) {
294 params->set (
"singTol", sing_tol);
296 sing_tol_ = sing_tol;
338 void project ( MV &X, Teuchos::RCP<MV> MX,
339 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
340 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
346 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
347 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
377 int normalize ( MV &X, Teuchos::RCP<MV> MX,
378 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B)
const;
383 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_; }
509 int max_ortho_steps_;
511 MagnitudeType blk_tol_;
513 MagnitudeType sing_tol_;
517#ifdef BELOS_TEUCHOS_TIME_MONITOR
518 Teuchos::RCP<Teuchos::Time> timerOrtho_, timerUpdate_, timerNorm_, timerInnerProd_;
522 mutable Teuchos::RCP<Teuchos::ParameterList> defaultParams_;
525 int findBasis(MV &X, Teuchos::RCP<MV> MX,
526 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > C,
527 bool completeBasis,
int howMany = -1 )
const;
530 bool blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
531 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
532 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
535 bool blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
536 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
537 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const;
552 int blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
553 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
554 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
555 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const;
559 template<
class ScalarType,
class MV,
class OP>
562 template<
class ScalarType,
class MV,
class OP>
563 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
565 = 10*Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::squareroot(
566 Teuchos::ScalarTraits<
typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps() );
568 template<
class ScalarType,
class MV,
class OP>
569 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
571 = 10*Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::eps();
573 template<
class ScalarType,
class MV,
class OP>
576 template<
class ScalarType,
class MV,
class OP>
577 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
579 = Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
581 template<
class ScalarType,
class MV,
class OP>
582 const typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType
584 = Teuchos::ScalarTraits<typename ICGSOrthoManager<ScalarType,MV,OP>::MagnitudeType>::zero();
588 template<
class ScalarType,
class MV,
class OP>
591 if (label != label_) {
593#ifdef BELOS_TEUCHOS_TIME_MONITOR
594 std::stringstream ss;
595 ss << label_ +
": ICGS[" << max_ortho_steps_ <<
"]";
597 std::string orthoLabel = ss.str() +
": Orthogonalization";
598 timerOrtho_ = Teuchos::TimeMonitor::getNewCounter(orthoLabel);
600 std::string updateLabel = ss.str() +
": Ortho (Update)";
601 timerUpdate_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
603 std::string normLabel = ss.str() +
": Ortho (Norm)";
604 timerNorm_ = Teuchos::TimeMonitor::getNewCounter(normLabel);
606 std::string ipLabel = ss.str() +
": Ortho (Inner Product)";
607 timerInnerProd_ = Teuchos::TimeMonitor::getNewCounter(ipLabel);
614 template<
class ScalarType,
class MV,
class OP>
615 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
617 const ScalarType ONE = SCT::one();
618 int rank = MVT::GetNumberVecs(X);
619 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
621 for (
int i=0; i<rank; i++) {
624 return xTx.normFrobenius();
629 template<
class ScalarType,
class MV,
class OP>
630 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
632 int r1 = MVT::GetNumberVecs(X1);
633 int r2 = MVT::GetNumberVecs(X2);
634 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r2,r1);
636 return xTx.normFrobenius();
641 template<
class ScalarType,
class MV,
class OP>
646 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
647 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
648 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
650 using Teuchos::Array;
652 using Teuchos::is_null;
655 using Teuchos::SerialDenseMatrix;
656 typedef SerialDenseMatrix< int, ScalarType > serial_dense_matrix_type;
657 typedef typename Array< RCP< const MV > >::size_type size_type;
659#ifdef BELOS_TEUCHOS_TIME_MONITOR
660 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
663 ScalarType ONE = SCT::one();
664 const MagnitudeType ZERO = MGT::zero();
667 int xc = MVT::GetNumberVecs( X );
668 ptrdiff_t xr = MVT::GetGlobalLength( X );
675 B = rcp (
new serial_dense_matrix_type (xc, xc));
685 for (size_type k = 0; k < nq; ++k)
687 const int numRows = MVT::GetNumberVecs (*Q[k]);
688 const int numCols = xc;
691 C[k] = rcp (
new serial_dense_matrix_type (numRows, numCols));
692 else if (C[k]->numRows() != numRows || C[k]->numCols() != numCols)
694 int err = C[k]->reshape (numRows, numCols);
695 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, std::runtime_error,
696 "IMGS orthogonalization: failed to reshape "
697 "C[" << k <<
"] (the array of block "
698 "coefficients resulting from projecting X "
699 "against Q[1:" << nq <<
"]).");
705 if (MX == Teuchos::null) {
707 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
708 OPT::Apply(*(this->_Op),X,*MX);
713 MX = Teuchos::rcp( &X,
false );
716 int mxc = MVT::GetNumberVecs( *MX );
717 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
720 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Belos::ICGSOrthoManager::projectAndNormalize(): X must be non-empty" );
723 for (
int i=0; i<nq; i++) {
724 numbas += MVT::GetNumberVecs( *Q[i] );
728 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
729 "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of B" );
731 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
732 "Belos::ICGSOrthoManager::projectAndNormalize(): MVT returned negative dimensions for X,MX" );
734 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
735 "Belos::ICGSOrthoManager::projectAndNormalize(): Size of X must be consistant with size of MX" );
741 bool dep_flg =
false;
747 dep_flg = blkOrtho1( X, MX, C, Q );
750 if ( B == Teuchos::null ) {
751 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
753 std::vector<ScalarType> diag(xc);
755#ifdef BELOS_TEUCHOS_TIME_MONITOR
756 Teuchos::TimeMonitor normTimer( *timerNorm_ );
758 MVT::MvDot( X, *MX, diag );
760 (*B)(0,0) = SCT::squareroot(SCT::magnitude(diag[0]));
762 if (SCT::magnitude((*B)(0,0)) > ZERO) {
764 MVT::MvScale( X, ONE/(*B)(0,0) );
767 MVT::MvScale( *MX, ONE/(*B)(0,0) );
774 Teuchos::RCP<MV> tmpX, tmpMX;
775 tmpX = MVT::CloneCopy(X);
777 tmpMX = MVT::CloneCopy(*MX);
781 dep_flg = blkOrtho( X, MX, C, Q );
787 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
790 MVT::Assign( *tmpX, X );
792 MVT::Assign( *tmpMX, *MX );
797 rank = findBasis( X, MX, B,
false );
802 rank = blkOrthoSing( *tmpX, tmpMX, C, B, Q );
805 MVT::Assign( *tmpX, X );
807 MVT::Assign( *tmpMX, *MX );
814 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
815 "Belos::ICGSOrthoManager::projectAndNormalize(): Debug error in rank variable." );
825 template<
class ScalarType,
class MV,
class OP>
827 MV &X, Teuchos::RCP<MV> MX,
828 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B )
const {
830#ifdef BELOS_TEUCHOS_TIME_MONITOR
831 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
835 return findBasis(X, MX, B,
true);
842 template<
class ScalarType,
class MV,
class OP>
844 MV &X, Teuchos::RCP<MV> MX,
845 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
846 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const {
862#ifdef BELOS_TEUCHOS_TIME_MONITOR
863 Teuchos::TimeMonitor orthotimer(*timerOrtho_);
866 int xc = MVT::GetNumberVecs( X );
867 ptrdiff_t xr = MVT::GetGlobalLength( X );
869 std::vector<int> qcs(nq);
871 if (nq == 0 || xc == 0 || xr == 0) {
874 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
883 if (MX == Teuchos::null) {
885 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
886 OPT::Apply(*(this->_Op),X,*MX);
891 MX = Teuchos::rcp( &X,
false );
893 int mxc = MVT::GetNumberVecs( *MX );
894 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
897 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
898 "Belos::ICGSOrthoManager::project(): MVT returned negative dimensions for X,MX" );
900 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
901 "Belos::ICGSOrthoManager::project(): Size of X not consistant with MX,Q" );
904 for (
int i=0; i<nq; i++) {
905 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
906 "Belos::ICGSOrthoManager::project(): Q lengths not mutually consistant" );
907 qcs[i] = MVT::GetNumberVecs( *Q[i] );
908 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
909 "Belos::ICGSOrthoManager::project(): Q has less rows than columns" );
912 if ( C[i] == Teuchos::null ) {
913 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
916 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
917 "Belos::ICGSOrthoManager::project(): Size of Q not consistant with size of C" );
922 blkOrtho( X, MX, C, Q );
929 template<
class ScalarType,
class MV,
class OP>
934 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
952 const ScalarType ONE = SCT::one ();
953 const MagnitudeType ZERO = SCT::magnitude (SCT::zero ());
955 const int xc = MVT::GetNumberVecs (X);
956 const ptrdiff_t xr = MVT::GetGlobalLength (X);
969 if (MX == Teuchos::null) {
971 MX = MVT::Clone(X,xc);
972 OPT::Apply(*(this->_Op),X,*MX);
979 if ( B == Teuchos::null ) {
980 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
983 const int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
984 const ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
987 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
988 "Belos::ICGSOrthoManager::findBasis(): X must be non-empty" );
989 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
990 "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of B" );
991 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
992 "Belos::ICGSOrthoManager::findBasis(): Size of X not consistant with size of MX" );
993 TEUCHOS_TEST_FOR_EXCEPTION( as<ptrdiff_t> (xc) > xr, std::invalid_argument,
994 "Belos::ICGSOrthoManager::findBasis(): Size of X not feasible for normalization" );
995 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::invalid_argument,
996 "Belos::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
1001 int xstart = xc - howMany;
1003 for (
int j = xstart; j < xc; j++) {
1009 bool addVec =
false;
1012 std::vector<int> index(1);
1014 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1015 Teuchos::RCP<MV> MXj;
1018 MXj = MVT::CloneViewNonConst( *MX, index );
1026 std::vector<int> prev_idx( numX );
1027 Teuchos::RCP<const MV> prevX, prevMX;
1028 Teuchos::RCP<MV> oldMXj;
1031 for (
int i=0; i<numX; i++) {
1034 prevX = MVT::CloneView( X, prev_idx );
1036 prevMX = MVT::CloneView( *MX, prev_idx );
1039 oldMXj = MVT::CloneCopy( *MXj );
1043 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1044 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1049#ifdef BELOS_TEUCHOS_TIME_MONITOR
1050 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1052 MVT::MvDot( *Xj, *MXj, oldDot );
1055 TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(oldDot[0]) < ZERO, OrthoError,
1056 "Belos::ICGSOrthoManager::findBasis(): Negative definiteness discovered in inner product" );
1060 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1062 for (
int i=0; i<max_ortho_steps_; ++i) {
1066#ifdef BELOS_TEUCHOS_TIME_MONITOR
1067 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1075#ifdef BELOS_TEUCHOS_TIME_MONITOR
1076 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1078 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1086#ifdef BELOS_TEUCHOS_TIME_MONITOR
1087 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1089 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1103#ifdef BELOS_TEUCHOS_TIME_MONITOR
1104 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1106 MVT::MvDot( *Xj, *oldMXj, newDot );
1109 newDot[0] = oldDot[0];
1113 if (completeBasis) {
1117 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(sing_tol_*oldDot[0]) ) {
1122 std::cout <<
"Belos::ICGSOrthoManager::findBasis() --> Random for column " << numX << std::endl;
1125 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1126 Teuchos::RCP<MV> tempMXj;
1127 MVT::MvRandom( *tempXj );
1129 tempMXj = MVT::Clone( X, 1 );
1130 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1136#ifdef BELOS_TEUCHOS_TIME_MONITOR
1137 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1139 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1142 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++){
1144#ifdef BELOS_TEUCHOS_TIME_MONITOR
1145 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1150#ifdef BELOS_TEUCHOS_TIME_MONITOR
1151 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1153 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *tempXj );
1156#ifdef BELOS_TEUCHOS_TIME_MONITOR
1157 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1159 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *tempMXj );
1164#ifdef BELOS_TEUCHOS_TIME_MONITOR
1165 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1167 MVT::MvDot( *tempXj, *tempMXj, newDot );
1170 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1172 MVT::Assign( *tempXj, *Xj );
1174 MVT::Assign( *tempMXj, *MXj );
1186 if ( SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*blk_tol_) ) {
1194 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1195 if (SCT::magnitude(diag) > ZERO) {
1196 MVT::MvScale( *Xj, ONE/diag );
1199 MVT::MvScale( *MXj, ONE/diag );
1213 for (
int i=0; i<numX; i++) {
1214 (*B)(i,j) = product(i,0);
1225 template<
class ScalarType,
class MV,
class OP>
1227 ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho1 ( MV &X, Teuchos::RCP<MV> MX,
1228 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1229 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
1232 int xc = MVT::GetNumberVecs( X );
1233 const ScalarType ONE = SCT::one();
1235 std::vector<int> qcs( nq );
1236 for (
int i=0; i<nq; i++) {
1237 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1242 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1244 for (
int i=0; i<nq; i++) {
1247#ifdef BELOS_TEUCHOS_TIME_MONITOR
1248 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1254#ifdef BELOS_TEUCHOS_TIME_MONITOR
1255 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1257 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1263 OPT::Apply( *(this->_Op), X, *MX);
1267 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1268 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1270#ifdef BELOS_TEUCHOS_TIME_MONITOR
1271 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1273 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1280 for (
int j = 1; j < max_ortho_steps_; ++j) {
1282 for (
int i=0; i<nq; i++) {
1283 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1287#ifdef BELOS_TEUCHOS_TIME_MONITOR
1288 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1294#ifdef BELOS_TEUCHOS_TIME_MONITOR
1295 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1297 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1303#ifdef BELOS_TEUCHOS_TIME_MONITOR
1304 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1307 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1309 else if (xc <= qcs[i]) {
1311 OPT::Apply( *(this->_Op), X, *MX);
1322 template<
class ScalarType,
class MV,
class OP>
1324 ICGSOrthoManager<ScalarType, MV, OP>::blkOrtho ( MV &X, Teuchos::RCP<MV> MX,
1325 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1326 Teuchos::ArrayView<Teuchos::RCP<const MV> > Q)
const
1329 int xc = MVT::GetNumberVecs( X );
1330 bool dep_flg =
false;
1331 const ScalarType ONE = SCT::one();
1333 std::vector<int> qcs( nq );
1334 for (
int i=0; i<nq; i++) {
1335 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1341 std::vector<ScalarType> oldDot( xc );
1343#ifdef BELOS_TEUCHOS_TIME_MONITOR
1344 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1346 MVT::MvDot( X, *MX, oldDot );
1349 Teuchos::Array<Teuchos::RCP<MV> > MQ(nq);
1351 for (
int i=0; i<nq; i++) {
1354#ifdef BELOS_TEUCHOS_TIME_MONITOR
1355 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1361#ifdef BELOS_TEUCHOS_TIME_MONITOR
1362 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1364 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
1369 OPT::Apply( *(this->_Op), X, *MX);
1373 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1374 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1376#ifdef BELOS_TEUCHOS_TIME_MONITOR
1377 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1379 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
1386 for (
int j = 1; j < max_ortho_steps_; ++j) {
1388 for (
int i=0; i<nq; i++) {
1389 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(),C[i]->numCols());
1393#ifdef BELOS_TEUCHOS_TIME_MONITOR
1394 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1400#ifdef BELOS_TEUCHOS_TIME_MONITOR
1401 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1403 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
1409#ifdef BELOS_TEUCHOS_TIME_MONITOR
1410 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1413 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
1415 else if (xc <= qcs[i]) {
1417 OPT::Apply( *(this->_Op), X, *MX);
1424 std::vector<ScalarType> newDot(xc);
1426#ifdef BELOS_TEUCHOS_TIME_MONITOR
1427 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1429 MVT::MvDot( X, *MX, newDot );
1433 for (
int i=0; i<xc; i++){
1434 if (SCT::magnitude(newDot[i]) < SCT::magnitude(oldDot[i] * blk_tol_)) {
1443 template<
class ScalarType,
class MV,
class OP>
1445 ICGSOrthoManager<ScalarType, MV, OP>::blkOrthoSing ( MV &X, Teuchos::RCP<MV> MX,
1446 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
1447 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
1448 Teuchos::ArrayView<Teuchos::RCP<const MV> > QQ)
const
1450 Teuchos::Array<Teuchos::RCP<const MV> > Q (QQ);
1452 const ScalarType ONE = SCT::one();
1453 const ScalarType ZERO = SCT::zero();
1456 int xc = MVT::GetNumberVecs( X );
1457 std::vector<int> indX( 1 );
1458 std::vector<ScalarType> oldDot( 1 ), newDot( 1 );
1460 std::vector<int> qcs( nq );
1461 for (
int i=0; i<nq; i++) {
1462 qcs[i] = MVT::GetNumberVecs( *Q[i] );
1466 Teuchos::RCP<const MV> lastQ;
1467 Teuchos::RCP<MV> Xj, MXj;
1468 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lastC;
1471 for (
int j=0; j<xc; j++) {
1473 bool dep_flg =
false;
1477 std::vector<int> index( j );
1478 for (
int ind=0; ind<j; ind++) {
1481 lastQ = MVT::CloneView( X, index );
1484 Q.push_back( lastQ );
1486 qcs.push_back( MVT::GetNumberVecs( *lastQ ) );
1491 Xj = MVT::CloneViewNonConst( X, indX );
1493 MXj = MVT::CloneViewNonConst( *MX, indX );
1501#ifdef BELOS_TEUCHOS_TIME_MONITOR
1502 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1504 MVT::MvDot( *Xj, *MXj, oldDot );
1507 Teuchos::Array<Teuchos::RCP<MV> > MQ(Q.size());
1509 for (
int i=0; i<Q.size(); i++) {
1512 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1516#ifdef BELOS_TEUCHOS_TIME_MONITOR
1517 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1522#ifdef BELOS_TEUCHOS_TIME_MONITOR
1523 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1526 MVT::MvTimesMatAddMv( -ONE, *Q[i], tempC, ONE, *Xj );
1531 OPT::Apply( *(this->_Op), *Xj, *MXj);
1535 MQ[i] = MVT::Clone( *Q[i], qcs[i] );
1536 OPT::Apply( *(this->_Op), *Q[i], *MQ[i] );
1538#ifdef BELOS_TEUCHOS_TIME_MONITOR
1539 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1541 MVT::MvTimesMatAddMv( -ONE, *MQ[i], tempC, ONE, *MXj );
1548 for (
int num_ortho_steps=1; num_ortho_steps < max_ortho_steps_; ++num_ortho_steps) {
1550 for (
int i=0; i<Q.size(); i++) {
1551 Teuchos::SerialDenseMatrix<int,ScalarType> tempC( Teuchos::View, *C[i], qcs[i], 1, 0, j );
1552 Teuchos::SerialDenseMatrix<int,ScalarType> C2( qcs[i], 1 );
1556#ifdef BELOS_TEUCHOS_TIME_MONITOR
1557 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1563#ifdef BELOS_TEUCHOS_TIME_MONITOR
1564 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1566 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, *Xj );
1573#ifdef BELOS_TEUCHOS_TIME_MONITOR
1574 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1576 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MXj );
1578 else if (xc <= qcs[i]) {
1580 OPT::Apply( *(this->_Op), *Xj, *MXj);
1589#ifdef BELOS_TEUCHOS_TIME_MONITOR
1590 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1592 MVT::MvDot( *Xj, *MXj, newDot );
1596 if (SCT::magnitude(newDot[0]) < SCT::magnitude(oldDot[0]*sing_tol_)) {
1602 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1604 MVT::MvScale( *Xj, ONE/diag );
1607 MVT::MvScale( *MXj, ONE/diag );
1615 Teuchos::RCP<MV> tempXj = MVT::Clone( X, 1 );
1616 Teuchos::RCP<MV> tempMXj;
1617 MVT::MvRandom( *tempXj );
1619 tempMXj = MVT::Clone( X, 1 );
1620 OPT::Apply( *(this->_Op), *tempXj, *tempMXj );
1626#ifdef BELOS_TEUCHOS_TIME_MONITOR
1627 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1629 MVT::MvDot( *tempXj, *tempMXj, oldDot );
1632 for (
int num_orth=0; num_orth<max_ortho_steps_; num_orth++) {
1634 for (
int i=0; i<Q.size(); i++) {
1635 Teuchos::SerialDenseMatrix<int,ScalarType> product( qcs[i], 1 );
1639#ifdef BELOS_TEUCHOS_TIME_MONITOR
1640 Teuchos::TimeMonitor innerProdTimer( *timerInnerProd_ );
1645#ifdef BELOS_TEUCHOS_TIME_MONITOR
1646 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1648 MVT::MvTimesMatAddMv( -ONE, *Q[i], product, ONE, *tempXj );
1654#ifdef BELOS_TEUCHOS_TIME_MONITOR
1655 Teuchos::TimeMonitor updateTimer( *timerUpdate_ );
1658 MVT::MvTimesMatAddMv( -ONE, *MQ[i], product, ONE, *tempMXj );
1660 else if (xc <= qcs[i]) {
1662 OPT::Apply( *(this->_Op), *tempXj, *tempMXj);
1671#ifdef BELOS_TEUCHOS_TIME_MONITOR
1672 Teuchos::TimeMonitor normTimer( *timerNorm_ );
1674 MVT::MvDot( *tempXj, *tempMXj, newDot );
1678 if ( SCT::magnitude(newDot[0]) >= SCT::magnitude(oldDot[0]*sing_tol_) ) {
1679 ScalarType diag = SCT::squareroot(SCT::magnitude(newDot[0]));
1685 MVT::MvAddMv( ONE/diag, *tempXj, ZERO, *tempXj, *Xj );
1687 MVT::MvAddMv( ONE/diag, *tempMXj, ZERO, *tempMXj, *MXj );
1707 template<
class ScalarType,
class MV,
class OP>
1710 using Teuchos::ParameterList;
1711 using Teuchos::parameterList;
1714 RCP<ParameterList> params = parameterList (
"ICGS");
1719 "Maximum number of orthogonalization passes (includes the "
1720 "first). Default is 2, since \"twice is enough\" for Krylov "
1723 "Block reorthogonalization threshold.");
1725 "Singular block detection threshold.");
1730 template<
class ScalarType,
class MV,
class OP>
1733 using Teuchos::ParameterList;
1736 RCP<ParameterList> params = getICGSDefaultParameters<ScalarType, MV, OP>();
1738 params->set (
"maxNumOrthogPasses",
1740 params->set (
"blkTol",
1742 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 getSingTol() const
Return parameter for singular block detection.
MagnitudeType getBlkTol() const
Return parameter for block re-orthogonalization threshhold.
~ICGSOrthoManager()
Destructor.
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 const MagnitudeType sing_tol_fast_
Singular block detection threshold (fast).
static const MagnitudeType blk_tol_fast_
Block reorthogonalization threshold (fast).
void project(MV &X, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C, Teuchos::ArrayView< Teuchos::RCP< const MV > > Q) const
This method calls project(X,Teuchos::null,C,Q); see documentation for that function.
void 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 int max_ortho_steps_default_
Max number of (re)orthogonalization steps, including the first (default).
void setBlkTol(const MagnitudeType blk_tol)
Set parameter for block re-orthogonalization threshhold.
void setParameterList(const Teuchos::RCP< Teuchos::ParameterList > &plist)
ICGSOrthoManager(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.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
void setSingTol(const MagnitudeType sing_tol)
Set parameter for singular block detection.
static const MagnitudeType sing_tol_default_
Singular block detection threshold (default).
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 ...
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 ,...
const std::string & getLabel() const
This method returns the label being used by the timers in the orthogonalization manager.
static const MagnitudeType blk_tol_default_
Block reorthogonalization threshold (default).
ICGSOrthoManager(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.
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...
Teuchos::RCP< const Teuchos::ParameterList > getFastParameters() const
"Fast" but possibly unsafe or less accurate parameters.
void setLabel(const std::string &label)
This method sets the label used by the timers in the orthogonalization manager.
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 .
static const int max_ortho_steps_fast_
Max number of (re)orthogonalization steps, including the first (fast).
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.
Teuchos::RCP< Teuchos::ParameterList > getICGSDefaultParameters()
"Default" parameters for robustness and accuracy.
Teuchos::RCP< Teuchos::ParameterList > getICGSFastParameters()
"Fast" but possibly unsafe or less accurate parameters.