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();
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
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();
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";
411 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 412 *out <<
"Allocating MX...\n";
414 if (MX == Teuchos::null) {
423 MX = Teuchos::rcpFromRef(X);
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++) {
439 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
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 );
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";
480 if (MQ[i] == Teuchos::null) {
481 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 482 *out <<
"Updating MX via M*X...\n";
488 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 489 *out <<
"Updating MX via M*Q...\n";
497 std::vector<ScalarType> newDot(xc);
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";
528 if (MQ[i] == Teuchos::null) {
529 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 530 *out <<
"Updating MX via M*X...\n";
536 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 537 *out <<
"Updating MX via M*Q...\n";
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 {
570 if (MX == Teuchos::null) {
580 if ( B == Teuchos::null ) {
581 B = Teuchos::rcp(
new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
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";
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";
648 MX = Teuchos::rcpFromRef(X);
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++) {
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";
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);
758 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 759 *out <<
"Applying operator to random vector.\n";
771 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
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());
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);
861 Teuchos::RCP<MV> MXj;
862 if ((this->_hasOp)) {
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;
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);
902 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 903 *out <<
"origNorm = " << origNorm[0] <<
"\n";
914 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 915 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
924 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 925 *out <<
"Updating MX[" << j <<
"]...\n";
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";
958 if ((this->_hasOp)) {
959 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 960 *out <<
"Updating MX[" << j <<
"]...\n";
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";
985 if ((this->_hasOp)) {
986 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 987 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
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";
1012 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1013 *out <<
"Normalizing M*X[" << j <<
"]...\n";
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";
1044 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 1045 *out <<
"Updating M*X[" << j <<
"]...\n";
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 static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Replace each element of the vectors in mv with alpha.
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 ...
Teuchos::ScalarTraits< ScalarType >::magnitudeType getKappa() const
Return parameter for re-orthogonalization threshold.
~BasicOrthoManager()
Destructor.
Declaration of basic traits for the multivector type.
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 MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Compute a vector b where the components are the individual dot-products of the i-th columns of A and ...
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...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
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 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...
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.
Traits class which defines basic operations on multivectors.
static Teuchos::RCP< MV > CloneCopy(const MV &mv)
Creates a new MV and copies contents of mv into the new vector (deep copy).
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...
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 ...
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.
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 ...
static int GetNumberVecs(const MV &mv)
Obtain the number of vectors in mv.
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)...
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.
static Teuchos::RCP< MV > Clone(const MV &mv, const int numvecs)
Creates a new empty MV containing numvecs columns.