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;
84 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";
Anasazi 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.
Virtual base class which defines basic traits for the operator type.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
void setNumIters(int numIters)
Set parameter for re-orthogonalization threshold.
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...
~ICGSOrthoManager()
Destructor.
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 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 getNumIters() const
Return parameter for re-orthogonalization threshold.
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 ...
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...
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 .
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...
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 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.
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.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
Exception thrown to signal error in an orthogonalization manager method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.