47#ifndef ANASAZI_BASIC_ORTHOMANAGER_HPP
48#define ANASAZI_BASIC_ORTHOMANAGER_HPP
61#include "Teuchos_TimeMonitor.hpp"
62#include "Teuchos_LAPACK.hpp"
63#include "Teuchos_BLAS.hpp"
64#ifdef ANASAZI_BASIC_ORTHO_DEBUG
65# include <Teuchos_FancyOStream.hpp>
70 template<
class ScalarType,
class MV,
class OP>
74 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
75 typedef Teuchos::ScalarTraits<ScalarType> SCT;
85 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa = 1.41421356 ,
86 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
87 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
139 Teuchos::Array<Teuchos::RCP<const MV> > Q,
140 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
141 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
142 Teuchos::RCP<MV> MX = Teuchos::null,
143 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
187 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
188 Teuchos::RCP<MV> MX = Teuchos::null
253 Teuchos::Array<Teuchos::RCP<const MV> > Q,
254 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
255 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
256 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
257 Teuchos::RCP<MV> MX = Teuchos::null,
258 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
270 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
277 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
278 orthogErrorMat(
const MV &X1,
const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2)
const;
286 void setKappa(
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa ) { kappa_ = kappa; }
289 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
getKappa()
const {
return kappa_; }
295 MagnitudeType kappa_;
300 int findBasis(MV &X, Teuchos::RCP<MV> MX,
301 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
302 bool completeBasis,
int howMany = -1 )
const;
307 Teuchos::RCP<Teuchos::Time> timerReortho_;
314 template<
class ScalarType,
class MV,
class OP>
316 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa,
317 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
318 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol ) :
319 MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol)
320#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
321 , timerReortho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi::BasicOrthoManager::Re-orthogonalization"))
324 TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
325 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"eps\" must be non-negative.");
327 Teuchos::LAPACK<int,MagnitudeType> lapack;
328 eps_ = lapack.LAMCH(
'E');
329 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.75);
331 TEUCHOS_TEST_FOR_EXCEPTION(
332 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
333 std::invalid_argument,
334 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"tol\" must be in [0,1].");
341 template<
class ScalarType,
class MV,
class OP>
342 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
344 const ScalarType ONE = SCT::one();
345 int rank = MVT::GetNumberVecs(X);
346 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
348 for (
int i=0; i<rank; i++) {
351 return xTx.normFrobenius();
358 template<
class ScalarType,
class MV,
class OP>
359 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
361 int r1 = MVT::GetNumberVecs(X1);
362 int r2 = MVT::GetNumberVecs(X2);
363 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
365 return xTx.normFrobenius();
371 template<
class ScalarType,
class MV,
class OP>
374 Teuchos::Array<Teuchos::RCP<const MV> > Q,
375 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
377 Teuchos::Array<Teuchos::RCP<const MV> > MQ
394#ifdef ANASAZI_BASIC_ORTHO_DEBUG
396 Teuchos::RCP<Teuchos::FancyOStream>
397 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
398 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
399 *out <<
"Entering Anasazi::BasicOrthoManager::projectMat(...)\n";
402 ScalarType ONE = SCT::one();
404 int xc = MVT::GetNumberVecs( X );
405 ptrdiff_t xr = MVT::GetGlobalLength( X );
407 std::vector<int> qcs(nq);
409 if (nq == 0 || xc == 0 || xr == 0) {
410#ifdef ANASAZI_BASIC_ORTHO_DEBUG
411 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
415 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
424#ifdef ANASAZI_BASIC_ORTHO_DEBUG
425 *out <<
"Allocating MX...\n";
427 if (MX == Teuchos::null) {
429 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
430 OPT::Apply(*(this->_Op),X,*MX);
431 this->_OpCounter += MVT::GetNumberVecs(X);
436 MX = Teuchos::rcpFromRef(X);
438 int mxc = MVT::GetNumberVecs( *MX );
439 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
442 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
443 "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" );
445 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
446 "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" );
449 for (
int i=0; i<nq; i++) {
450 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
451 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
452 qcs[i] = MVT::GetNumberVecs( *Q[i] );
453 TEUCHOS_TEST_FOR_EXCEPTION( qr <
static_cast<ptrdiff_t
>(qcs[i]), std::invalid_argument,
454 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
456 if ( C[i] == Teuchos::null ) {
457 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
460 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
461 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
468 std::vector<ScalarType> oldDot( xc );
469 MVT::MvDot( X, *MX, oldDot );
470#ifdef ANASAZI_BASIC_ORTHO_DEBUG
471 *out <<
"oldDot = { ";
472 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
478 for (
int i=0; i<nq; i++) {
482#ifdef ANASAZI_BASIC_ORTHO_DEBUG
483 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
485 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
490 if (MQ[i] == Teuchos::null) {
491#ifdef ANASAZI_BASIC_ORTHO_DEBUG
492 *out <<
"Updating MX via M*X...\n";
494 OPT::Apply( *(this->_Op), X, *MX);
495 this->_OpCounter += MVT::GetNumberVecs(X);
498#ifdef ANASAZI_BASIC_ORTHO_DEBUG
499 *out <<
"Updating MX via M*Q...\n";
501 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
507 std::vector<ScalarType> newDot(xc);
508 MVT::MvDot( X, *MX, newDot );
509#ifdef ANASAZI_BASIC_ORTHO_DEBUG
510 *out <<
"newDot = { ";
511 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
516 for (
int j = 0; j < xc; ++j) {
518 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
519#ifdef ANASAZI_BASIC_ORTHO_DEBUG
520 *out <<
"kappa_*newDot[" <<j<<
"] == " << kappa_*newDot[j] <<
"... another step of Gram-Schmidt.\n";
522#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
523 Teuchos::TimeMonitor lcltimer( *timerReortho_ );
525 for (
int i=0; i<nq; i++) {
526 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
531#ifdef ANASAZI_BASIC_ORTHO_DEBUG
532 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
534 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
538 if (MQ[i] == Teuchos::null) {
539#ifdef ANASAZI_BASIC_ORTHO_DEBUG
540 *out <<
"Updating MX via M*X...\n";
542 OPT::Apply( *(this->_Op), X, *MX);
543 this->_OpCounter += MVT::GetNumberVecs(X);
546#ifdef ANASAZI_BASIC_ORTHO_DEBUG
547 *out <<
"Updating MX via M*Q...\n";
549 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
557#ifdef ANASAZI_BASIC_ORTHO_DEBUG
558 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
566 template<
class ScalarType,
class MV,
class OP>
569 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
570 Teuchos::RCP<MV> MX)
const {
574 int xc = MVT::GetNumberVecs(X);
575 ptrdiff_t xr = MVT::GetGlobalLength(X);
580 if (MX == Teuchos::null) {
582 MX = MVT::Clone(X,xc);
583 OPT::Apply(*(this->_Op),X,*MX);
584 this->_OpCounter += MVT::GetNumberVecs(X);
590 if ( B == Teuchos::null ) {
591 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
594 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
595 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
598 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
599 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
600 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
601 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
602 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
603 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
604 TEUCHOS_TEST_FOR_EXCEPTION(
static_cast<ptrdiff_t
>(xc) > xr, std::invalid_argument,
605 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
607 return findBasis(X, MX, *B,
true );
614 template<
class ScalarType,
class MV,
class OP>
617 Teuchos::Array<Teuchos::RCP<const MV> > Q,
618 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
619 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
621 Teuchos::Array<Teuchos::RCP<const MV> > MQ
624#ifdef ANASAZI_BASIC_ORTHO_DEBUG
626 Teuchos::RCP<Teuchos::FancyOStream>
627 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
628 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
629 *out <<
"Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
633 int xc = MVT::GetNumberVecs( X );
634 ptrdiff_t xr = MVT::GetGlobalLength( X );
640 if ( B == Teuchos::null ) {
641 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
646 if (MX == Teuchos::null) {
648#ifdef ANASAZI_BASIC_ORTHO_DEBUG
649 *out <<
"Allocating MX...\n";
651 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
652 OPT::Apply(*(this->_Op),X,*MX);
653 this->_OpCounter += MVT::GetNumberVecs(X);
658 MX = Teuchos::rcpFromRef(X);
661 int mxc = MVT::GetNumberVecs( *MX );
662 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
664 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
666 ptrdiff_t numbas = 0;
667 for (
int i=0; i<nq; i++) {
668 numbas += MVT::GetNumberVecs( *Q[i] );
672 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
673 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
675 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
676 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
678 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
679 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
681 TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
682 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
685#ifdef ANASAZI_BASIC_ORTHO_DEBUG
686 *out <<
"Orthogonalizing X against Q...\n";
688 projectMat(X,Q,C,MX,MQ);
690 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
696 int curxsize = xc - rank;
701#ifdef ANASAZI_BASIC_ORTHO_DEBUG
702 *out <<
"Attempting to find orthonormal basis for X...\n";
704 rank = findBasis(X,MX,*B,
false,curxsize);
706 if (oldrank != -1 && rank != oldrank) {
712 for (
int i=0; i<xc; i++) {
713 (*B)(i,oldrank) = oldCoeff(i,0);
718 if (rank != oldrank) {
726 for (
int i=0; i<xc; i++) {
727 oldCoeff(i,0) = (*B)(i,rank);
734#ifdef ANASAZI_BASIC_ORTHO_DEBUG
735 *out <<
"Finished computing basis.\n";
740 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank,
OrthoError,
741 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
743 if (rank != oldrank) {
759#ifdef ANASAZI_BASIC_ORTHO_DEBUG
760 *out <<
"Randomizing X[" << rank <<
"]...\n";
762 Teuchos::RCP<MV> curX, curMX;
763 std::vector<int> ind(1);
765 curX = MVT::CloneViewNonConst(X,ind);
766 MVT::MvRandom(*curX);
768#ifdef ANASAZI_BASIC_ORTHO_DEBUG
769 *out <<
"Applying operator to random vector.\n";
771 curMX = MVT::CloneViewNonConst(*MX,ind);
772 OPT::Apply( *(this->_Op), *curX, *curMX );
773 this->_OpCounter += MVT::GetNumberVecs(*curX);
781 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
782 projectMat(*curX,Q,dummyC,curMX,MQ);
788 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
789 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
791#ifdef ANASAZI_BASIC_ORTHO_DEBUG
792 *out <<
"Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
803 template<
class ScalarType,
class MV,
class OP>
805 MV &X, Teuchos::RCP<MV> MX,
806 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
807 bool completeBasis,
int howMany )
const {
824#ifdef ANASAZI_BASIC_ORTHO_DEBUG
826 Teuchos::RCP<Teuchos::FancyOStream>
827 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
828 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
829 *out <<
"Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
832 const ScalarType ONE = SCT::one();
833 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
835 int xc = MVT::GetNumberVecs( X );
844 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp ==
true && MX == Teuchos::null, std::logic_error,
845 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
846 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
847 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
852 int xstart = xc - howMany;
854 for (
int j = xstart; j < xc; j++) {
863 for (
int i=j+1; i<xc; ++i) {
868 std::vector<int> index(1);
870 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
871 Teuchos::RCP<MV> MXj;
872 if ((this->_hasOp)) {
874 MXj = MVT::CloneViewNonConst( *MX, index );
882 std::vector<int> prev_idx( numX );
883 Teuchos::RCP<const MV> prevX, prevMX;
886 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
887 prevX = MVT::CloneViewNonConst( X, prev_idx );
889 prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
898 for (
int numTrials = 0; numTrials < 10; numTrials++) {
899#ifdef ANASAZI_BASIC_ORTHO_DEBUG
900 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
904 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
905 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
910 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
912#ifdef ANASAZI_BASIC_ORTHO_DEBUG
913 *out <<
"origNorm = " << origNorm[0] <<
"\n";
924#ifdef ANASAZI_BASIC_ORTHO_DEBUG
925 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
927 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
934#ifdef ANASAZI_BASIC_ORTHO_DEBUG
935 *out <<
"Updating MX[" << j <<
"]...\n";
937 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
942 MagnitudeType product_norm = product.normOne();
944#ifdef ANASAZI_BASIC_ORTHO_DEBUG
945 *out <<
"newNorm = " << newNorm[0] <<
"\n";
946 *out <<
"prodoct_norm = " << product_norm <<
"\n";
950 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
951#ifdef ANASAZI_BASIC_ORTHO_DEBUG
952 if (product_norm/newNorm[0] >= tol_) {
953 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
956 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
961 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
964#ifdef ANASAZI_BASIC_ORTHO_DEBUG
965 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
967 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
968 if ((this->_hasOp)) {
969#ifdef ANASAZI_BASIC_ORTHO_DEBUG
970 *out <<
"Updating MX[" << j <<
"]...\n";
972 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
976 product_norm = P2.normOne();
977#ifdef ANASAZI_BASIC_ORTHO_DEBUG
978 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
979 *out <<
"product_norm = " << product_norm <<
"\n";
981 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
983#ifdef ANASAZI_BASIC_ORTHO_DEBUG
984 if (product_norm/newNorm2[0] >= tol_) {
985 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
987 else if (newNorm[0] < newNorm2[0]) {
988 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
991 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
994 MVT::MvInit(*Xj,ZERO);
995 if ((this->_hasOp)) {
996#ifdef ANASAZI_BASIC_ORTHO_DEBUG
997 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
999 MVT::MvInit(*MXj,ZERO);
1006 if (numTrials == 0) {
1007 for (
int i=0; i<numX; i++) {
1008 B(i,j) = product(i,0);
1014 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1015#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1016 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1020 MVT::MvScale( *Xj, ONE/newNorm[0]);
1022#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1023 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1026 MVT::MvScale( *MXj, ONE/newNorm[0]);
1030 if (numTrials == 0) {
1031 B(j,j) = newNorm[0];
1039#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1040 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1047 if (completeBasis) {
1049#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1050 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1052 MVT::MvRandom( *Xj );
1054#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1055 *out <<
"Updating M*X[" << j <<
"]...\n";
1057 OPT::Apply( *(this->_Op), *Xj, *MXj );
1058 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1069 if (rankDef ==
true) {
1070 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1071 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
1072#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1073 *out <<
"Returning early, rank " << j <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1080#ifdef ANASAZI_BASIC_ORTHO_DEBUG
1081 *out <<
"Returning " << xc <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
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.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
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 getKappa() const
Return parameter for re-orthogonalization threshold.
~BasicOrthoManager()
Destructor.
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 setKappa(typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa)
Set parameter for re-orthogonalization threshold.
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 ...
Teuchos::ScalarTraits< ScalarType >::magnitudeType orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP< const MV > MX1, Teuchos::RCP< const MV > MX2) const
This method computes the error in orthogonality of two multivectors, measured as the Frobenius norm o...
BasicOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa=1.41421356, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType eps=0.0, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType tol=0.20)
Constructor specifying re-orthogonalization tolerance.
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.
Exception thrown to signal error in an orthogonalization manager method.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.