34 #ifndef ANASAZI_BASIC_ORTHOMANAGER_HPP 35 #define ANASAZI_BASIC_ORTHOMANAGER_HPP 48 #include "Teuchos_TimeMonitor.hpp" 49 #include "Teuchos_LAPACK.hpp" 50 #include "Teuchos_BLAS.hpp" 51 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 52 # include <Teuchos_FancyOStream.hpp> 57 template<
class ScalarType,
class MV,
class OP>
61 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
62 typedef Teuchos::ScalarTraits<ScalarType> SCT;
72 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa = 1.41421356 ,
73 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
74 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
126 Teuchos::Array<Teuchos::RCP<const MV> > Q,
127 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
128 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
129 Teuchos::RCP<MV> MX = Teuchos::null,
130 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
174 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
175 Teuchos::RCP<MV> MX = Teuchos::null
240 Teuchos::Array<Teuchos::RCP<const MV> > Q,
241 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
242 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
243 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
244 Teuchos::RCP<MV> MX = Teuchos::null,
245 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
257 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
264 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
265 orthogErrorMat(
const MV &X1,
const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2)
const;
273 void setKappa(
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa ) { kappa_ = kappa; }
276 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
getKappa()
const {
return kappa_; }
282 MagnitudeType kappa_;
287 int findBasis(MV &X, Teuchos::RCP<MV> MX,
288 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
289 bool completeBasis,
int howMany = -1 )
const;
294 Teuchos::RCP<Teuchos::Time> timerReortho_;
301 template<
class ScalarType,
class MV,
class OP>
303 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa,
304 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
305 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol ) :
306 MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol)
307 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR
308 , timerReortho_(Teuchos::TimeMonitor::getNewTimer(
"Anasazi::BasicOrthoManager::Re-orthogonalization"))
311 TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
312 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"eps\" must be non-negative.");
314 Teuchos::LAPACK<int,MagnitudeType> lapack;
315 eps_ = lapack.LAMCH(
'E');
316 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.75);
318 TEUCHOS_TEST_FOR_EXCEPTION(
319 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
320 std::invalid_argument,
321 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"tol\" must be in [0,1].");
328 template<
class ScalarType,
class MV,
class OP>
329 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
331 const ScalarType ONE = SCT::one();
332 int rank = MVT::GetNumberVecs(X);
333 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
335 for (
int i=0; i<rank; i++) {
338 return xTx.normFrobenius();
345 template<
class ScalarType,
class MV,
class OP>
346 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
348 int r1 = MVT::GetNumberVecs(X1);
349 int r2 = MVT::GetNumberVecs(X2);
350 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
352 return xTx.normFrobenius();
358 template<
class ScalarType,
class MV,
class OP>
361 Teuchos::Array<Teuchos::RCP<const MV> > Q,
362 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
364 Teuchos::Array<Teuchos::RCP<const MV> > MQ
381 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 383 Teuchos::RCP<Teuchos::FancyOStream>
384 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
385 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
386 *out <<
"Entering Anasazi::BasicOrthoManager::projectMat(...)\n";
389 ScalarType ONE = SCT::one();
391 int xc = MVT::GetNumberVecs( X );
392 ptrdiff_t xr = MVT::GetGlobalLength( X );
394 std::vector<int> qcs(nq);
396 if (nq == 0 || xc == 0 || xr == 0) {
397 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 398 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
402 ptrdiff_t qr = MVT::GetGlobalLength ( *Q[0] );
411 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 412 *out <<
"Allocating MX...\n";
414 if (MX == Teuchos::null) {
416 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
417 OPT::Apply(*(this->_Op),X,*MX);
418 this->_OpCounter += MVT::GetNumberVecs(X);
423 MX = Teuchos::rcpFromRef(X);
425 int mxc = MVT::GetNumberVecs( *MX );
426 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
429 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
430 "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" );
432 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
433 "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" );
437 for (
int i=0; i<nq; i++) {
438 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
439 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
440 qcs[i] = MVT::GetNumberVecs( *Q[i] );
441 TEUCHOS_TEST_FOR_EXCEPTION( qr < static_cast<ptrdiff_t>(qcs[i]), std::invalid_argument,
442 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
446 if ( C[i] == Teuchos::null ) {
447 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
450 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
451 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
458 std::vector<ScalarType> oldDot( xc );
459 MVT::MvDot( X, *MX, oldDot );
460 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 461 *out <<
"oldDot = { ";
462 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
468 for (
int i=0; i<nq; i++) {
472 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 473 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
475 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
480 if (MQ[i] == Teuchos::null) {
481 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 482 *out <<
"Updating MX via M*X...\n";
484 OPT::Apply( *(this->_Op), X, *MX);
485 this->_OpCounter += MVT::GetNumberVecs(X);
488 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 489 *out <<
"Updating MX via M*Q...\n";
491 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
497 std::vector<ScalarType> newDot(xc);
498 MVT::MvDot( X, *MX, newDot );
499 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 500 *out <<
"newDot = { ";
501 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
506 for (
int j = 0; j < xc; ++j) {
508 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
509 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 510 *out <<
"kappa_*newDot[" <<j<<
"] == " << kappa_*newDot[j] <<
"... another step of Gram-Schmidt.\n";
512 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 513 Teuchos::TimeMonitor lcltimer( *timerReortho_ );
515 for (
int i=0; i<nq; i++) {
516 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
521 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 522 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
524 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
528 if (MQ[i] == Teuchos::null) {
529 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 530 *out <<
"Updating MX via M*X...\n";
532 OPT::Apply( *(this->_Op), X, *MX);
533 this->_OpCounter += MVT::GetNumberVecs(X);
536 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 537 *out <<
"Updating MX via M*Q...\n";
539 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
547 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 548 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
556 template<
class ScalarType,
class MV,
class OP>
559 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
560 Teuchos::RCP<MV> MX)
const {
564 int xc = MVT::GetNumberVecs(X);
565 ptrdiff_t xr = MVT::GetGlobalLength(X);
570 if (MX == Teuchos::null) {
572 MX = MVT::Clone(X,xc);
573 OPT::Apply(*(this->_Op),X,*MX);
574 this->_OpCounter += MVT::GetNumberVecs(X);
580 if ( B == Teuchos::null ) {
581 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
584 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
585 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
588 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
589 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
590 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
591 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
592 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
593 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
594 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
595 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
597 return findBasis(X, MX, *B,
true );
604 template<
class ScalarType,
class MV,
class OP>
607 Teuchos::Array<Teuchos::RCP<const MV> > Q,
608 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
609 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
611 Teuchos::Array<Teuchos::RCP<const MV> > MQ
614 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 616 Teuchos::RCP<Teuchos::FancyOStream>
617 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
618 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
619 *out <<
"Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
623 int xc = MVT::GetNumberVecs( X );
624 ptrdiff_t xr = MVT::GetGlobalLength( X );
630 if ( B == Teuchos::null ) {
631 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
636 if (MX == Teuchos::null) {
638 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 639 *out <<
"Allocating MX...\n";
641 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
642 OPT::Apply(*(this->_Op),X,*MX);
643 this->_OpCounter += MVT::GetNumberVecs(X);
648 MX = Teuchos::rcpFromRef(X);
651 int mxc = MVT::GetNumberVecs( *MX );
652 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
654 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
656 ptrdiff_t numbas = 0;
657 for (
int i=0; i<nq; i++) {
658 numbas += MVT::GetNumberVecs( *Q[i] );
662 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
663 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
665 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
666 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
668 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
669 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
671 TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
672 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
675 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 676 *out <<
"Orthogonalizing X against Q...\n";
678 projectMat(X,Q,C,MX,MQ);
680 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
686 int curxsize = xc - rank;
691 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 692 *out <<
"Attempting to find orthonormal basis for X...\n";
694 rank = findBasis(X,MX,*B,
false,curxsize);
696 if (oldrank != -1 && rank != oldrank) {
702 for (
int i=0; i<xc; i++) {
703 (*B)(i,oldrank) = oldCoeff(i,0);
708 if (rank != oldrank) {
716 for (
int i=0; i<xc; i++) {
717 oldCoeff(i,0) = (*B)(i,rank);
724 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 725 *out <<
"Finished computing basis.\n";
730 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank,
OrthoError,
731 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
733 if (rank != oldrank) {
749 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 750 *out <<
"Randomizing X[" << rank <<
"]...\n";
752 Teuchos::RCP<MV> curX, curMX;
753 std::vector<int> ind(1);
755 curX = MVT::CloneViewNonConst(X,ind);
756 MVT::MvRandom(*curX);
758 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 759 *out <<
"Applying operator to random vector.\n";
761 curMX = MVT::CloneViewNonConst(*MX,ind);
762 OPT::Apply( *(this->_Op), *curX, *curMX );
763 this->_OpCounter += MVT::GetNumberVecs(*curX);
771 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
772 projectMat(*curX,Q,dummyC,curMX,MQ);
778 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
779 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
781 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 782 *out <<
"Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
793 template<
class ScalarType,
class MV,
class OP>
795 MV &X, Teuchos::RCP<MV> MX,
796 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
797 bool completeBasis,
int howMany )
const {
814 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 816 Teuchos::RCP<Teuchos::FancyOStream>
817 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
818 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
819 *out <<
"Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
822 const ScalarType ONE = SCT::one();
823 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
825 int xc = MVT::GetNumberVecs( X );
834 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp ==
true && MX == Teuchos::null, std::logic_error,
835 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
836 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
837 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
842 int xstart = xc - howMany;
844 for (
int j = xstart; j < xc; j++) {
853 for (
int i=j+1; i<xc; ++i) {
858 std::vector<int> index(1);
860 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
861 Teuchos::RCP<MV> MXj;
862 if ((this->_hasOp)) {
864 MXj = MVT::CloneViewNonConst( *MX, index );
872 std::vector<int> prev_idx( numX );
873 Teuchos::RCP<const MV> prevX, prevMX;
876 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
877 prevX = MVT::CloneViewNonConst( X, prev_idx );
879 prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
888 for (
int numTrials = 0; numTrials < 10; numTrials++) {
889 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 890 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
894 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
895 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
900 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
902 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 903 *out <<
"origNorm = " << origNorm[0] <<
"\n";
914 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 915 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
917 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
924 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 925 *out <<
"Updating MX[" << j <<
"]...\n";
927 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
932 MagnitudeType product_norm = product.normOne();
934 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 935 *out <<
"newNorm = " << newNorm[0] <<
"\n";
936 *out <<
"prodoct_norm = " << product_norm <<
"\n";
940 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
941 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 942 if (product_norm/newNorm[0] >= tol_) {
943 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
946 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
951 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
954 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 955 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
957 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
958 if ((this->_hasOp)) {
959 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 960 *out <<
"Updating MX[" << j <<
"]...\n";
962 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
966 product_norm = P2.normOne();
967 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 968 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
969 *out <<
"product_norm = " << product_norm <<
"\n";
971 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
973 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 974 if (product_norm/newNorm2[0] >= tol_) {
975 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
977 else if (newNorm[0] < newNorm2[0]) {
978 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
981 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
984 MVT::MvInit(*Xj,ZERO);
985 if ((this->_hasOp)) {
986 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 987 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
989 MVT::MvInit(*MXj,ZERO);
996 if (numTrials == 0) {
997 for (
int i=0; i<numX; i++) {
998 B(i,j) = product(i,0);
1004 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1005 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1006 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1010 MVT::MvScale( *Xj, ONE/newNorm[0]);
1012 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1013 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1016 MVT::MvScale( *MXj, ONE/newNorm[0]);
1020 if (numTrials == 0) {
1021 B(j,j) = newNorm[0];
1029 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1030 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1037 if (completeBasis) {
1039 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1040 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1042 MVT::MvRandom( *Xj );
1044 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1045 *out <<
"Updating M*X[" << j <<
"]...\n";
1047 OPT::Apply( *(this->_Op), *Xj, *MXj );
1048 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1059 if (rankDef ==
true) {
1060 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1061 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
1062 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1063 *out <<
"Returning early, rank " << j <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1070 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1071 *out <<
"Returning " << xc <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1078 #endif // ANASAZI_BASIC_ORTHOMANAGER_HPP Teuchos::ScalarTraits< ScalarType >::magnitudeType getKappa() const
Return parameter for re-orthogonalization threshold.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
~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 ...
Declaration of basic traits for the multivector type.
Virtual base class which defines basic traits for the operator 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...
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 ...
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...
void setKappa(typename Teuchos::ScalarTraits< ScalarType >::magnitudeType kappa)
Set parameter for re-orthogonalization threshold.
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.
Exception thrown to signal error in an orthogonalization manager method.
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 ...
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.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
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...