47 #ifndef ANASAZI_ICSG_ORTHOMANAGER_HPP 48 #define ANASAZI_ICSG_ORTHOMANAGER_HPP 61 #include "Teuchos_TimeMonitor.hpp" 62 #include "Teuchos_LAPACK.hpp" 63 #include "Teuchos_BLAS.hpp" 64 #ifdef ANASAZI_ICGS_DEBUG 65 # include <Teuchos_FancyOStream.hpp> 70 template<
class ScalarType,
class MV,
class OP>
74 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
75 typedef Teuchos::ScalarTraits<ScalarType> SCT;
83 ICGSOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null,
int numIters = 2,
85 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
86 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
180 Teuchos::Array<Teuchos::RCP<const MV> > X,
181 Teuchos::Array<Teuchos::RCP<const MV> > Y,
183 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
184 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
185 Teuchos::RCP<MV> MS = Teuchos::null,
186 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
187 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
284 Teuchos::Array<Teuchos::RCP<const MV> > X,
285 Teuchos::Array<Teuchos::RCP<const MV> > Y,
287 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
288 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
289 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
290 Teuchos::RCP<MV> MS = Teuchos::null,
291 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
292 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
316 Teuchos::Array<Teuchos::RCP<const MV> > Q,
317 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
318 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
319 Teuchos::RCP<MV> MX = Teuchos::null,
320 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
334 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
335 Teuchos::RCP<MV> MX = Teuchos::null
349 Teuchos::Array<Teuchos::RCP<const MV> > Q,
350 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
351 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
352 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
353 Teuchos::RCP<MV> MX = Teuchos::null,
354 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
366 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
373 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
374 orthogErrorMat(
const MV &X1,
const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2)
const;
383 numIters_ = numIters;
384 TEUCHOS_TEST_FOR_EXCEPTION(numIters_ < 1,std::invalid_argument,
385 "Anasazi::ICGSOrthoManager::setNumIters(): input must be >= 1.");
401 int findBasis(MV &X, Teuchos::RCP<MV> MX,
402 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
403 bool completeBasis,
int howMany = -1)
const;
410 template<
class ScalarType,
class MV,
class OP>
413 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
414 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol) :
418 TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
419 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative.");
421 Teuchos::LAPACK<int,MagnitudeType> lapack;
422 eps_ = lapack.LAMCH(
'E');
423 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.50);
425 TEUCHOS_TEST_FOR_EXCEPTION(
426 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
427 std::invalid_argument,
428 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1].");
435 template<
class ScalarType,
class MV,
class OP>
436 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
438 const ScalarType ONE = SCT::one();
439 int rank = MVT::GetNumberVecs(X);
440 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
442 for (
int i=0; i<rank; i++) {
445 return xTx.normFrobenius();
452 template<
class ScalarType,
class MV,
class OP>
453 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
455 int r1 = MVT::GetNumberVecs(X1);
456 int r2 = MVT::GetNumberVecs(X2);
457 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
459 return xTx.normFrobenius();
465 template<
class ScalarType,
class MV,
class OP>
468 Teuchos::Array<Teuchos::RCP<const MV> > Q,
469 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
471 Teuchos::Array<Teuchos::RCP<const MV> > MQ
474 projectGen(X,Q,Q,
true,C,MX,MQ,MQ);
481 template<
class ScalarType,
class MV,
class OP>
484 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
485 Teuchos::RCP<MV> MX)
const 490 int xc = MVT::GetNumberVecs(X);
491 ptrdiff_t xr = MVT::GetGlobalLength(X);
496 if (MX == Teuchos::null) {
498 MX = MVT::Clone(X,xc);
499 OPT::Apply(*(this->_Op),X,*MX);
500 this->_OpCounter += MVT::GetNumberVecs(X);
506 if ( B == Teuchos::null ) {
507 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
510 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
511 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
514 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
515 "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" );
516 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
517 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
518 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
519 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
520 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
521 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
523 return findBasis(X, MX, *B,
true );
530 template<
class ScalarType,
class MV,
class OP>
533 Teuchos::Array<Teuchos::RCP<const MV> > Q,
534 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
535 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
537 Teuchos::Array<Teuchos::RCP<const MV> > MQ
540 return projectAndNormalizeGen(X,Q,Q,
true,C,B,MX,MQ,MQ);
546 template<
class ScalarType,
class MV,
class OP>
549 Teuchos::Array<Teuchos::RCP<const MV> > X,
550 Teuchos::Array<Teuchos::RCP<const MV> > Y,
552 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
554 Teuchos::Array<Teuchos::RCP<const MV> > MX,
555 Teuchos::Array<Teuchos::RCP<const MV> > MY
573 #ifdef ANASAZI_ICGS_DEBUG 575 Teuchos::RCP<Teuchos::FancyOStream>
576 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
577 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
578 *out <<
"Entering Anasazi::ICGSOrthoManager::projectGen(...)\n";
581 const ScalarType ONE = SCT::one();
582 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
583 Teuchos::LAPACK<int,ScalarType> lapack;
584 Teuchos::BLAS<int,ScalarType> blas;
586 int sc = MVT::GetNumberVecs( S );
587 ptrdiff_t sr = MVT::GetGlobalLength( S );
588 int numxy = X.length();
589 TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument,
590 "Anasazi::ICGSOrthoManager::projectGen(): X and Y must contain the same number of multivectors.");
591 std::vector<int> xyc(numxy);
593 if (numxy == 0 || sc == 0 || sr == 0) {
594 #ifdef ANASAZI_ICGS_DEBUG 595 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
609 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument,
610 "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." );
613 if (this->_hasOp ==
true) {
614 if (MS != Teuchos::null) {
615 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MS) != sr, std::invalid_argument,
616 "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." );
617 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
618 "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." );
623 ptrdiff_t sumxyc = 0;
626 for (
int i=0; i<numxy; i++) {
627 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
628 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*X[i]) != sr, std::invalid_argument,
629 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"] length not consistent with S." );
630 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*Y[i]) != sr, std::invalid_argument,
631 "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i <<
"] length not consistent with S." );
632 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
633 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"] and Y[" << i <<
"] widths not consistent." );
635 xyc[i] = MVT::GetNumberVecs( *X[i] );
636 TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument,
637 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"],Y[" << i <<
"] have less rows than columns, and therefore cannot be full rank." );
641 if ( C[i] == Teuchos::null ) {
642 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
645 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument,
646 "Anasazi::ICGSOrthoManager::projectGen(): Size of Q not consistent with size of C." );
650 if (MX[i] != Teuchos::null) {
651 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
652 "Anasazi::ICGSOrthoManager::projectGen(): Size of MX[" << i <<
"] not consistent with size of X[" << i <<
"]." );
657 if (MY[i] != Teuchos::null) {
658 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
659 "Anasazi::ICGSOrthoManager::projectGen(): Size of MY[" << i <<
"] not consistent with size of Y[" << i <<
"]." );
667 TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
668 "Anasazi::ICGSOrthoManager::projectGen(): " 669 << (X[i] == Teuchos::null ?
"Y[" :
"X[") << i <<
"] was provided but " 670 << (X[i] == Teuchos::null ?
"X[" :
"Y[") << i <<
"] was not.");
675 TEUCHOS_TEST_FOR_EXCEPTION(sumxyc > sr, std::invalid_argument,
676 "Anasazi::ICGSOrthoManager::projectGen(): dimension of all X[i],Y[i] is " 677 << sumxyc <<
", but length of vectors is only " << sr <<
". This is infeasible.");
712 if (MXmissing == 0) {
715 if (MYmissing == 0) {
722 switch (whichAlloc) {
738 if (MS == Teuchos::null) {
739 #ifdef ANASAZI_ICGS_DEBUG 740 *out <<
"Allocating MS...\n";
742 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
743 OPT::Apply(*(this->_Op),S,*MS);
744 this->_OpCounter += MVT::GetNumberVecs(S);
748 if (whichAlloc == 0) {
750 for (
int i=0; i<numxy; i++) {
751 if (MX[i] == Teuchos::null) {
752 #ifdef ANASAZI_ICGS_DEBUG 753 *out <<
"Allocating MX[" << i <<
"]...\n";
755 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
756 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
758 this->_OpCounter += xyc[i];
765 MS = Teuchos::rcpFromRef(S);
768 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
769 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
778 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > YMX(numxy);
779 if (isBiortho ==
false) {
780 for (
int i=0; i<numxy; i++) {
781 #ifdef ANASAZI_ICGS_DEBUG 782 *out <<
"Computing YMX[" << i <<
"] and its Cholesky factorization...\n";
784 YMX[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],xyc[i]) );
786 #ifdef ANASAZI_ICGS_DEBUG 793 MagnitudeType err = ZERO;
794 for (
int jj=0; jj<YMX[i]->numCols(); jj++) {
795 err =+ SCT::magnitude(SCT::imag((*YMX[i])(jj,jj)));
796 for (
int ii=jj; ii<YMX[i]->numRows(); ii++) {
797 err += SCT::magnitude( (*YMX[i])(ii,jj) - SCT::conjugate((*YMX[i])(jj,ii)) );
800 *out <<
"Symmetry error in YMX[" << i <<
"] == " << err <<
"\n";
805 lapack.POTRF(
'U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info);
806 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
807 "Anasazi::ICGSOrthoManager::projectGen(): Error computing Cholesky factorization of Y[i]^T * M * X[i] using POTRF: returned info " << info);
812 #ifdef ANASAZI_ICGS_DEBUG 813 std::vector<MagnitudeType> oldNorms(sc);
815 *out <<
"oldNorms = { ";
816 std::copy(oldNorms.begin(), oldNorms.end(), std::ostream_iterator<MagnitudeType>(*out,
" "));
822 Teuchos::Array<Teuchos::SerialDenseMatrix<int,ScalarType> > Ccur(numxy);
823 for (
int i=0; i<numxy; i++) {
824 C[i]->putScalar(ZERO);
825 Ccur[i].reshape(C[i]->numRows(),C[i]->numCols());
828 for (
int iter=0; iter < numIters_; iter++) {
829 #ifdef ANASAZI_ICGS_DEBUG 830 *out <<
"beginning iteration " << iter+1 <<
"\n";
834 for (
int i=0; i<numxy; i++) {
837 if (isBiortho ==
false) {
840 lapack.POTRS(
'U',YMX[i]->numCols(),Ccur[i].numCols(),
841 YMX[i]->values(),YMX[i]->stride(),
842 Ccur[i].values(),Ccur[i].stride(), &info);
843 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
844 "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info <<
" from lapack::POTRS." );
848 #ifdef ANASAZI_ICGS_DEBUG 849 *out <<
"Applying projector P_{X[" << i <<
"],Y[" << i <<
"]}...\n";
851 MVT::MvTimesMatAddMv( -ONE, *X[i], Ccur[i], ONE, S );
858 #ifdef ANASAZI_ICGS_DEBUG 859 *out <<
"Updating MS...\n";
861 OPT::Apply( *(this->_Op), S, *MS);
862 this->_OpCounter += sc;
864 else if (updateMS == 2) {
865 #ifdef ANASAZI_ICGS_DEBUG 866 *out <<
"Updating MS...\n";
868 MVT::MvTimesMatAddMv( -ONE, *MX[i], Ccur[i], ONE, *MS );
873 #ifdef ANASAZI_ICGS_DEBUG 874 std::vector<MagnitudeType> newNorms(sc);
876 *out <<
"newNorms = { ";
877 std::copy(newNorms.begin(), newNorms.end(), std::ostream_iterator<MagnitudeType>(*out,
" "));
882 #ifdef ANASAZI_ICGS_DEBUG 883 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
891 template<
class ScalarType,
class MV,
class OP>
894 Teuchos::Array<Teuchos::RCP<const MV> > X,
895 Teuchos::Array<Teuchos::RCP<const MV> > Y,
897 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
898 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
900 Teuchos::Array<Teuchos::RCP<const MV> > MX,
901 Teuchos::Array<Teuchos::RCP<const MV> > MY
919 #ifdef ANASAZI_ICGS_DEBUG 921 Teuchos::RCP<Teuchos::FancyOStream>
922 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
923 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
924 *out <<
"Entering Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
928 int sc = MVT::GetNumberVecs( S );
929 ptrdiff_t sr = MVT::GetGlobalLength( S );
930 int numxy = X.length();
931 TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument,
932 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X and Y must contain the same number of multivectors.");
933 std::vector<int> xyc(numxy);
935 if (sc == 0 || sr == 0) {
936 #ifdef ANASAZI_ICGS_DEBUG 937 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
951 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument,
952 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." );
955 if (this->_hasOp ==
true) {
956 if (MS != Teuchos::null) {
957 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MS) != sr, std::invalid_argument,
958 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." );
959 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
960 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." );
965 ptrdiff_t sumxyc = 0;
968 for (
int i=0; i<numxy; i++) {
969 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
970 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*X[i]) != sr, std::invalid_argument,
971 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"] length not consistent with S." );
972 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*Y[i]) != sr, std::invalid_argument,
973 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i <<
"] length not consistent with S." );
974 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
975 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"] and Y[" << i <<
"] widths not consistent." );
977 xyc[i] = MVT::GetNumberVecs( *X[i] );
978 TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument,
979 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"],Y[" << i <<
"] have less rows than columns, and therefore cannot be full rank." );
983 if ( C[i] == Teuchos::null ) {
984 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
987 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument,
988 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of Q not consistent with size of C." );
992 if (MX[i] != Teuchos::null) {
993 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
994 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MX[" << i <<
"] not consistent with size of X[" << i <<
"]." );
999 if (MY[i] != Teuchos::null) {
1000 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
1001 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MY[" << i <<
"] not consistent with size of Y[" << i <<
"]." );
1004 MYmissing += xyc[i];
1009 TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
1010 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): " 1011 << (X[i] == Teuchos::null ?
"Y[" :
"X[") << i <<
"] was provided but " 1012 << (X[i] == Teuchos::null ?
"X[" :
"Y[") << i <<
"] was not.");
1017 TEUCHOS_TEST_FOR_EXCEPTION(sumxyc + sc > sr, std::invalid_argument,
1018 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): dimension of all X[i],Y[i] is " 1019 << sumxyc <<
" and requested " << sc <<
"-dimensional basis, but length of vectors is only " 1020 << sr <<
". This is infeasible.");
1055 if (MXmissing == 0) {
1058 if (MYmissing == 0) {
1065 switch (whichAlloc) {
1081 if (MS == Teuchos::null) {
1082 #ifdef ANASAZI_ICGS_DEBUG 1083 *out <<
"Allocating MS...\n";
1085 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
1086 OPT::Apply(*(this->_Op),S,*MS);
1087 this->_OpCounter += MVT::GetNumberVecs(S);
1091 if (whichAlloc == 0) {
1093 for (
int i=0; i<numxy; i++) {
1094 if (MX[i] == Teuchos::null) {
1095 #ifdef ANASAZI_ICGS_DEBUG 1096 *out <<
"Allocating MX[" << i <<
"]...\n";
1098 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
1099 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
1101 this->_OpCounter += xyc[i];
1108 MS = Teuchos::rcpFromRef(S);
1111 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
1112 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
1117 if ( B == Teuchos::null ) {
1118 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(sc,sc) );
1122 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != sc || B->numCols() != sc, std::invalid_argument,
1123 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of S must be consistent with size of B" );
1128 projectGen(S,X,Y,isBiortho,C,MS,MX,MY);
1130 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(sc,1);
1136 int curssize = sc - rank;
1141 #ifdef ANASAZI_ICGS_DEBUG 1142 *out <<
"Attempting to find orthonormal basis for X...\n";
1144 rank = findBasis(S,MS,*B,
false,curssize);
1146 if (oldrank != -1 && rank != oldrank) {
1152 for (
int i=0; i<sc; i++) {
1153 (*B)(i,oldrank) = oldCoeff(i,0);
1158 if (rank != oldrank) {
1166 for (
int i=0; i<sc; i++) {
1167 oldCoeff(i,0) = (*B)(i,rank);
1174 #ifdef ANASAZI_ICGS_DEBUG 1175 *out <<
"Finished computing basis.\n";
1180 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank,
OrthoError,
1181 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen");
1183 if (rank != oldrank) {
1191 if (numTries <= 0) {
1199 #ifdef ANASAZI_ICGS_DEBUG 1200 *out <<
"Inserting random vector in X[" << rank <<
"]. Attempt " << 10-numTries <<
".\n";
1202 Teuchos::RCP<MV> curS, curMS;
1203 std::vector<int> ind(1);
1205 curS = MVT::CloneViewNonConst(S,ind);
1206 MVT::MvRandom(*curS);
1208 #ifdef ANASAZI_ICGS_DEBUG 1209 *out <<
"Applying operator to random vector.\n";
1211 curMS = MVT::CloneViewNonConst(*MS,ind);
1212 OPT::Apply( *(this->_Op), *curS, *curMS );
1213 this->_OpCounter += MVT::GetNumberVecs(*curS);
1221 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
1222 projectGen(*curS,X,Y,isBiortho,dummyC,curMS,MX,MY);
1228 TEUCHOS_TEST_FOR_EXCEPTION( rank > sc || rank < 0, std::logic_error,
1229 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Debug error in rank variable." );
1231 #ifdef ANASAZI_ICGS_DEBUG 1232 *out <<
"Returning " << rank <<
" from Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
1243 template<
class ScalarType,
class MV,
class OP>
1245 MV &X, Teuchos::RCP<MV> MX,
1246 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
1247 bool completeBasis,
int howMany )
const {
1264 #ifdef ANASAZI_ICGS_DEBUG 1266 Teuchos::RCP<Teuchos::FancyOStream>
1267 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
1268 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
1269 *out <<
"Entering Anasazi::ICGSOrthoManager::findBasis(...)\n";
1272 const ScalarType ONE = SCT::one();
1273 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
1275 int xc = MVT::GetNumberVecs( X );
1277 if (howMany == -1) {
1284 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp ==
true && MX == Teuchos::null, std::logic_error,
1285 "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS.");
1286 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
1287 "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
1292 int xstart = xc - howMany;
1294 for (
int j = xstart; j < xc; j++) {
1303 for (
int i=j+1; i<xc; ++i) {
1308 std::vector<int> index(1);
1310 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1311 Teuchos::RCP<MV> MXj;
1312 if ((this->_hasOp)) {
1314 MXj = MVT::CloneViewNonConst( *MX, index );
1322 std::vector<int> prev_idx( numX );
1323 Teuchos::RCP<const MV> prevX, prevMX;
1326 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
1327 prevX = MVT::CloneView( X, prev_idx );
1329 prevMX = MVT::CloneView( *MX, prev_idx );
1333 bool rankDef =
true;
1338 for (
int numTrials = 0; numTrials < 10; numTrials++) {
1339 #ifdef ANASAZI_ICGS_DEBUG 1340 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
1344 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1345 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
1350 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
1352 #ifdef ANASAZI_ICGS_DEBUG 1353 *out <<
"origNorm = " << origNorm[0] <<
"\n";
1364 #ifdef ANASAZI_ICGS_DEBUG 1365 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
1367 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1374 #ifdef ANASAZI_ICGS_DEBUG 1375 *out <<
"Updating MX[" << j <<
"]...\n";
1377 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1382 MagnitudeType product_norm = product.normOne();
1384 #ifdef ANASAZI_ICGS_DEBUG 1385 *out <<
"newNorm = " << newNorm[0] <<
"\n";
1386 *out <<
"prodoct_norm = " << product_norm <<
"\n";
1390 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
1391 #ifdef ANASAZI_ICGS_DEBUG 1392 if (product_norm/newNorm[0] >= tol_) {
1393 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
1396 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
1401 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1404 #ifdef ANASAZI_ICGS_DEBUG 1405 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
1407 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1408 if ((this->_hasOp)) {
1409 #ifdef ANASAZI_ICGS_DEBUG 1410 *out <<
"Updating MX[" << j <<
"]...\n";
1412 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1416 product_norm = P2.normOne();
1417 #ifdef ANASAZI_ICGS_DEBUG 1418 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
1419 *out <<
"product_norm = " << product_norm <<
"\n";
1421 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
1423 #ifdef ANASAZI_ICGS_DEBUG 1424 if (product_norm/newNorm2[0] >= tol_) {
1425 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
1427 else if (newNorm[0] < newNorm2[0]) {
1428 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
1431 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
1434 MVT::MvInit(*Xj,ZERO);
1435 if ((this->_hasOp)) {
1436 #ifdef ANASAZI_ICGS_DEBUG 1437 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
1439 MVT::MvInit(*MXj,ZERO);
1446 if (numTrials == 0) {
1447 for (
int i=0; i<numX; i++) {
1448 B(i,j) = product(i,0);
1454 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1455 #ifdef ANASAZI_ICGS_DEBUG 1456 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1460 MVT::MvScale( *Xj, ONE/newNorm[0]);
1462 #ifdef ANASAZI_ICGS_DEBUG 1463 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1466 MVT::MvScale( *MXj, ONE/newNorm[0]);
1470 if (numTrials == 0) {
1471 B(j,j) = newNorm[0];
1479 #ifdef ANASAZI_ICGS_DEBUG 1480 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1487 if (completeBasis) {
1489 #ifdef ANASAZI_ICGS_DEBUG 1490 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1492 MVT::MvRandom( *Xj );
1494 #ifdef ANASAZI_ICGS_DEBUG 1495 *out <<
"Updating M*X[" << j <<
"]...\n";
1497 OPT::Apply( *(this->_Op), *Xj, *MXj );
1498 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1509 if (rankDef ==
true) {
1510 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1511 "Anasazi::ICGSOrthoManager::findBasis(): Unable to complete basis" );
1512 #ifdef ANASAZI_ICGS_DEBUG 1513 *out <<
"Returning early, rank " << j <<
" from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1520 #ifdef ANASAZI_ICGS_DEBUG 1521 *out <<
"Returning " << xc <<
" from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1528 #endif // ANASAZI_ICSG_ORTHOMANAGER_HPP int getNumIters() const
Return parameter for re-orthogonalization threshold.
Declaration of basic traits for the multivector type.
~ICGSOrthoManager()
Destructor.
Virtual base class which defines basic traits for the operator type.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthonormErrorMat(const MV &X, Teuchos::RCP< const MV > MX=Teuchos::null) const
This method computes the error in orthonormality of a multivector, measured as the Frobenius norm of ...
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
ICGSOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, int numIters=2, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType eps=0.0, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType tol=0.20)
Constructor specifying the operator defining the inner product as well as the number of orthogonaliza...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
int projectAndNormalizeGen(MV &S, Teuchos::Array< Teuchos::RCP< const MV > > X, Teuchos::Array< Teuchos::RCP< const MV > > Y, bool isBiOrtho, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MS=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MX=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null)), Teuchos::Array< Teuchos::RCP< const MV > > MY=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Applies a series of generic projectors and returns an orthonormal basis for the residual data...
void projectMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Given a list of mutually orthogonal and internally orthonormal bases Q, this method projects a multiv...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
void normMat(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec, Teuchos::RCP< const MV > MX=Teuchos::null) const
Provides the norm induced by the matrix-based inner product.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
int projectAndNormalizeMat(MV &X, Teuchos::Array< Teuchos::RCP< const MV > > Q, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MQ=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Given a set of bases Q[i] and a multivector X, this method computes an orthonormal basis for ...
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
void innerProdMat(const MV &X, const MV &Y, Teuchos::SerialDenseMatrix< int, ScalarType > &Z, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const
Provides a matrix-based inner product.
Exception thrown to signal error in an orthogonalization manager method.
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP< const MV > MX1, Teuchos::RCP< const MV > MX2) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...
void projectGen(MV &S, Teuchos::Array< Teuchos::RCP< const MV > > X, Teuchos::Array< Teuchos::RCP< const MV > > Y, bool isBiOrtho, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< MV > MS=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MX=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null)), Teuchos::Array< Teuchos::RCP< const MV > > MY=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Applies a series of generic projectors.
int normalizeMat(MV &X, Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B=Teuchos::null, Teuchos::RCP< MV > MX=Teuchos::null) const
This method takes a multivector X and attempts to compute an orthonormal basis for ...
void setNumIters(int numIters)
Set parameter for re-orthogonalization threshold.