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" );
450 for (
int i=0; i<nq; i++) {
451 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength( *Q[i] ) != qr, std::invalid_argument,
452 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
453 qcs[i] = MVT::GetNumberVecs( *Q[i] );
454 TEUCHOS_TEST_FOR_EXCEPTION( qr < static_cast<ptrdiff_t>(qcs[i]), std::invalid_argument,
455 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
459 if ( C[i] == Teuchos::null ) {
460 C[i] = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
463 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
464 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
471 std::vector<ScalarType> oldDot( xc );
472 MVT::MvDot( X, *MX, oldDot );
473 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 474 *out <<
"oldDot = { ";
475 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
481 for (
int i=0; i<nq; i++) {
485 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 486 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
488 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
493 if (MQ[i] == Teuchos::null) {
494 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 495 *out <<
"Updating MX via M*X...\n";
497 OPT::Apply( *(this->_Op), X, *MX);
498 this->_OpCounter += MVT::GetNumberVecs(X);
501 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 502 *out <<
"Updating MX via M*Q...\n";
504 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
510 std::vector<ScalarType> newDot(xc);
511 MVT::MvDot( X, *MX, newDot );
512 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 513 *out <<
"newDot = { ";
514 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out,
" "));
519 for (
int j = 0; j < xc; ++j) {
521 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
522 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 523 *out <<
"kappa_*newDot[" <<j<<
"] == " << kappa_*newDot[j] <<
"... another step of Gram-Schmidt.\n";
525 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 526 Teuchos::TimeMonitor lcltimer( *timerReortho_ );
528 for (
int i=0; i<nq; i++) {
529 Teuchos::SerialDenseMatrix<int,ScalarType> C2(C[i]->numRows(), C[i]->numCols());
534 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 535 *out <<
"Applying projector P_Q[" << i <<
"]...\n";
537 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
541 if (MQ[i] == Teuchos::null) {
542 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 543 *out <<
"Updating MX via M*X...\n";
545 OPT::Apply( *(this->_Op), X, *MX);
546 this->_OpCounter += MVT::GetNumberVecs(X);
549 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 550 *out <<
"Updating MX via M*Q...\n";
552 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
560 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 561 *out <<
"Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
569 template<
class ScalarType,
class MV,
class OP>
572 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
573 Teuchos::RCP<MV> MX)
const {
577 int xc = MVT::GetNumberVecs(X);
578 ptrdiff_t xr = MVT::GetGlobalLength(X);
583 if (MX == Teuchos::null) {
585 MX = MVT::Clone(X,xc);
586 OPT::Apply(*(this->_Op),X,*MX);
587 this->_OpCounter += MVT::GetNumberVecs(X);
593 if ( B == Teuchos::null ) {
594 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
597 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
598 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
601 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
602 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
603 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
604 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
605 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
606 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
607 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
608 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
610 return findBasis(X, MX, *B,
true );
617 template<
class ScalarType,
class MV,
class OP>
620 Teuchos::Array<Teuchos::RCP<const MV> > Q,
621 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
622 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
624 Teuchos::Array<Teuchos::RCP<const MV> > MQ
627 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 629 Teuchos::RCP<Teuchos::FancyOStream>
630 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
631 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
632 *out <<
"Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
636 int xc = MVT::GetNumberVecs( X );
637 ptrdiff_t xr = MVT::GetGlobalLength( X );
643 if ( B == Teuchos::null ) {
644 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
649 if (MX == Teuchos::null) {
651 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 652 *out <<
"Allocating MX...\n";
654 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
655 OPT::Apply(*(this->_Op),X,*MX);
656 this->_OpCounter += MVT::GetNumberVecs(X);
661 MX = Teuchos::rcpFromRef(X);
664 int mxc = MVT::GetNumberVecs( *MX );
665 ptrdiff_t mxr = MVT::GetGlobalLength( *MX );
667 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
"Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
669 ptrdiff_t numbas = 0;
670 for (
int i=0; i<nq; i++) {
671 numbas += MVT::GetNumberVecs( *Q[i] );
675 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
676 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
678 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
679 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
681 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
682 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
684 TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
685 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
688 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 689 *out <<
"Orthogonalizing X against Q...\n";
691 projectMat(X,Q,C,MX,MQ);
693 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
699 int curxsize = xc - rank;
704 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 705 *out <<
"Attempting to find orthonormal basis for X...\n";
707 rank = findBasis(X,MX,*B,
false,curxsize);
709 if (oldrank != -1 && rank != oldrank) {
715 for (
int i=0; i<xc; i++) {
716 (*B)(i,oldrank) = oldCoeff(i,0);
721 if (rank != oldrank) {
729 for (
int i=0; i<xc; i++) {
730 oldCoeff(i,0) = (*B)(i,rank);
737 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 738 *out <<
"Finished computing basis.\n";
743 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank,
OrthoError,
744 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
746 if (rank != oldrank) {
762 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 763 *out <<
"Randomizing X[" << rank <<
"]...\n";
765 Teuchos::RCP<MV> curX, curMX;
766 std::vector<int> ind(1);
768 curX = MVT::CloneViewNonConst(X,ind);
769 MVT::MvRandom(*curX);
771 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 772 *out <<
"Applying operator to random vector.\n";
774 curMX = MVT::CloneViewNonConst(*MX,ind);
775 OPT::Apply( *(this->_Op), *curX, *curMX );
776 this->_OpCounter += MVT::GetNumberVecs(*curX);
784 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
785 projectMat(*curX,Q,dummyC,curMX,MQ);
791 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
792 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
794 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 795 *out <<
"Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
806 template<
class ScalarType,
class MV,
class OP>
808 MV &X, Teuchos::RCP<MV> MX,
809 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
810 bool completeBasis,
int howMany )
const {
827 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 829 Teuchos::RCP<Teuchos::FancyOStream>
830 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
831 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
832 *out <<
"Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
835 const ScalarType ONE = SCT::one();
836 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
838 int xc = MVT::GetNumberVecs( X );
847 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp ==
true && MX == Teuchos::null, std::logic_error,
848 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
849 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
850 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
855 int xstart = xc - howMany;
857 for (
int j = xstart; j < xc; j++) {
866 for (
int i=j+1; i<xc; ++i) {
871 std::vector<int> index(1);
873 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index );
874 Teuchos::RCP<MV> MXj;
875 if ((this->_hasOp)) {
877 MXj = MVT::CloneViewNonConst( *MX, index );
885 std::vector<int> prev_idx( numX );
886 Teuchos::RCP<const MV> prevX, prevMX;
889 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
890 prevX = MVT::CloneViewNonConst( X, prev_idx );
892 prevMX = MVT::CloneViewNonConst( *MX, prev_idx );
901 for (
int numTrials = 0; numTrials < 10; numTrials++) {
902 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 903 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
907 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
908 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
913 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
915 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 916 *out <<
"origNorm = " << origNorm[0] <<
"\n";
927 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 928 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
930 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
937 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 938 *out <<
"Updating MX[" << j <<
"]...\n";
940 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
945 MagnitudeType product_norm = product.normOne();
947 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 948 *out <<
"newNorm = " << newNorm[0] <<
"\n";
949 *out <<
"prodoct_norm = " << product_norm <<
"\n";
953 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
954 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 955 if (product_norm/newNorm[0] >= tol_) {
956 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
959 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
964 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
967 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 968 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
970 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
971 if ((this->_hasOp)) {
972 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 973 *out <<
"Updating MX[" << j <<
"]...\n";
975 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
979 product_norm = P2.normOne();
980 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 981 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
982 *out <<
"product_norm = " << product_norm <<
"\n";
984 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
986 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 987 if (product_norm/newNorm2[0] >= tol_) {
988 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
990 else if (newNorm[0] < newNorm2[0]) {
991 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
994 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
997 MVT::MvInit(*Xj,ZERO);
998 if ((this->_hasOp)) {
999 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1000 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
1002 MVT::MvInit(*MXj,ZERO);
1009 if (numTrials == 0) {
1010 for (
int i=0; i<numX; i++) {
1011 B(i,j) = product(i,0);
1017 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1018 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1019 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1023 MVT::MvScale( *Xj, ONE/newNorm[0]);
1025 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1026 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1029 MVT::MvScale( *MXj, ONE/newNorm[0]);
1033 if (numTrials == 0) {
1034 B(j,j) = newNorm[0];
1042 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1043 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1050 if (completeBasis) {
1052 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1053 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1055 MVT::MvRandom( *Xj );
1057 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1058 *out <<
"Updating M*X[" << j <<
"]...\n";
1060 OPT::Apply( *(this->_Op), *Xj, *MXj );
1061 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1072 if (rankDef ==
true) {
1073 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
1074 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
1075 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1076 *out <<
"Returning early, rank " << j <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1083 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1084 *out <<
"Returning " << xc <<
" from Anasazi::BasicOrthoManager::findBasis(...)\n";
1091 #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...