48#ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP
49#define ANASAZI_SVQB_ORTHOMANAGER_HPP
64#include "Teuchos_LAPACK.hpp"
68 template<
class ScalarType,
class MV,
class OP>
72 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
73 typedef Teuchos::ScalarTraits<ScalarType> SCT;
74 typedef Teuchos::ScalarTraits<MagnitudeType> SCTM;
85 SVQBOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null,
bool debug =
false );
138 Teuchos::Array<Teuchos::RCP<const MV> > Q,
139 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
140 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
141 Teuchos::RCP<MV> MX = Teuchos::null,
142 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
186 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
187 Teuchos::RCP<MV> MX = Teuchos::null
252 Teuchos::Array<Teuchos::RCP<const MV> > Q,
253 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
254 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
255 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
256 Teuchos::RCP<MV> MX = Teuchos::null,
257 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
269 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
276 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
280 Teuchos::RCP<const MV> MX = Teuchos::null,
281 Teuchos::RCP<const MV> MY = Teuchos::null
292 int findBasis(MV &X, Teuchos::RCP<MV> MX,
293 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
294 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
295 Teuchos::Array<Teuchos::RCP<const MV> > Q,
296 bool normalize_in )
const;
302 template<
class ScalarType,
class MV,
class OP>
304 :
MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(
" *** "), debug_(debug) {
306 Teuchos::LAPACK<int,MagnitudeType> lapack;
307 eps_ = lapack.LAMCH(
'E');
309 std::cout <<
"eps_ == " << eps_ << std::endl;
316 template<
class ScalarType,
class MV,
class OP>
317 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
319 const ScalarType ONE = SCT::one();
320 int rank = MVT::GetNumberVecs(X);
321 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
323 for (
int i=0; i<rank; i++) {
326 return xTx.normFrobenius();
331 template<
class ScalarType,
class MV,
class OP>
332 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
336 Teuchos::RCP<const MV> MX,
337 Teuchos::RCP<const MV> MY
339 int r1 = MVT::GetNumberVecs(X);
340 int r2 = MVT::GetNumberVecs(Y);
341 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
343 return xTx.normFrobenius();
349 template<
class ScalarType,
class MV,
class OP>
352 Teuchos::Array<Teuchos::RCP<const MV> > Q,
353 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
355 Teuchos::Array<Teuchos::RCP<const MV> > MQ)
const {
357 findBasis(X,MX,C,Teuchos::null,Q,
false);
364 template<
class ScalarType,
class MV,
class OP>
367 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
368 Teuchos::RCP<MV> MX)
const {
369 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C;
370 Teuchos::Array<Teuchos::RCP<const MV> > Q;
371 return findBasis(X,MX,C,B,Q,
true);
378 template<
class ScalarType,
class MV,
class OP>
381 Teuchos::Array<Teuchos::RCP<const MV> > Q,
382 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
383 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
385 Teuchos::Array<Teuchos::RCP<const MV> > MQ)
const {
387 return findBasis(X,MX,C,B,Q,
true);
419 template<
class ScalarType,
class MV,
class OP>
421 MV &X, Teuchos::RCP<MV> MX,
422 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
423 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
424 Teuchos::Array<Teuchos::RCP<const MV> > Q,
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) {
451 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
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)) {
462 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
463 "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
467 TEUCHOS_TEST_FOR_EXCEPTION( qsize != 0 && qr != xr , std::invalid_argument,
468 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
475 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > newC(nq);
477 for (
int i=0; i<nq; i++) {
479 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
480 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
481 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
482 "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
484 if ( C[i] == Teuchos::null ) {
485 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
488 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc, std::invalid_argument,
489 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
492 C[i]->putScalar(ZERO);
493 newC[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(C[i]->numRows(),C[i]->numCols()) );
502 if (normalize_in ==
true) {
503 if ( B == Teuchos::null ) {
504 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
506 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
507 "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
510 for (
int i=0; i<xc; i++) {
522 Teuchos::RCP<MV> workX;
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);
538 Teuchos::SerialDenseMatrix<int,ScalarType> XtMX(xc,xc), workU(1,1);
539 Teuchos::LAPACK<int,ScalarType> lapack;
544 std::vector<ScalarType> work;
545 std::vector<MagnitudeType> lambda, lambdahi, rwork;
548 int lwork = lapack.ILAENV(1,
"hetrd",
"VU",xc,-1,-1,-1);
551 TEUCHOS_TEST_FOR_EXCEPTION( lwork < 0, OrthoError,
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;
568 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
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++) {
664 info = C[i]->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,*newC[i],*B,ONE);
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";
712 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OrthoError,
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);
788 info = B->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,XtMX,workU,ZERO);
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++) {
808 Teuchos::RCP<MV> Xnull,MXnull;
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) : "
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator type.
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
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.
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.
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
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 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 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 ...
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...
~SVQBOrthoManager()
Destructor.
SVQBOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, bool debug=false)
Constructor specifying re-orthogonalization tolerance.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.