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();
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
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();
430 std::vector<int> qcs(nq);
431 for (
int i=0; i<nq; 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++) {
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;
514 if (MX == Teuchos::null) {
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);
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];
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++) {
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;}
741 for (
int j=0; j<xc; j++) {
742 for (
int i=0; i<xc; i++) {
743 workU(i,j) *= Dhi[i]*lambdahi[j];
753 if (maxLambda >= tolerance * minLambda) {
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;
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 static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvRandom(MV &mv)
Replace the vectors in mv with random vectors.
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const Teuchos::SerialDenseMatrix< int, ScalarType > &B, const ScalarType beta, MV &mv)
Update mv with .
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
void norm(const MV &X, std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > &normvec) const
Implements the interface OrthoManager::norm().
Declaration of basic traits for the multivector type.
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 ...
Virtual base class which defines basic traits for the operator type.
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.
static void Apply(const OP &Op, const MV &x, MV &y)
Application method which performs operation y = Op*x. An OperatorError exception is thrown if there i...
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...
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 ...
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
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...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
static void MvScale(MV &mv, const ScalarType alpha)
Scale each element of the vectors in mv with alpha.
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Creates a new MV that shares the selected contents of mv (shallow copy).
Exception thrown to signal error in an orthogonalization manager method.
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Copy the vectors in A to a set of vectors in mv indicated by the indices given in index...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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()
Destructor.
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...
SVQBOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, bool debug=false)
Constructor specifying re-orthogonalization tolerance.
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.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.