34 #ifndef ANASAZI_BASIC_ORTHOMANAGER_HPP 35 #define ANASAZI_BASIC_ORTHOMANAGER_HPP 51 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 52 # include <Teuchos_FancyOStream.hpp> 57 template<
class ScalarType,
class MV,
class OP>
282 MagnitudeType kappa_;
289 bool completeBasis,
int howMany = -1 )
const;
301 template<
class ScalarType,
class MV,
class OP>
306 MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol)
307 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
308 , timerReortho_(
Teuchos::TimeMonitor::getNewTimer(
"Anasazi::BasicOrthoManager::Re-orthogonalization"))
312 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"eps\" must be non-negative.");
315 eps_ = lapack.
LAMCH(
'E');
320 std::invalid_argument,
321 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"tol\" must be in [0,1].");
328 template<
class ScalarType,
class MV,
class OP>
331 const ScalarType ONE = SCT::one();
332 int rank = MVT::GetNumberVecs(X);
335 for (
int i=0; i<rank; i++) {
345 template<
class ScalarType,
class MV,
class OP>
348 int r1 = MVT::GetNumberVecs(X1);
349 int r2 = MVT::GetNumberVecs(X2);
358 template<
class ScalarType,
class MV,
class OP>
381 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 384 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
386 *out <<
"Entering Anasazi::BasicOrthoManager::projectMat(...)\n";
389 ScalarType ONE = SCT::one();
391 int xc = MVT::GetNumberVecs( X );
392 ptrdiff_t xr = MVT::GetGlobalLength( X );
394 std::vector<int> qcs(nq);
396 if (nq == 0 || xc == 0 || xr == 0) {
397 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 398 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
402 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
411 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 412 *out <<
"Allocating MX...\n";
414 if (MX == Teuchos::null) {
416 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
417 OPT::Apply(*(this->_Op),X,*MX);
418 this->_OpCounter += MVT::GetNumberVecs(X);
423 MX = Teuchos::rcpFromRef(X);
425 int mxc = MVT::GetNumberVecs( *MX );
426 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
430 "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" );
433 "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" );
437 for (
int i=0; i<nq; i++) {
439 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
440 qcs[i] = MVT::GetNumberVecs( *Q[i] );
442 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
446 if ( C[i] == Teuchos::null ) {
451 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
458 std::vector<ScalarType> oldDot( xc );
459 MVT::MvDot( X, *MX, oldDot );
460 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 461 *out <<
"oldDot = { ";
462 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
468 for (
int i=0; i<nq; i++) {
472 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 473 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
475 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
480 if (MQ[i] == Teuchos::null) {
481 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 482 *out <<
"Updating MX via M*X...\n";
484 OPT::Apply( *(this->_Op), X, *MX);
485 this->_OpCounter += MVT::GetNumberVecs(X);
488 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 489 *out <<
"Updating MX via M*Q...\n";
491 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
497 std::vector<ScalarType> newDot(xc);
498 MVT::MvDot( X, *MX, newDot );
499 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 500 *out <<
"newDot = { ";
501 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
506 for (
int j = 0; j < xc; ++j) {
508 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
509 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 510 *out <<
"kappa_*newDot[" <<j<<
"] == " << kappa_*newDot[j] <<
"... another step of Gram-Schmidt.\n";
512 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 515 for (
int i=0; i<nq; i++) {
521 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 522 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
524 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
528 if (MQ[i] == Teuchos::null) {
529 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 530 *out <<
"Updating MX via M*X...\n";
532 OPT::Apply( *(this->_Op), X, *MX);
533 this->_OpCounter += MVT::GetNumberVecs(X);
536 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 537 *out <<
"Updating MX via M*Q...\n";
539 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
547 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 548 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
556 template<
class ScalarType,
class MV,
class OP>
564 int xc = MVT::GetNumberVecs(X);
565 ptrdiff_t xr = MVT::GetGlobalLength(X);
570 if (MX == Teuchos::null) {
572 MX = MVT::Clone(X,xc);
573 OPT::Apply(*(this->_Op),X,*MX);
574 this->_OpCounter += MVT::GetNumberVecs(X);
580 if ( B == Teuchos::null ) {
584 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
585 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
589 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
591 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
593 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
595 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
597 return findBasis(X, MX, *B,
true );
604 template<
class ScalarType,
class MV,
class OP>
614 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 617 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
619 *out <<
"Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
623 int xc = MVT::GetNumberVecs( X );
624 ptrdiff_t xr = MVT::GetGlobalLength( X );
630 if ( B == Teuchos::null ) {
636 if (MX == Teuchos::null) {
638 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 639 *out <<
"Allocating MX...\n";
641 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
642 OPT::Apply(*(this->_Op),X,*MX);
643 this->_OpCounter += MVT::GetNumberVecs(X);
648 MX = Teuchos::rcpFromRef(X);
651 int mxc = MVT::GetNumberVecs( *MX );
652 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
654 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
656 ptrdiff_t numbas = 0;
657 for (
int i=0; i<nq; i++) {
658 numbas += MVT::GetNumberVecs( *Q[i] );
663 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
666 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
669 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
672 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
675 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 676 *out <<
"Orthogonalizing X against Q...\n";
678 projectMat(X,Q,C,MX,MQ);
686 int curxsize = xc - rank;
691 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 692 *out <<
"Attempting to find orthonormal basis for X...\n";
694 rank = findBasis(X,MX,*B,
false,curxsize);
696 if (oldrank != -1 && rank != oldrank) {
702 for (
int i=0; i<xc; i++) {
703 (*B)(i,oldrank) = oldCoeff(i,0);
708 if (rank != oldrank) {
716 for (
int i=0; i<xc; i++) {
717 oldCoeff(i,0) = (*B)(i,rank);
724 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 725 *out <<
"Finished computing basis.\n";
731 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
733 if (rank != oldrank) {
749 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 750 *out <<
"Randomizing X[" << rank <<
"]...\n";
753 std::vector<int> ind(1);
755 curX = MVT::CloneViewNonConst(X,ind);
756 MVT::MvRandom(*curX);
758 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 759 *out <<
"Applying operator to random vector.\n";
761 curMX = MVT::CloneViewNonConst(*MX,ind);
762 OPT::Apply( *(this->_Op), *curX, *curMX );
763 this->_OpCounter += MVT::GetNumberVecs(*curX);
772 projectMat(*curX,Q,dummyC,curMX,MQ);
779 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
781 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 782 *out <<
"Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
793 template<
class ScalarType,
class MV,
class OP>
797 bool completeBasis,
int howMany )
const {
814 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 817 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
819 *out <<
"Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
822 const ScalarType ONE = SCT::one();
823 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
825 int xc = MVT::GetNumberVecs( X );
835 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
837 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
842 int xstart = xc - howMany;
844 for (
int j = xstart; j < xc; j++) {
853 for (
int i=j+1; i<xc; ++i) {
858 std::vector<int> index(1);
862 if ((this->_hasOp)) {
864 MXj = MVT::CloneViewNonConst( *MX, index );
872 std::vector<int> prev_idx( numX );
876 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
877 prevX = MVT::CloneViewNonConst( X, prev_idx );
879 prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
888 for (
int numTrials = 0; numTrials < 10; numTrials++) {
889 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 890 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
895 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
902 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 903 *out <<
"origNorm = " << origNorm[0] <<
"\n";
914 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 915 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
917 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
924 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 925 *out <<
"Updating MX[" << j <<
"]...\n";
927 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
932 MagnitudeType product_norm = product.normOne();
934 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 935 *out <<
"newNorm = " << newNorm[0] <<
"\n";
936 *out <<
"prodoct_norm = " << product_norm <<
"\n";
940 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
941 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 942 if (product_norm/newNorm[0] >= tol_) {
943 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
946 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
954 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 955 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
957 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
958 if ((this->_hasOp)) {
959 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 960 *out <<
"Updating MX[" << j <<
"]...\n";
962 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
966 product_norm = P2.normOne();
967 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 968 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
969 *out <<
"product_norm = " << product_norm <<
"\n";
971 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
973 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 974 if (product_norm/newNorm2[0] >= tol_) {
975 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
977 else if (newNorm[0] < newNorm2[0]) {
978 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
981 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
984 MVT::MvInit(*Xj,ZERO);
985 if ((this->_hasOp)) {
986 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 987 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
989 MVT::MvInit(*MXj,ZERO);
996 if (numTrials == 0) {
997 for (
int i=0; i<numX; i++) {
998 B(i,j) = product(i,0);
1004 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1005 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1006 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1010 MVT::MvScale( *Xj, ONE/newNorm[0]);
1012 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1013 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1016 MVT::MvScale( *MXj, ONE/newNorm[0]);
1020 if (numTrials == 0) {
1021 B(j,j) = newNorm[0];
1029 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1030 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1037 if (completeBasis) {
1039 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1040 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1042 MVT::MvRandom( *Xj );
1044 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1045 *out <<
"Updating M*X[" << j <<
"]...\n";
1047 OPT::Apply( *(this->_Op), *Xj, *MXj );
1048 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1059 if (rankDef ==
true) {
1061 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
1062 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1063 *out <<
"Returning early, rank " << j <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1070 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1071 *out <<
"Returning " << xc <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1078 #endif // ANASAZI_BASIC_ORTHOMANAGER_HPP ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Teuchos::ScalarTraits< ScalarType >::magnitudeType getKappa() const
Return parameter for re-orthogonalization threshold.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
~BasicOrthoManager()
Destructor.
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 ...
Declaration of basic traits for the multivector type.
basic_FancyOStream & setShowProcRank(const bool showProcRank)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
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 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 ...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
ScalarType LAMCH(const char CMACH) const
void setKappa(typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa)
Set parameter for re-orthogonalization threshold.
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.
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
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.
static magnitudeType magnitude(ScalarType a)
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 ...
BasicOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa=1.41421356, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType eps=0.0, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType tol=0.20)
Constructor specifying re-orthogonalization tolerance.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
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...