34 #ifndef ANASAZI_ICSG_ORTHOMANAGER_HPP 35 #define ANASAZI_ICSG_ORTHOMANAGER_HPP 48 #include "Teuchos_TimeMonitor.hpp" 49 #include "Teuchos_LAPACK.hpp" 50 #include "Teuchos_BLAS.hpp" 51 #ifdef ANASAZI_ICGS_DEBUG 52 # include <Teuchos_FancyOStream.hpp> 57 template<
class ScalarType,
class MV,
class OP>
61 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
62 typedef Teuchos::ScalarTraits<ScalarType> SCT;
70 ICGSOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null,
int numIters = 2,
72 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
73 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
167 Teuchos::Array<Teuchos::RCP<const MV> > X,
168 Teuchos::Array<Teuchos::RCP<const MV> > Y,
170 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
171 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
172 Teuchos::RCP<MV> MS = Teuchos::null,
173 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
174 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
271 Teuchos::Array<Teuchos::RCP<const MV> > X,
272 Teuchos::Array<Teuchos::RCP<const MV> > Y,
274 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
275 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
276 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
277 Teuchos::RCP<MV> MS = Teuchos::null,
278 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
279 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
303 Teuchos::Array<Teuchos::RCP<const MV> > Q,
304 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
305 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
306 Teuchos::RCP<MV> MX = Teuchos::null,
307 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
321 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
322 Teuchos::RCP<MV> MX = Teuchos::null
336 Teuchos::Array<Teuchos::RCP<const MV> > Q,
337 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
338 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
339 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
340 Teuchos::RCP<MV> MX = Teuchos::null,
341 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
353 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
360 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
361 orthogErrorMat(
const MV &X1,
const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2)
const;
370 numIters_ = numIters;
371 TEUCHOS_TEST_FOR_EXCEPTION(numIters_ < 1,std::invalid_argument,
372 "Anasazi::ICGSOrthoManager::setNumIters(): input must be >= 1.");
388 int findBasis(MV &X, Teuchos::RCP<MV> MX,
389 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
390 bool completeBasis,
int howMany = -1)
const;
397 template<
class ScalarType,
class MV,
class OP>
400 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
401 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol) :
405 TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
406 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative.");
408 Teuchos::LAPACK<int,MagnitudeType> lapack;
409 eps_ = lapack.LAMCH(
'E');
410 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.50);
412 TEUCHOS_TEST_FOR_EXCEPTION(
413 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
414 std::invalid_argument,
415 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1].");
422 template<
class ScalarType,
class MV,
class OP>
423 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
425 const ScalarType ONE = SCT::one();
426 int rank = MVT::GetNumberVecs(X);
427 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
429 for (
int i=0; i<rank; i++) {
432 return xTx.normFrobenius();
439 template<
class ScalarType,
class MV,
class OP>
440 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
442 int r1 = MVT::GetNumberVecs(X1);
443 int r2 = MVT::GetNumberVecs(X2);
444 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
446 return xTx.normFrobenius();
452 template<
class ScalarType,
class MV,
class OP>
455 Teuchos::Array<Teuchos::RCP<const MV> > Q,
456 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
458 Teuchos::Array<Teuchos::RCP<const MV> > MQ
461 projectGen(X,Q,Q,
true,C,MX,MQ,MQ);
468 template<
class ScalarType,
class MV,
class OP>
471 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
472 Teuchos::RCP<MV> MX)
const 477 int xc = MVT::GetNumberVecs(X);
478 ptrdiff_t xr = MVT::GetGlobalLength(X);
483 if (MX == Teuchos::null) {
485 MX = MVT::Clone(X,xc);
486 OPT::Apply(*(this->_Op),X,*MX);
487 this->_OpCounter += MVT::GetNumberVecs(X);
493 if ( B == Teuchos::null ) {
494 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
497 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
498 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
501 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
502 "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" );
503 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
504 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
505 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
506 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
507 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
508 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
510 return findBasis(X, MX, *B,
true );
517 template<
class ScalarType,
class MV,
class OP>
520 Teuchos::Array<Teuchos::RCP<const MV> > Q,
521 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
522 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
524 Teuchos::Array<Teuchos::RCP<const MV> > MQ
527 return projectAndNormalizeGen(X,Q,Q,
true,C,B,MX,MQ,MQ);
533 template<
class ScalarType,
class MV,
class OP>
536 Teuchos::Array<Teuchos::RCP<const MV> > X,
537 Teuchos::Array<Teuchos::RCP<const MV> > Y,
539 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
541 Teuchos::Array<Teuchos::RCP<const MV> > MX,
542 Teuchos::Array<Teuchos::RCP<const MV> > MY
560 #ifdef ANASAZI_ICGS_DEBUG 562 Teuchos::RCP<Teuchos::FancyOStream>
563 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
564 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
565 *out <<
"Entering Anasazi::ICGSOrthoManager::projectGen(...)\n";
568 const ScalarType ONE = SCT::one();
569 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
570 Teuchos::LAPACK<int,ScalarType> lapack;
571 Teuchos::BLAS<int,ScalarType> blas;
573 int sc = MVT::GetNumberVecs( S );
574 ptrdiff_t sr = MVT::GetGlobalLength( S );
575 int numxy = X.length();
576 TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument,
577 "Anasazi::ICGSOrthoManager::projectGen(): X and Y must contain the same number of multivectors.");
578 std::vector<int> xyc(numxy);
580 if (numxy == 0 || sc == 0 || sr == 0) {
581 #ifdef ANASAZI_ICGS_DEBUG 582 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
596 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument,
597 "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." );
600 if (this->_hasOp ==
true) {
601 if (MS != Teuchos::null) {
602 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MS) != sr, std::invalid_argument,
603 "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." );
604 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
605 "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." );
610 ptrdiff_t sumxyc = 0;
613 for (
int i=0; i<numxy; i++) {
614 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
615 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*X[i]) != sr, std::invalid_argument,
616 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"] length not consistent with S." );
617 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*Y[i]) != sr, std::invalid_argument,
618 "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i <<
"] length not consistent with S." );
619 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
620 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"] and Y[" << i <<
"] widths not consistent." );
622 xyc[i] = MVT::GetNumberVecs( *X[i] );
623 TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument,
624 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"],Y[" << i <<
"] have less rows than columns, and therefore cannot be full rank." );
628 if ( C[i] == Teuchos::null ) {
629 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
632 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument,
633 "Anasazi::ICGSOrthoManager::projectGen(): Size of Q not consistent with size of C." );
637 if (MX[i] != Teuchos::null) {
638 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
639 "Anasazi::ICGSOrthoManager::projectGen(): Size of MX[" << i <<
"] not consistent with size of X[" << i <<
"]." );
644 if (MY[i] != Teuchos::null) {
645 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
646 "Anasazi::ICGSOrthoManager::projectGen(): Size of MY[" << i <<
"] not consistent with size of Y[" << i <<
"]." );
654 TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
655 "Anasazi::ICGSOrthoManager::projectGen(): " 656 << (X[i] == Teuchos::null ?
"Y[" :
"X[") << i <<
"] was provided but " 657 << (X[i] == Teuchos::null ?
"X[" :
"Y[") << i <<
"] was not.");
662 TEUCHOS_TEST_FOR_EXCEPTION(sumxyc > sr, std::invalid_argument,
663 "Anasazi::ICGSOrthoManager::projectGen(): dimension of all X[i],Y[i] is " 664 << sumxyc <<
", but length of vectors is only " << sr <<
". This is infeasible.");
699 if (MXmissing == 0) {
702 if (MYmissing == 0) {
709 switch (whichAlloc) {
725 if (MS == Teuchos::null) {
726 #ifdef ANASAZI_ICGS_DEBUG 727 *out <<
"Allocating MS...\n";
729 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
730 OPT::Apply(*(this->_Op),S,*MS);
731 this->_OpCounter += MVT::GetNumberVecs(S);
735 if (whichAlloc == 0) {
737 for (
int i=0; i<numxy; i++) {
738 if (MX[i] == Teuchos::null) {
739 #ifdef ANASAZI_ICGS_DEBUG 740 *out <<
"Allocating MX[" << i <<
"]...\n";
742 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
743 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
745 this->_OpCounter += xyc[i];
752 MS = Teuchos::rcpFromRef(S);
755 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
756 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
765 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > YMX(numxy);
766 if (isBiortho ==
false) {
767 for (
int i=0; i<numxy; i++) {
768 #ifdef ANASAZI_ICGS_DEBUG 769 *out <<
"Computing YMX[" << i <<
"] and its Cholesky factorization...\n";
771 YMX[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],xyc[i]) );
773 #ifdef ANASAZI_ICGS_DEBUG 780 MagnitudeType err = ZERO;
781 for (
int jj=0; jj<YMX[i]->numCols(); jj++) {
782 err =+ SCT::magnitude(SCT::imag((*YMX[i])(jj,jj)));
783 for (
int ii=jj; ii<YMX[i]->numRows(); ii++) {
784 err += SCT::magnitude( (*YMX[i])(ii,jj) - SCT::conjugate((*YMX[i])(jj,ii)) );
787 *out <<
"Symmetry error in YMX[" << i <<
"] == " << err <<
"\n";
792 lapack.POTRF(
'U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info);
793 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
794 "Anasazi::ICGSOrthoManager::projectGen(): Error computing Cholesky factorization of Y[i]^T * M * X[i] using POTRF: returned info " << info);
799 #ifdef ANASAZI_ICGS_DEBUG 800 std::vector<MagnitudeType> oldNorms(sc);
802 *out <<
"oldNorms = { ";
803 std::copy(oldNorms.begin(), oldNorms.end(), std::ostream_iterator<MagnitudeType>(*out,
" "));
809 Teuchos::Array<Teuchos::SerialDenseMatrix<int,ScalarType> > Ccur(numxy);
810 for (
int i=0; i<numxy; i++) {
811 C[i]->putScalar(ZERO);
812 Ccur[i].reshape(C[i]->numRows(),C[i]->numCols());
815 for (
int iter=0; iter < numIters_; iter++) {
816 #ifdef ANASAZI_ICGS_DEBUG 817 *out <<
"beginning iteration " << iter+1 <<
"\n";
821 for (
int i=0; i<numxy; i++) {
824 if (isBiortho ==
false) {
827 lapack.POTRS(
'U',YMX[i]->numCols(),Ccur[i].numCols(),
828 YMX[i]->values(),YMX[i]->stride(),
829 Ccur[i].values(),Ccur[i].stride(), &info);
830 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
831 "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info <<
" from lapack::POTRS." );
835 #ifdef ANASAZI_ICGS_DEBUG 836 *out <<
"Applying projector P_{X[" << i <<
"],Y[" << i <<
"]}...\n";
838 MVT::MvTimesMatAddMv( -ONE, *X[i], Ccur[i], ONE, S );
845 #ifdef ANASAZI_ICGS_DEBUG 846 *out <<
"Updating MS...\n";
848 OPT::Apply( *(this->_Op), S, *MS);
849 this->_OpCounter += sc;
851 else if (updateMS == 2) {
852 #ifdef ANASAZI_ICGS_DEBUG 853 *out <<
"Updating MS...\n";
855 MVT::MvTimesMatAddMv( -ONE, *MX[i], Ccur[i], ONE, *MS );
860 #ifdef ANASAZI_ICGS_DEBUG 861 std::vector<MagnitudeType> newNorms(sc);
863 *out <<
"newNorms = { ";
864 std::copy(newNorms.begin(), newNorms.end(), std::ostream_iterator<MagnitudeType>(*out,
" "));
869 #ifdef ANASAZI_ICGS_DEBUG 870 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
878 template<
class ScalarType,
class MV,
class OP>
881 Teuchos::Array<Teuchos::RCP<const MV> > X,
882 Teuchos::Array<Teuchos::RCP<const MV> > Y,
884 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
885 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
887 Teuchos::Array<Teuchos::RCP<const MV> > MX,
888 Teuchos::Array<Teuchos::RCP<const MV> > MY
906 #ifdef ANASAZI_ICGS_DEBUG 908 Teuchos::RCP<Teuchos::FancyOStream>
909 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
910 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
911 *out <<
"Entering Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
915 int sc = MVT::GetNumberVecs( S );
916 ptrdiff_t sr = MVT::GetGlobalLength( S );
917 int numxy = X.length();
918 TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument,
919 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X and Y must contain the same number of multivectors.");
920 std::vector<int> xyc(numxy);
922 if (sc == 0 || sr == 0) {
923 #ifdef ANASAZI_ICGS_DEBUG 924 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
938 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument,
939 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." );
942 if (this->_hasOp ==
true) {
943 if (MS != Teuchos::null) {
944 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MS) != sr, std::invalid_argument,
945 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." );
946 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
947 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." );
952 ptrdiff_t sumxyc = 0;
955 for (
int i=0; i<numxy; i++) {
956 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
957 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*X[i]) != sr, std::invalid_argument,
958 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"] length not consistent with S." );
959 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*Y[i]) != sr, std::invalid_argument,
960 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i <<
"] length not consistent with S." );
961 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
962 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"] and Y[" << i <<
"] widths not consistent." );
964 xyc[i] = MVT::GetNumberVecs( *X[i] );
965 TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument,
966 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"],Y[" << i <<
"] have less rows than columns, and therefore cannot be full rank." );
970 if ( C[i] == Teuchos::null ) {
971 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
974 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument,
975 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of Q not consistent with size of C." );
979 if (MX[i] != Teuchos::null) {
980 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
981 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MX[" << i <<
"] not consistent with size of X[" << i <<
"]." );
986 if (MY[i] != Teuchos::null) {
987 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
988 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MY[" << i <<
"] not consistent with size of Y[" << i <<
"]." );
996 TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
997 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): " 998 << (X[i] == Teuchos::null ?
"Y[" :
"X[") << i <<
"] was provided but " 999 << (X[i] == Teuchos::null ?
"X[" :
"Y[") << i <<
"] was not.");
1004 TEUCHOS_TEST_FOR_EXCEPTION(sumxyc + sc > sr, std::invalid_argument,
1005 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): dimension of all X[i],Y[i] is " 1006 << sumxyc <<
" and requested " << sc <<
"-dimensional basis, but length of vectors is only " 1007 << sr <<
". This is infeasible.");
1042 if (MXmissing == 0) {
1045 if (MYmissing == 0) {
1052 switch (whichAlloc) {
1068 if (MS == Teuchos::null) {
1069 #ifdef ANASAZI_ICGS_DEBUG 1070 *out <<
"Allocating MS...\n";
1072 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
1073 OPT::Apply(*(this->_Op),S,*MS);
1074 this->_OpCounter += MVT::GetNumberVecs(S);
1078 if (whichAlloc == 0) {
1080 for (
int i=0; i<numxy; i++) {
1081 if (MX[i] == Teuchos::null) {
1082 #ifdef ANASAZI_ICGS_DEBUG 1083 *out <<
"Allocating MX[" << i <<
"]...\n";
1085 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
1086 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
1088 this->_OpCounter += xyc[i];
1095 MS = Teuchos::rcpFromRef(S);
1098 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
1099 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
1104 if ( B == Teuchos::null ) {
1105 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(sc,sc) );
1109 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != sc || B->numCols() != sc, std::invalid_argument,
1110 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of S must be consistent with size of B" );
1115 projectGen(S,X,Y,isBiortho,C,MS,MX,MY);
1117 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(sc,1);
1123 int curssize = sc - rank;
1128 #ifdef ANASAZI_ICGS_DEBUG 1129 *out <<
"Attempting to find orthonormal basis for X...\n";
1131 rank = findBasis(S,MS,*B,
false,curssize);
1133 if (oldrank != -1 && rank != oldrank) {
1139 for (
int i=0; i<sc; i++) {
1140 (*B)(i,oldrank) = oldCoeff(i,0);
1145 if (rank != oldrank) {
1153 for (
int i=0; i<sc; i++) {
1154 oldCoeff(i,0) = (*B)(i,rank);
1161 #ifdef ANASAZI_ICGS_DEBUG 1162 *out <<
"Finished computing basis.\n";
1167 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank,
OrthoError,
1168 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen");
1170 if (rank != oldrank) {
1178 if (numTries <= 0) {
1186 #ifdef ANASAZI_ICGS_DEBUG 1187 *out <<
"Inserting random vector in X[" << rank <<
"]. Attempt " << 10-numTries <<
".\n";
1189 Teuchos::RCP<MV> curS, curMS;
1190 std::vector<int> ind(1);
1192 curS = MVT::CloneViewNonConst(S,ind);
1193 MVT::MvRandom(*curS);
1195 #ifdef ANASAZI_ICGS_DEBUG 1196 *out <<
"Applying operator to random vector.\n";
1198 curMS = MVT::CloneViewNonConst(*MS,ind);
1199 OPT::Apply( *(this->_Op), *curS, *curMS );
1200 this->_OpCounter += MVT::GetNumberVecs(*curS);
1208 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
1209 projectGen(*curS,X,Y,isBiortho,dummyC,curMS,MX,MY);
1215 TEUCHOS_TEST_FOR_EXCEPTION( rank > sc || rank < 0, std::logic_error,
1216 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Debug error in rank variable." );
1218 #ifdef ANASAZI_ICGS_DEBUG 1219 *out <<
"Returning " << rank <<
" from Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
1230 template<
class ScalarType,
class MV,
class OP>
1232 MV &X, Teuchos::RCP<MV> MX,
1233 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
1234 bool completeBasis,
int howMany )
const {
1251 #ifdef ANASAZI_ICGS_DEBUG 1253 Teuchos::RCP<Teuchos::FancyOStream>
1254 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
1255 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
1256 *out <<
"Entering Anasazi::ICGSOrthoManager::findBasis(...)\n";
1259 const ScalarType ONE = SCT::one();
1260 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
1262 int xc = MVT::GetNumberVecs( X );
1264 if (howMany == -1) {
1271 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp ==
true && MX == Teuchos::null, std::logic_error,
1272 "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS.");
1273 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
1274 "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
1279 int xstart = xc - howMany;
1281 for (
int j = xstart; j < xc; j++) {
1290 for (
int i=j+1; i<xc; ++i) {
1295 std::vector<int> index(1);
1297 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
1298 Teuchos::RCP<MV> MXj;
1299 if ((this->_hasOp)) {
1301 MXj = MVT::CloneViewNonConst( *MX, index );
1309 std::vector<int> prev_idx( numX );
1310 Teuchos::RCP<const MV> prevX, prevMX;
1313 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
1314 prevX = MVT::CloneView( X, prev_idx );
1316 prevMX = MVT::CloneView( *MX, prev_idx );
1320 bool rankDef =
true;
1325 for (
int numTrials = 0; numTrials < 10; numTrials++) {
1326 #ifdef ANASAZI_ICGS_DEBUG 1327 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
1331 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
1332 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
1337 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
1339 #ifdef ANASAZI_ICGS_DEBUG 1340 *out <<
"origNorm = " << origNorm[0] <<
"\n";
1351 #ifdef ANASAZI_ICGS_DEBUG 1352 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
1354 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1361 #ifdef ANASAZI_ICGS_DEBUG 1362 *out <<
"Updating MX[" << j <<
"]...\n";
1364 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1369 MagnitudeType product_norm = product.normOne();
1371 #ifdef ANASAZI_ICGS_DEBUG 1372 *out <<
"newNorm = " << newNorm[0] <<
"\n";
1373 *out <<
"prodoct_norm = " << product_norm <<
"\n";
1377 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
1378 #ifdef ANASAZI_ICGS_DEBUG 1379 if (product_norm/newNorm[0] >= tol_) {
1380 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
1383 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
1388 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
1391 #ifdef ANASAZI_ICGS_DEBUG 1392 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
1394 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1395 if ((this->_hasOp)) {
1396 #ifdef ANASAZI_ICGS_DEBUG 1397 *out <<
"Updating MX[" << j <<
"]...\n";
1399 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1403 product_norm = P2.normOne();
1404 #ifdef ANASAZI_ICGS_DEBUG 1405 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
1406 *out <<
"product_norm = " << product_norm <<
"\n";
1408 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
1410 #ifdef ANASAZI_ICGS_DEBUG 1411 if (product_norm/newNorm2[0] >= tol_) {
1412 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
1414 else if (newNorm[0] < newNorm2[0]) {
1415 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
1418 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
1421 MVT::MvInit(*Xj,ZERO);
1422 if ((this->_hasOp)) {
1423 #ifdef ANASAZI_ICGS_DEBUG 1424 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
1426 MVT::MvInit(*MXj,ZERO);
1433 if (numTrials == 0) {
1434 for (
int i=0; i<numX; i++) {
1435 B(i,j) = product(i,0);
1441 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1442 #ifdef ANASAZI_ICGS_DEBUG 1443 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1447 MVT::MvScale( *Xj, ONE/newNorm[0]);
1449 #ifdef ANASAZI_ICGS_DEBUG 1450 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1453 MVT::MvScale( *MXj, ONE/newNorm[0]);
1457 if (numTrials == 0) {
1458 B(j,j) = newNorm[0];
1466 #ifdef ANASAZI_ICGS_DEBUG 1467 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1474 if (completeBasis) {
1476 #ifdef ANASAZI_ICGS_DEBUG 1477 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1479 MVT::MvRandom( *Xj );
1481 #ifdef ANASAZI_ICGS_DEBUG 1482 *out <<
"Updating M*X[" << j <<
"]...\n";
1484 OPT::Apply( *(this->_Op), *Xj, *MXj );
1485 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1496 if (rankDef ==
true) {
1497 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1498 "Anasazi::ICGSOrthoManager::findBasis(): Unable to complete basis" );
1499 #ifdef ANASAZI_ICGS_DEBUG 1500 *out <<
"Returning early, rank " << j <<
" from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1507 #ifdef ANASAZI_ICGS_DEBUG 1508 *out <<
"Returning " << xc <<
" from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1515 #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.