35 #ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP 36 #define ANASAZI_SVQB_ORTHOMANAGER_HPP 51 #include "Teuchos_LAPACK.hpp" 55 template<
class ScalarType,
class MV,
class OP>
59 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
60 typedef Teuchos::ScalarTraits<ScalarType> SCT;
61 typedef Teuchos::ScalarTraits<MagnitudeType> SCTM;
71 SVQBOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null,
bool debug =
false );
125 Teuchos::Array<Teuchos::RCP<const MV> > Q,
126 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
127 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
128 Teuchos::RCP<MV> MX = Teuchos::null,
129 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
173 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
174 Teuchos::RCP<MV> MX = Teuchos::null
239 Teuchos::Array<Teuchos::RCP<const MV> > Q,
240 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
241 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
242 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
243 Teuchos::RCP<MV> MX = Teuchos::null,
244 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
256 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
263 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
267 Teuchos::RCP<const MV> MX = Teuchos::null,
268 Teuchos::RCP<const MV> MY = Teuchos::null
279 int findBasis(MV &X, Teuchos::RCP<MV> MX,
280 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
281 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
282 Teuchos::Array<Teuchos::RCP<const MV> > Q,
283 bool normalize_in )
const;
289 template<
class ScalarType,
class MV,
class OP>
291 :
MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(
" *** "), debug_(debug) {
293 Teuchos::LAPACK<int,MagnitudeType> lapack;
294 eps_ = lapack.LAMCH(
'E');
296 std::cout <<
"eps_ == " << eps_ << std::endl;
303 template<
class ScalarType,
class MV,
class OP>
304 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
306 const ScalarType ONE = SCT::one();
307 int rank = MVT::GetNumberVecs(X);
308 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
310 for (
int i=0; i<rank; i++) {
313 return xTx.normFrobenius();
318 template<
class ScalarType,
class MV,
class OP>
319 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
323 Teuchos::RCP<const MV> MX,
324 Teuchos::RCP<const MV> MY
326 int r1 = MVT::GetNumberVecs(X);
327 int r2 = MVT::GetNumberVecs(Y);
328 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
330 return xTx.normFrobenius();
336 template<
class ScalarType,
class MV,
class OP>
339 Teuchos::Array<Teuchos::RCP<const MV> > Q,
340 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
342 Teuchos::Array<Teuchos::RCP<const MV> > MQ)
const {
344 findBasis(X,MX,C,Teuchos::null,Q,
false);
351 template<
class ScalarType,
class MV,
class OP>
354 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
355 Teuchos::RCP<MV> MX)
const {
356 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C;
357 Teuchos::Array<Teuchos::RCP<const MV> > Q;
358 return findBasis(X,MX,C,B,Q,
true);
365 template<
class ScalarType,
class MV,
class OP>
368 Teuchos::Array<Teuchos::RCP<const MV> > Q,
369 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
370 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
372 Teuchos::Array<Teuchos::RCP<const MV> > MQ)
const {
374 return findBasis(X,MX,C,B,Q,
true);
406 template<
class ScalarType,
class MV,
class OP>
408 MV &X, Teuchos::RCP<MV> MX,
409 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
410 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
411 Teuchos::Array<Teuchos::RCP<const MV> > Q,
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) {
438 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
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)) {
449 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
450 "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
454 TEUCHOS_TEST_FOR_EXCEPTION( qsize != 0 && qr != xr , std::invalid_argument,
455 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
462 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > newC(nq);
464 for (
int i=0; i<nq; i++) {
466 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
467 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
468 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
469 "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
471 if ( C[i] == Teuchos::null ) {
472 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
475 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc, std::invalid_argument,
476 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
479 C[i]->putScalar(ZERO);
480 newC[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(C[i]->numRows(),C[i]->numCols()) );
489 if (normalize_in ==
true) {
490 if ( B == Teuchos::null ) {
491 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
493 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
494 "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
497 for (
int i=0; i<xc; i++) {
509 Teuchos::RCP<MV> workX;
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);
525 Teuchos::SerialDenseMatrix<int,ScalarType> XtMX(xc,xc), workU(1,1);
526 Teuchos::LAPACK<int,ScalarType> lapack;
531 std::vector<ScalarType> work;
532 std::vector<MagnitudeType> lambda, lambdahi, rwork;
535 int lwork = lapack.ILAENV(1,
"hetrd",
"VU",xc,-1,-1,-1);
538 TEUCHOS_TEST_FOR_EXCEPTION( lwork < 0, OrthoError,
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;
555 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
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++) {
651 info = C[i]->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,*newC[i],*B,ONE);
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";
699 TEUCHOS_TEST_FOR_EXCEPTION( info != 0, OrthoError,
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);
775 info = B->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,XtMX,workU,ZERO);
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++) {
795 Teuchos::RCP<MV> Xnull,MXnull;
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 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...
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...
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 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...