47 #ifndef ANASAZI_BASIC_ORTHOMANAGER_HPP
48 #define ANASAZI_BASIC_ORTHOMANAGER_HPP
64 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
65 # include <Teuchos_FancyOStream.hpp>
70 template<
class ScalarType,
class MV,
class OP>
295 MagnitudeType kappa_;
302 bool completeBasis,
int howMany = -1 )
const;
314 template<
class ScalarType,
class MV,
class OP>
319 MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol)
320 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
321 , timerReortho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi::BasicOrthoManager::Re-orthogonalization"))
325 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"eps\" must be non-negative.");
328 eps_ = lapack.
LAMCH(
'E');
333 std::invalid_argument,
334 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"tol\" must be in [0,1].");
341 template<
class ScalarType,
class MV,
class OP>
344 const ScalarType ONE = SCT::one();
345 int rank = MVT::GetNumberVecs(X);
348 for (
int i=0; i<rank; i++) {
358 template<
class ScalarType,
class MV,
class OP>
361 int r1 = MVT::GetNumberVecs(X1);
362 int r2 = MVT::GetNumberVecs(X2);
371 template<
class ScalarType,
class MV,
class OP>
394 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
397 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
399 *out <<
"Entering Anasazi::BasicOrthoManager::projectMat(...)\n";
402 ScalarType ONE = SCT::one();
404 int xc = MVT::GetNumberVecs( X );
405 ptrdiff_t xr = MVT::GetGlobalLength( X );
407 std::vector<int> qcs(nq);
409 if (nq == 0 || xc == 0 || xr == 0) {
410 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
411 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
415 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
424 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
425 *out <<
"Allocating MX...\n";
427 if (MX == Teuchos::null) {
429 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
430 OPT::Apply(*(this->_Op),X,*MX);
431 this->_OpCounter += MVT::GetNumberVecs(X);
436 MX = Teuchos::rcpFromRef(X);
438 int mxc = MVT::GetNumberVecs( *MX );
439 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
443 "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" );
446 "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" );
449 for (
int i=0; i<nq; i++) {
451 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
452 qcs[i] = MVT::GetNumberVecs( *Q[i] );
454 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
456 if ( C[i] == Teuchos::null ) {
461 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
468 std::vector<ScalarType> oldDot( xc );
469 MVT::MvDot( X, *MX, oldDot );
470 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
471 *out <<
"oldDot = { ";
472 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
478 for (
int i=0; i<nq; i++) {
482 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
483 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
485 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
490 if (MQ[i] == Teuchos::null) {
491 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
492 *out <<
"Updating MX via M*X...\n";
494 OPT::Apply( *(this->_Op), X, *MX);
495 this->_OpCounter += MVT::GetNumberVecs(X);
498 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
499 *out <<
"Updating MX via M*Q...\n";
501 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
507 std::vector<ScalarType> newDot(xc);
508 MVT::MvDot( X, *MX, newDot );
509 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
510 *out <<
"newDot = { ";
511 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
516 for (
int j = 0; j < xc; ++j) {
518 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
519 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
520 *out <<
"kappa_*newDot[" <<j<<
"] == " << kappa_*newDot[j] <<
"... another step of Gram-Schmidt.\n";
522 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
525 for (
int i=0; i<nq; i++) {
531 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
532 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
534 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
538 if (MQ[i] == Teuchos::null) {
539 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
540 *out <<
"Updating MX via M*X...\n";
542 OPT::Apply( *(this->_Op), X, *MX);
543 this->_OpCounter += MVT::GetNumberVecs(X);
546 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
547 *out <<
"Updating MX via M*Q...\n";
549 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
557 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
558 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
566 template<
class ScalarType,
class MV,
class OP>
574 int xc = MVT::GetNumberVecs(X);
575 ptrdiff_t xr = MVT::GetGlobalLength(X);
580 if (MX == Teuchos::null) {
582 MX = MVT::Clone(X,xc);
583 OPT::Apply(*(this->_Op),X,*MX);
584 this->_OpCounter += MVT::GetNumberVecs(X);
590 if ( B == Teuchos::null ) {
594 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
595 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
599 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
601 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
603 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
604 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
605 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
607 return findBasis(X, MX, *B,
true );
614 template<
class ScalarType,
class MV,
class OP>
624 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
627 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
629 *out <<
"Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
633 int xc = MVT::GetNumberVecs( X );
634 ptrdiff_t xr = MVT::GetGlobalLength( X );
640 if ( B == Teuchos::null ) {
646 if (MX == Teuchos::null) {
648 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
649 *out <<
"Allocating MX...\n";
651 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
652 OPT::Apply(*(this->_Op),X,*MX);
653 this->_OpCounter += MVT::GetNumberVecs(X);
658 MX = Teuchos::rcpFromRef(X);
661 int mxc = MVT::GetNumberVecs( *MX );
662 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
664 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
666 ptrdiff_t numbas = 0;
667 for (
int i=0; i<nq; i++) {
668 numbas += MVT::GetNumberVecs( *Q[i] );
673 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
676 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
678 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
679 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
681 TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
682 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
685 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
686 *out <<
"Orthogonalizing X against Q...\n";
688 projectMat(X,Q,C,MX,MQ);
696 int curxsize = xc - rank;
701 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
702 *out <<
"Attempting to find orthonormal basis for X...\n";
704 rank = findBasis(X,MX,*B,
false,curxsize);
706 if (oldrank != -1 && rank != oldrank) {
712 for (
int i=0; i<xc; i++) {
713 (*B)(i,oldrank) = oldCoeff(i,0);
718 if (rank != oldrank) {
726 for (
int i=0; i<xc; i++) {
727 oldCoeff(i,0) = (*B)(i,rank);
734 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
735 *out <<
"Finished computing basis.\n";
740 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank,
OrthoError,
741 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
743 if (rank != oldrank) {
759 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
760 *out <<
"Randomizing X[" << rank <<
"]...\n";
763 std::vector<int> ind(1);
765 curX = MVT::CloneViewNonConst(X,ind);
766 MVT::MvRandom(*curX);
768 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
769 *out <<
"Applying operator to random vector.\n";
771 curMX = MVT::CloneViewNonConst(*MX,ind);
772 OPT::Apply( *(this->_Op), *curX, *curMX );
773 this->_OpCounter += MVT::GetNumberVecs(*curX);
782 projectMat(*curX,Q,dummyC,curMX,MQ);
788 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
789 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
791 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
792 *out <<
"Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
803 template<
class ScalarType,
class MV,
class OP>
807 bool completeBasis,
int howMany )
const {
824 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
827 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
829 *out <<
"Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
832 const ScalarType ONE = SCT::one();
833 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
835 int xc = MVT::GetNumberVecs( X );
845 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
847 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
852 int xstart = xc - howMany;
854 for (
int j = xstart; j < xc; j++) {
863 for (
int i=j+1; i<xc; ++i) {
868 std::vector<int> index(1);
872 if ((this->_hasOp)) {
874 MXj = MVT::CloneViewNonConst( *MX, index );
882 std::vector<int> prev_idx( numX );
886 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
887 prevX = MVT::CloneViewNonConst( X, prev_idx );
889 prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
898 for (
int numTrials = 0; numTrials < 10; numTrials++) {
899 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
900 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
905 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
912 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
913 *out <<
"origNorm = " << origNorm[0] <<
"\n";
924 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
925 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
927 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
934 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
935 *out <<
"Updating MX[" << j <<
"]...\n";
937 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
942 MagnitudeType product_norm = product.normOne();
944 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
945 *out <<
"newNorm = " << newNorm[0] <<
"\n";
946 *out <<
"prodoct_norm = " << product_norm <<
"\n";
950 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
951 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
952 if (product_norm/newNorm[0] >= tol_) {
953 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
956 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
964 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
965 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
967 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
968 if ((this->_hasOp)) {
969 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
970 *out <<
"Updating MX[" << j <<
"]...\n";
972 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
976 product_norm = P2.normOne();
977 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
978 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
979 *out <<
"product_norm = " << product_norm <<
"\n";
981 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
983 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
984 if (product_norm/newNorm2[0] >= tol_) {
985 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
987 else if (newNorm[0] < newNorm2[0]) {
988 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
991 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
994 MVT::MvInit(*Xj,ZERO);
995 if ((this->_hasOp)) {
996 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
997 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
999 MVT::MvInit(*MXj,ZERO);
1006 if (numTrials == 0) {
1007 for (
int i=0; i<numX; i++) {
1008 B(i,j) = product(i,0);
1014 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1015 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1016 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1020 MVT::MvScale( *Xj, ONE/newNorm[0]);
1022 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1023 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1026 MVT::MvScale( *MXj, ONE/newNorm[0]);
1030 if (numTrials == 0) {
1031 B(j,j) = newNorm[0];
1039 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1040 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1047 if (completeBasis) {
1049 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1050 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1052 MVT::MvRandom( *Xj );
1054 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1055 *out <<
"Updating M*X[" << j <<
"]...\n";
1057 OPT::Apply( *(this->_Op), *Xj, *MXj );
1058 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1069 if (rankDef ==
true) {
1071 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
1072 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1073 *out <<
"Returning early, rank " << j <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1080 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
1081 *out <<
"Returning " << xc <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1088 #endif // ANASAZI_BASIC_ORTHOMANAGER_HPP
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Teuchos::ScalarTraits< ScalarType >::magnitudeType getKappa() const
Return parameter for re-orthogonalization threshold.
~BasicOrthoManager()
Destructor.
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 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.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
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...
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...
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)
void setKappa(typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa)
Set parameter for re-orthogonalization threshold.
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...
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 ...
Exception thrown to signal error in an orthogonalization manager method.
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 ...
static magnitudeType magnitude(ScalarTypea)
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)...
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 ...
ScalarType LAMCH(const char &CMACH) const
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.