48 #ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP 49 #define ANASAZI_SVQB_ORTHOMANAGER_HPP 68 template<
class ScalarType,
class MV,
class OP>
296 bool normalize_in )
const;
302 template<
class ScalarType,
class MV,
class OP>
304 :
MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(
" *** "), debug_(debug) {
307 eps_ = lapack.
LAMCH(
'E');
309 std::cout <<
"eps_ == " << eps_ << std::endl;
316 template<
class ScalarType,
class MV,
class OP>
319 const ScalarType ONE = SCT::one();
320 int rank = MVT::GetNumberVecs(X);
323 for (
int i=0; i<rank; i++) {
331 template<
class ScalarType,
class MV,
class OP>
339 int r1 = MVT::GetNumberVecs(X);
340 int r2 = MVT::GetNumberVecs(Y);
349 template<
class ScalarType,
class MV,
class OP>
357 findBasis(X,MX,C,Teuchos::null,Q,
false);
364 template<
class ScalarType,
class MV,
class OP>
371 return findBasis(X,MX,C,B,Q,
true);
378 template<
class ScalarType,
class MV,
class OP>
387 return findBasis(X,MX,C,B,Q,
true);
419 template<
class ScalarType,
class MV,
class OP>
425 bool normalize_in)
const {
427 const ScalarType ONE = SCT::one();
428 const MagnitudeType MONE = SCTM::one();
429 const MagnitudeType ZERO = SCTM::zero();
436 int xc = MVT::GetNumberVecs(X);
437 ptrdiff_t xr = MVT::GetGlobalLength( X );
441 ptrdiff_t qr = (nq == 0) ? 0 : MVT::GetGlobalLength(*Q[0]);
443 std::vector<int> qcs(nq);
444 for (
int i=0; i<nq; i++) {
445 qcs[i] = MVT::GetNumberVecs(*Q[i]);
449 if (normalize_in ==
true && qsize + xc > xr) {
452 "Anasazi::SVQBOrthoManager::findBasis(): Orthogonalization constraints not feasible" );
456 if (normalize_in ==
false && (qsize == 0 || xc == 0)) {
460 else if (normalize_in ==
true && (xc == 0 || xr == 0)) {
463 "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
468 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
477 for (
int i=0; i<nq; i++) {
480 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
482 "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
484 if ( C[i] == Teuchos::null ) {
489 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
492 C[i]->putScalar(ZERO);
502 if (normalize_in ==
true) {
503 if ( B == Teuchos::null ) {
507 "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
510 for (
int i=0; i<xc; i++) {
524 workX = MVT::Clone(X,xc);
527 if (MX == Teuchos::null) {
529 MX = MVT::Clone(X,xc);
530 OPT::Apply(*(this->_Op),X,*MX);
531 this->_OpCounter += MVT::GetNumberVecs(X);
535 MX = Teuchos::rcpFromRef(X);
537 std::vector<ScalarType> normX(xc), invnormX(xc);
544 std::vector<ScalarType> work;
545 std::vector<MagnitudeType> lambda, lambdahi, rwork;
548 int lwork = lapack.
ILAENV(1,
"hetrd",
"VU",xc,-1,-1,-1);
552 "Anasazi::SVQBOrthoManager::findBasis(): Error code from ILAENV" );
554 lwork = (lwork+2)*xc;
557 lwork = (3*xc-2 > 1) ? 3*xc - 2 : 1;
562 workU.reshape(xc,xc);
566 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
567 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
569 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with MX" );
572 bool doGramSchmidt =
true;
574 MagnitudeType tolerance = MONE/SCTM::squareroot(eps_);
577 while (doGramSchmidt) {
587 std::vector<MagnitudeType> normX_mag(xc);
589 for (
int i=0; i<xc; ++i) {
590 normX[i] = normX_mag[i];
591 invnormX[i] = (normX_mag[i] == ZERO) ? ZERO : MONE/normX_mag[i];
595 MVT::MvScale(X,invnormX);
597 MVT::MvScale(*MX,invnormX);
601 std::vector<MagnitudeType> nrm2(xc);
602 std::cout << dbgstr <<
"max post-scale norm: (with/without MX) : ";
603 MagnitudeType maxpsnw = ZERO, maxpsnwo = ZERO;
605 for (
int i=0; i<xc; i++) {
606 maxpsnw = (nrm2[i] > maxpsnw ? nrm2[i] : maxpsnw);
609 for (
int i=0; i<xc; i++) {
610 maxpsnwo = (nrm2[i] > maxpsnwo ? nrm2[i] : maxpsnwo);
612 std::cout <<
"(" << maxpsnw <<
"," << maxpsnwo <<
")" << std::endl;
615 for (
int i=0; i<nq; i++) {
619 for (
int i=0; i<nq; i++) {
620 MVT::MvTimesMatAddMv(-ONE,*Q[i],*newC[i],ONE,X);
623 MVT::MvScale(X,normX);
626 OPT::Apply(*(this->_Op),X,*MX);
627 this->_OpCounter += MVT::GetNumberVecs(X);
635 MagnitudeType maxNorm = ZERO;
636 for (
int j=0; j<xc; j++) {
637 MagnitudeType sum = ZERO;
638 for (
int k=0; k<nq; k++) {
639 for (
int i=0; i<qcs[k]; i++) {
640 sum += SCT::magnitude((*newC[k])(i,j))*SCT::magnitude((*newC[k])(i,j));
643 maxNorm = (sum > maxNorm) ? sum : maxNorm;
647 if (maxNorm < 0.36) {
648 doGramSchmidt =
false;
652 for (
int k=0; k<nq; k++) {
653 for (
int j=0; j<xc; j++) {
654 for (
int i=0; i<qcs[k]; i++) {
655 (*newC[k])(i,j) *= normX[j];
663 for (
int i=0; i<nq; i++) {
665 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
670 for (
int i=0; i<nq; i++) {
677 doGramSchmidt =
false;
685 MagnitudeType condT = tolerance;
687 while (condT >= tolerance) {
695 std::vector<MagnitudeType> Dh(xc), Dhi(xc);
696 for (
int i=0; i<xc; i++) {
697 Dh[i] = SCT::magnitude(SCT::squareroot(XtMX(i,i)));
698 Dhi[i] = (Dh[i] == ZERO ? ZERO : MONE/Dh[i]);
701 for (
int i=0; i<xc; i++) {
702 for (
int j=0; j<xc; j++) {
703 XtMX(i,j) *= Dhi[i]*Dhi[j];
709 lapack.
HEEV(
'V',
'U', xc, XtMX.values(), XtMX.stride(), &lambda[0], &work[0], work.size(), &rwork[0], &info);
710 std::stringstream os;
711 os <<
"Anasazi::SVQBOrthoManager::findBasis(): Error code " << info <<
" from HEEV";
715 std::cout << dbgstr <<
"eigenvalues of XtMX: (";
716 for (
int i=0; i<xc-1; i++) {
717 std::cout << lambda[i] <<
",";
719 std::cout << lambda[xc-1] <<
")" << std::endl;
724 MagnitudeType maxLambda = lambda[xc-1],
725 minLambda = lambda[0];
727 for (
int i=0; i<xc; i++) {
728 if (lambda[i] <= 10*eps_*maxLambda) {
740 lambda[i] = SCTM::squareroot(lambda[i]);
741 lambdahi[i] = MONE/lambda[i];
748 std::vector<int> ind(xc);
749 for (
int i=0; i<xc; i++) {ind[i] = i;}
750 MVT::SetBlock(X,ind,*workX);
754 for (
int j=0; j<xc; j++) {
755 for (
int i=0; i<xc; i++) {
756 workU(i,j) *= Dhi[i]*lambdahi[j];
760 MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,X);
766 if (maxLambda >= tolerance * minLambda) {
768 OPT::Apply(*(this->_Op),X,*MX);
769 this->_OpCounter += MVT::GetNumberVecs(X);
774 MVT::SetBlock(*MX,ind,*workX);
777 MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,*MX);
783 for (
int j=0; j<xc; j++) {
784 for (
int i=0; i<xc; i++) {
785 workU(i,j) = Dh[i] * (*B)(i,j);
789 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
790 for (
int j=0; j<xc ;j++) {
791 for (
int i=0; i<xc; i++) {
792 (*B)(i,j) *= lambda[i];
799 std::cout << dbgstr <<
"augmenting multivec with " << iZeroMax+1 <<
" random directions" << std::endl;
804 std::vector<int> ind2(iZeroMax+1);
805 for (
int i=0; i<iZeroMax+1; i++) {
809 Xnull = MVT::CloneViewNonConst(X,ind2);
810 MVT::MvRandom(*Xnull);
812 MXnull = MVT::CloneViewNonConst(*MX,ind2);
813 OPT::Apply(*(this->_Op),*Xnull,*MXnull);
814 this->_OpCounter += MVT::GetNumberVecs(*Xnull);
815 MXnull = Teuchos::null;
817 Xnull = Teuchos::null;
819 doGramSchmidt =
true;
823 condT = SCTM::magnitude(maxLambda / minLambda);
825 std::cout << dbgstr <<
"condT: " << condT <<
", tolerance = " << tolerance << std::endl;
829 if ((doGramSchmidt ==
false) && (condT > SCTM::squareroot(tolerance))) {
830 doGramSchmidt =
true;
840 std::cout << dbgstr <<
"(numGS,numSVQB,numRand) : " 852 #endif // ANASAZI_SVQB_ORTHOMANAGER_HPP ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector 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...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
OrdinalType ILAENV(const OrdinalType &ispec, const std::string &NAME, const std::string &OPTS, const OrdinalType &N1=-1, const OrdinalType &N2=-1, const OrdinalType &N3=-1, const OrdinalType &N4=-1) const
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
ScalarType LAMCH(const char &CMACH) const
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)
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 ...
void HEEV(const char &JOBZ, const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, MagnitudeType *W, ScalarType *WORK, const OrdinalType &lwork, MagnitudeType *RWORK, OrdinalType *info) 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.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
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.
~SVQBOrthoManager()
Destructor.
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 ...
SVQBOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, bool debug=false)
Constructor specifying re-orthogonalization tolerance.
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 ...
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogErrorMat(const MV &X, const MV &Y, Teuchos::RCP< const MV > MX=Teuchos::null, Teuchos::RCP< const MV > MY=Teuchos::null) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...