16#ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP
17#define ANASAZI_SVQB_ORTHOMANAGER_HPP
32#include "Teuchos_LAPACK.hpp"
36 template<
class ScalarType,
class MV,
class OP>
40 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
41 typedef Teuchos::ScalarTraits<ScalarType> SCT;
42 typedef Teuchos::ScalarTraits<MagnitudeType> SCTM;
53 SVQBOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null,
bool debug =
false );
106 Teuchos::Array<Teuchos::RCP<const MV> > Q,
107 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
108 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
109 Teuchos::RCP<MV> MX = Teuchos::null,
110 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
154 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
155 Teuchos::RCP<MV> MX = Teuchos::null
220 Teuchos::Array<Teuchos::RCP<const MV> > Q,
221 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
222 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
223 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
224 Teuchos::RCP<MV> MX = Teuchos::null,
225 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
237 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
244 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
248 Teuchos::RCP<const MV> MX = Teuchos::null,
249 Teuchos::RCP<const MV> MY = Teuchos::null
260 int findBasis(MV &X, Teuchos::RCP<MV> MX,
261 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
262 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
263 Teuchos::Array<Teuchos::RCP<const MV> > Q,
264 bool normalize_in )
const;
270 template<
class ScalarType,
class MV,
class OP>
272 :
MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(
" *** "), debug_(debug) {
274 Teuchos::LAPACK<int,MagnitudeType> lapack;
275 eps_ = lapack.LAMCH(
'E');
277 std::cout <<
"eps_ == " << eps_ << std::endl;
284 template<
class ScalarType,
class MV,
class OP>
285 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
287 const ScalarType ONE = SCT::one();
288 int rank = MVT::GetNumberVecs(X);
289 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
291 for (
int i=0; i<rank; i++) {
294 return xTx.normFrobenius();
299 template<
class ScalarType,
class MV,
class OP>
300 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
304 Teuchos::RCP<const MV> MX,
305 Teuchos::RCP<const MV> MY
307 int r1 = MVT::GetNumberVecs(X);
308 int r2 = MVT::GetNumberVecs(Y);
309 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
311 return xTx.normFrobenius();
317 template<
class ScalarType,
class MV,
class OP>
320 Teuchos::Array<Teuchos::RCP<const MV> > Q,
321 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
323 Teuchos::Array<Teuchos::RCP<const MV> > MQ)
const {
325 findBasis(X,MX,C,Teuchos::null,Q,
false);
332 template<
class ScalarType,
class MV,
class OP>
335 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
336 Teuchos::RCP<MV> MX)
const {
337 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C;
338 Teuchos::Array<Teuchos::RCP<const MV> > Q;
339 return findBasis(X,MX,C,B,Q,
true);
346 template<
class ScalarType,
class MV,
class OP>
349 Teuchos::Array<Teuchos::RCP<const MV> > Q,
350 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
351 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
353 Teuchos::Array<Teuchos::RCP<const MV> > MQ)
const {
355 return findBasis(X,MX,C,B,Q,
true);
387 template<
class ScalarType,
class MV,
class OP>
389 MV &X, Teuchos::RCP<MV> MX,
390 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
391 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
392 Teuchos::Array<Teuchos::RCP<const MV> > Q,
393 bool normalize_in)
const {
395 const ScalarType ONE = SCT::one();
396 const MagnitudeType MONE = SCTM::one();
397 const MagnitudeType ZERO = SCTM::zero();
404 int xc = MVT::GetNumberVecs(X);
405 ptrdiff_t xr = MVT::GetGlobalLength( X );
409 ptrdiff_t qr = (nq == 0) ? 0 : MVT::GetGlobalLength(*Q[0]);
411 std::vector<int> qcs(nq);
412 for (
int i=0; i<nq; i++) {
413 qcs[i] = MVT::GetNumberVecs(*Q[i]);
417 if (normalize_in ==
true && qsize + xc > xr) {
419 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
420 "Anasazi::SVQBOrthoManager::findBasis(): Orthogonalization constraints not feasible" );
424 if (normalize_in ==
false && (qsize == 0 || xc == 0)) {
428 else if (normalize_in ==
true && (xc == 0 || xr == 0)) {
430 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
431 "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
435 TEUCHOS_TEST_FOR_EXCEPTION( qsize != 0 && qr != xr , std::invalid_argument,
436 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
443 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > newC(nq);
445 for (
int i=0; i<nq; i++) {
447 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
448 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
449 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
450 "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
452 if ( C[i] == Teuchos::null ) {
453 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
456 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc, std::invalid_argument,
457 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
460 C[i]->putScalar(ZERO);
461 newC[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(C[i]->numRows(),C[i]->numCols()) );
470 if (normalize_in ==
true) {
471 if ( B == Teuchos::null ) {
472 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
474 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
475 "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
478 for (
int i=0; i<xc; i++) {
490 Teuchos::RCP<MV> workX;
492 workX = MVT::Clone(X,xc);
495 if (MX == Teuchos::null) {
497 MX = MVT::Clone(X,xc);
498 OPT::Apply(*(this->_Op),X,*MX);
499 this->_OpCounter += MVT::GetNumberVecs(X);
503 MX = Teuchos::rcpFromRef(X);
505 std::vector<ScalarType> normX(xc), invnormX(xc);
506 Teuchos::SerialDenseMatrix<int,ScalarType> XtMX(xc,xc), workU(1,1);
507 Teuchos::LAPACK<int,ScalarType> lapack;
512 std::vector<ScalarType> work;
513 std::vector<MagnitudeType> lambda, lambdahi, rwork;
516 int lwork = lapack.ILAENV(1,
"hetrd",
"VU",xc,-1,-1,-1);
519 TEUCHOS_TEST_FOR_EXCEPTION( lwork < 0, OrthoError,
520 "Anasazi::SVQBOrthoManager::findBasis(): Error code from ILAENV" );
522 lwork = (lwork+2)*xc;
525 lwork = (3*xc-2 > 1) ? 3*xc - 2 : 1;
530 workU.reshape(xc,xc);
534 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
535 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
536 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
537 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with MX" );
540 bool doGramSchmidt =
true;
542 MagnitudeType tolerance = MONE/SCTM::squareroot(eps_);
545 while (doGramSchmidt) {
555 std::vector<MagnitudeType> normX_mag(xc);
557 for (
int i=0; i<xc; ++i) {
558 normX[i] = normX_mag[i];
559 invnormX[i] = (normX_mag[i] == ZERO) ? ZERO : MONE/normX_mag[i];
563 MVT::MvScale(X,invnormX);
565 MVT::MvScale(*MX,invnormX);
569 std::vector<MagnitudeType> nrm2(xc);
570 std::cout << dbgstr <<
"max post-scale norm: (with/without MX) : ";
571 MagnitudeType maxpsnw = ZERO, maxpsnwo = ZERO;
573 for (
int i=0; i<xc; i++) {
574 maxpsnw = (nrm2[i] > maxpsnw ? nrm2[i] : maxpsnw);
577 for (
int i=0; i<xc; i++) {
578 maxpsnwo = (nrm2[i] > maxpsnwo ? nrm2[i] : maxpsnwo);
580 std::cout <<
"(" << maxpsnw <<
"," << maxpsnwo <<
")" << std::endl;
583 for (
int i=0; i<nq; i++) {
587 for (
int i=0; i<nq; i++) {
588 MVT::MvTimesMatAddMv(-ONE,*Q[i],*newC[i],ONE,X);
591 MVT::MvScale(X,normX);
594 OPT::Apply(*(this->_Op),X,*MX);
595 this->_OpCounter += MVT::GetNumberVecs(X);
603 MagnitudeType maxNorm = ZERO;
604 for (
int j=0; j<xc; j++) {
605 MagnitudeType sum = ZERO;
606 for (
int k=0; k<nq; k++) {
607 for (
int i=0; i<qcs[k]; i++) {
608 sum += SCT::magnitude((*newC[k])(i,j))*SCT::magnitude((*newC[k])(i,j));
611 maxNorm = (sum > maxNorm) ? sum : maxNorm;
615 if (maxNorm < 0.36) {
616 doGramSchmidt =
false;
620 for (
int k=0; k<nq; k++) {
621 for (
int j=0; j<xc; j++) {
622 for (
int i=0; i<qcs[k]; i++) {
623 (*newC[k])(i,j) *= normX[j];
631 for (
int i=0; i<nq; i++) {
632 info = C[i]->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,*newC[i],*B,ONE);
633 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
638 for (
int i=0; i<nq; i++) {
645 doGramSchmidt =
false;
653 MagnitudeType condT = tolerance;
655 while (condT >= tolerance) {
663 std::vector<MagnitudeType> Dh(xc), Dhi(xc);
664 for (
int i=0; i<xc; i++) {
665 Dh[i] = SCT::magnitude(SCT::squareroot(XtMX(i,i)));
666 Dhi[i] = (Dh[i] == ZERO ? ZERO : MONE/Dh[i]);
669 for (
int i=0; i<xc; i++) {
670 for (
int j=0; j<xc; j++) {
671 XtMX(i,j) *= Dhi[i]*Dhi[j];
677 lapack.HEEV(
'V',
'U', xc, XtMX.values(), XtMX.stride(), &lambda[0], &work[0], work.size(), &rwork[0], &info);
678 std::stringstream os;
679 os <<
"Anasazi::SVQBOrthoManager::findBasis(): Error code " << info <<
" from HEEV";
680 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OrthoError,
683 std::cout << dbgstr <<
"eigenvalues of XtMX: (";
684 for (
int i=0; i<xc-1; i++) {
685 std::cout << lambda[i] <<
",";
687 std::cout << lambda[xc-1] <<
")" << std::endl;
692 MagnitudeType maxLambda = lambda[xc-1],
693 minLambda = lambda[0];
695 for (
int i=0; i<xc; i++) {
696 if (lambda[i] <= 10*eps_*maxLambda) {
708 lambda[i] = SCTM::squareroot(lambda[i]);
709 lambdahi[i] = MONE/lambda[i];
716 std::vector<int> ind(xc);
717 for (
int i=0; i<xc; i++) {ind[i] = i;}
718 MVT::SetBlock(X,ind,*workX);
722 for (
int j=0; j<xc; j++) {
723 for (
int i=0; i<xc; i++) {
724 workU(i,j) *= Dhi[i]*lambdahi[j];
728 MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,X);
734 if (maxLambda >= tolerance * minLambda) {
736 OPT::Apply(*(this->_Op),X,*MX);
737 this->_OpCounter += MVT::GetNumberVecs(X);
742 MVT::SetBlock(*MX,ind,*workX);
745 MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,*MX);
751 for (
int j=0; j<xc; j++) {
752 for (
int i=0; i<xc; i++) {
753 workU(i,j) = Dh[i] * (*B)(i,j);
756 info = B->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,XtMX,workU,ZERO);
757 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error,
"Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
758 for (
int j=0; j<xc ;j++) {
759 for (
int i=0; i<xc; i++) {
760 (*B)(i,j) *= lambda[i];
767 std::cout << dbgstr <<
"augmenting multivec with " << iZeroMax+1 <<
" random directions" << std::endl;
772 std::vector<int> ind2(iZeroMax+1);
773 for (
int i=0; i<iZeroMax+1; i++) {
776 Teuchos::RCP<MV> Xnull,MXnull;
777 Xnull = MVT::CloneViewNonConst(X,ind2);
778 MVT::MvRandom(*Xnull);
780 MXnull = MVT::CloneViewNonConst(*MX,ind2);
781 OPT::Apply(*(this->_Op),*Xnull,*MXnull);
782 this->_OpCounter += MVT::GetNumberVecs(*Xnull);
783 MXnull = Teuchos::null;
785 Xnull = Teuchos::null;
787 doGramSchmidt =
true;
791 condT = SCTM::magnitude(maxLambda / minLambda);
793 std::cout << dbgstr <<
"condT: " << condT <<
", tolerance = " << tolerance << std::endl;
797 if ((doGramSchmidt ==
false) && (condT > SCTM::squareroot(tolerance))) {
798 doGramSchmidt =
true;
808 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.