35 #ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP 36 #define ANASAZI_SVQB_ORTHOMANAGER_HPP 55 template<
class ScalarType,
class MV,
class OP>
283 bool normalize_in )
const;
289 template<
class ScalarType,
class MV,
class OP>
291 :
MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(
" *** "), debug_(debug) {
294 eps_ = lapack.
LAMCH(
'E');
296 std::cout <<
"eps_ == " << eps_ << std::endl;
303 template<
class ScalarType,
class MV,
class OP>
306 const ScalarType ONE = SCT::one();
307 int rank = MVT::GetNumberVecs(X);
310 for (
int i=0; i<rank; i++) {
318 template<
class ScalarType,
class MV,
class OP>
326 int r1 = MVT::GetNumberVecs(X);
327 int r2 = MVT::GetNumberVecs(Y);
336 template<
class ScalarType,
class MV,
class OP>
344 findBasis(X,MX,C,Teuchos::null,Q,
false);
351 template<
class ScalarType,
class MV,
class OP>
358 return findBasis(X,MX,C,B,Q,
true);
365 template<
class ScalarType,
class MV,
class OP>
374 return findBasis(X,MX,C,B,Q,
true);
406 template<
class ScalarType,
class MV,
class OP>
412 bool normalize_in)
const {
414 const ScalarType ONE = SCT::one();
415 const MagnitudeType MONE = SCTM::one();
416 const MagnitudeType ZERO = SCTM::zero();
423 int xc = MVT::GetNumberVecs(X);
424 ptrdiff_t xr = MVT::GetGlobalLength( X );
428 ptrdiff_t qr = (nq == 0) ? 0 : MVT::GetGlobalLength(*Q[0]);
430 std::vector<int> qcs(nq);
431 for (
int i=0; i<nq; i++) {
432 qcs[i] = MVT::GetNumberVecs(*Q[i]);
436 if (normalize_in ==
true && qsize + xc > xr) {
439 "Anasazi::SVQBOrthoManager::findBasis(): Orthogonalization constraints not feasible" );
443 if (normalize_in ==
false && (qsize == 0 || xc == 0)) {
447 else if (normalize_in ==
true && (xc == 0 || xr == 0)) {
450 "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
455 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
464 for (
int i=0; i<nq; i++) {
467 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
469 "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
471 if ( C[i] == Teuchos::null ) {
476 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
479 C[i]->putScalar(ZERO);
489 if (normalize_in ==
true) {
490 if ( B == Teuchos::null ) {
494 "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
497 for (
int i=0; i<xc; i++) {
511 workX = MVT::Clone(X,xc);
514 if (MX == Teuchos::null) {
516 MX = MVT::Clone(X,xc);
517 OPT::Apply(*(this->_Op),X,*MX);
518 this->_OpCounter += MVT::GetNumberVecs(X);
522 MX = Teuchos::rcpFromRef(X);
524 std::vector<ScalarType> normX(xc), invnormX(xc);
531 std::vector<ScalarType> work;
532 std::vector<MagnitudeType> lambda, lambdahi, rwork;
535 int lwork = lapack.
ILAENV(1,
"hetrd",
"VU",xc,-1,-1,-1);
539 "Anasazi::SVQBOrthoManager::findBasis(): Error code from ILAENV" );
541 lwork = (lwork+2)*xc;
544 lwork = (3*xc-2 > 1) ? 3*xc - 2 : 1;
549 workU.reshape(xc,xc);
553 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
554 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
556 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with MX" );
559 bool doGramSchmidt =
true;
561 MagnitudeType tolerance = MONE/SCTM::squareroot(eps_);
564 while (doGramSchmidt) {
574 std::vector<MagnitudeType> normX_mag(xc);
576 for (
int i=0; i<xc; ++i) {
577 normX[i] = normX_mag[i];
578 invnormX[i] = (normX_mag[i] == ZERO) ? ZERO : MONE/normX_mag[i];
582 MVT::MvScale(X,invnormX);
584 MVT::MvScale(*MX,invnormX);
588 std::vector<MagnitudeType> nrm2(xc);
589 std::cout << dbgstr <<
"max post-scale norm: (with/without MX) : ";
590 MagnitudeType maxpsnw = ZERO, maxpsnwo = ZERO;
592 for (
int i=0; i<xc; i++) {
593 maxpsnw = (nrm2[i] > maxpsnw ? nrm2[i] : maxpsnw);
596 for (
int i=0; i<xc; i++) {
597 maxpsnwo = (nrm2[i] > maxpsnwo ? nrm2[i] : maxpsnwo);
599 std::cout <<
"(" << maxpsnw <<
"," << maxpsnwo <<
")" << std::endl;
602 for (
int i=0; i<nq; i++) {
606 for (
int i=0; i<nq; i++) {
607 MVT::MvTimesMatAddMv(-ONE,*Q[i],*newC[i],ONE,X);
610 MVT::MvScale(X,normX);
613 OPT::Apply(*(this->_Op),X,*MX);
614 this->_OpCounter += MVT::GetNumberVecs(X);
622 MagnitudeType maxNorm = ZERO;
623 for (
int j=0; j<xc; j++) {
624 MagnitudeType sum = ZERO;
625 for (
int k=0; k<nq; k++) {
626 for (
int i=0; i<qcs[k]; i++) {
627 sum += SCT::magnitude((*newC[k])(i,j))*SCT::magnitude((*newC[k])(i,j));
630 maxNorm = (sum > maxNorm) ? sum : maxNorm;
634 if (maxNorm < 0.36) {
635 doGramSchmidt =
false;
639 for (
int k=0; k<nq; k++) {
640 for (
int j=0; j<xc; j++) {
641 for (
int i=0; i<qcs[k]; i++) {
642 (*newC[k])(i,j) *= normX[j];
650 for (
int i=0; i<nq; i++) {
652 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
657 for (
int i=0; i<nq; i++) {
664 doGramSchmidt =
false;
672 MagnitudeType condT = tolerance;
674 while (condT >= tolerance) {
682 std::vector<MagnitudeType> Dh(xc), Dhi(xc);
683 for (
int i=0; i<xc; i++) {
684 Dh[i] = SCT::magnitude(SCT::squareroot(XtMX(i,i)));
685 Dhi[i] = (Dh[i] == ZERO ? ZERO : MONE/Dh[i]);
688 for (
int i=0; i<xc; i++) {
689 for (
int j=0; j<xc; j++) {
690 XtMX(i,j) *= Dhi[i]*Dhi[j];
696 lapack.
HEEV(
'V',
'U', xc, XtMX.values(), XtMX.stride(), &lambda[0], &work[0], work.size(), &rwork[0], &info);
697 std::stringstream os;
698 os <<
"Anasazi::SVQBOrthoManager::findBasis(): Error code " << info <<
" from HEEV";
702 std::cout << dbgstr <<
"eigenvalues of XtMX: (";
703 for (
int i=0; i<xc-1; i++) {
704 std::cout << lambda[i] <<
",";
706 std::cout << lambda[xc-1] <<
")" << std::endl;
711 MagnitudeType maxLambda = lambda[xc-1],
712 minLambda = lambda[0];
714 for (
int i=0; i<xc; i++) {
715 if (lambda[i] <= 10*eps_*maxLambda) {
727 lambda[i] = SCTM::squareroot(lambda[i]);
728 lambdahi[i] = MONE/lambda[i];
735 std::vector<int> ind(xc);
736 for (
int i=0; i<xc; i++) {ind[i] = i;}
737 MVT::SetBlock(X,ind,*workX);
741 for (
int j=0; j<xc; j++) {
742 for (
int i=0; i<xc; i++) {
743 workU(i,j) *= Dhi[i]*lambdahi[j];
747 MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,X);
753 if (maxLambda >= tolerance * minLambda) {
755 OPT::Apply(*(this->_Op),X,*MX);
756 this->_OpCounter += MVT::GetNumberVecs(X);
761 MVT::SetBlock(*MX,ind,*workX);
764 MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,*MX);
770 for (
int j=0; j<xc; j++) {
771 for (
int i=0; i<xc; i++) {
772 workU(i,j) = Dh[i] * (*B)(i,j);
776 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
777 for (
int j=0; j<xc ;j++) {
778 for (
int i=0; i<xc; i++) {
779 (*B)(i,j) *= lambda[i];
786 std::cout << dbgstr <<
"augmenting multivec with " << iZeroMax+1 <<
" random directions" << std::endl;
791 std::vector<int> ind2(iZeroMax+1);
792 for (
int i=0; i<iZeroMax+1; i++) {
796 Xnull = MVT::CloneViewNonConst(X,ind2);
797 MVT::MvRandom(*Xnull);
799 MXnull = MVT::CloneViewNonConst(*MX,ind2);
800 OPT::Apply(*(this->_Op),*Xnull,*MXnull);
801 this->_OpCounter += MVT::GetNumberVecs(*Xnull);
802 MXnull = Teuchos::null;
804 Xnull = Teuchos::null;
806 doGramSchmidt =
true;
810 condT = SCTM::magnitude(maxLambda / minLambda);
812 std::cout << dbgstr <<
"condT: " << condT <<
", tolerance = " << tolerance << std::endl;
816 if ((doGramSchmidt ==
false) && (condT > SCTM::squareroot(tolerance))) {
817 doGramSchmidt =
true;
827 std::cout << dbgstr <<
"(numGS,numSVQB,numRand) : " 839 #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...
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
#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...
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...
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
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 ...
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.
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...