34 #ifndef ANASAZI_ICSG_ORTHOMANAGER_HPP 35 #define ANASAZI_ICSG_ORTHOMANAGER_HPP 51 #ifdef ANASAZI_ICGS_DEBUG 52 # include <Teuchos_FancyOStream.hpp> 57 template<
class ScalarType,
class MV,
class OP>
370 numIters_ = numIters;
372 "Anasazi::ICGSOrthoManager::setNumIters(): input must be >= 1.");
390 bool completeBasis,
int howMany = -1)
const;
397 template<
class ScalarType,
class MV,
class OP>
406 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative.");
409 eps_ = lapack.
LAMCH(
'E');
414 std::invalid_argument,
415 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1].");
422 template<
class ScalarType,
class MV,
class OP>
425 const ScalarType ONE = SCT::one();
426 int rank = MVT::GetNumberVecs(X);
429 for (
int i=0; i<rank; i++) {
439 template<
class ScalarType,
class MV,
class OP>
442 int r1 = MVT::GetNumberVecs(X1);
443 int r2 = MVT::GetNumberVecs(X2);
452 template<
class ScalarType,
class MV,
class OP>
461 projectGen(X,Q,Q,
true,C,MX,MQ,MQ);
468 template<
class ScalarType,
class MV,
class OP>
477 int xc = MVT::GetNumberVecs(X);
478 ptrdiff_t xr = MVT::GetGlobalLength(X);
483 if (MX == Teuchos::null) {
485 MX = MVT::Clone(X,xc);
486 OPT::Apply(*(this->_Op),X,*MX);
487 this->_OpCounter += MVT::GetNumberVecs(X);
493 if ( B == Teuchos::null ) {
497 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
498 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
502 "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" );
504 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
506 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
508 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
510 return findBasis(X, MX, *B,
true );
517 template<
class ScalarType,
class MV,
class OP>
527 return projectAndNormalizeGen(X,Q,Q,
true,C,B,MX,MQ,MQ);
533 template<
class ScalarType,
class MV,
class OP>
560 #ifdef ANASAZI_ICGS_DEBUG 563 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
565 *out <<
"Entering Anasazi::ICGSOrthoManager::projectGen(...)\n";
568 const ScalarType ONE = SCT::one();
569 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
573 int sc = MVT::GetNumberVecs( S );
574 ptrdiff_t sr = MVT::GetGlobalLength( S );
575 int numxy = X.length();
577 "Anasazi::ICGSOrthoManager::projectGen(): X and Y must contain the same number of multivectors.");
578 std::vector<int> xyc(numxy);
580 if (numxy == 0 || sc == 0 || sr == 0) {
581 #ifdef ANASAZI_ICGS_DEBUG 582 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
597 "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." );
600 if (this->_hasOp ==
true) {
601 if (MS != Teuchos::null) {
603 "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." );
605 "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." );
610 ptrdiff_t sumxyc = 0;
613 for (
int i=0; i<numxy; i++) {
614 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
616 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"] length not consistent with S." );
618 "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i <<
"] length not consistent with S." );
620 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"] and Y[" << i <<
"] widths not consistent." );
622 xyc[i] = MVT::GetNumberVecs( *X[i] );
624 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"],Y[" << i <<
"] have less rows than columns, and therefore cannot be full rank." );
628 if ( C[i] == Teuchos::null ) {
633 "Anasazi::ICGSOrthoManager::projectGen(): Size of Q not consistent with size of C." );
637 if (MX[i] != Teuchos::null) {
638 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
639 "Anasazi::ICGSOrthoManager::projectGen(): Size of MX[" << i <<
"] not consistent with size of X[" << i <<
"]." );
644 if (MY[i] != Teuchos::null) {
645 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
646 "Anasazi::ICGSOrthoManager::projectGen(): Size of MY[" << i <<
"] not consistent with size of Y[" << i <<
"]." );
655 "Anasazi::ICGSOrthoManager::projectGen(): " 656 << (X[i] == Teuchos::null ?
"Y[" :
"X[") << i <<
"] was provided but " 657 << (X[i] == Teuchos::null ?
"X[" :
"Y[") << i <<
"] was not.");
663 "Anasazi::ICGSOrthoManager::projectGen(): dimension of all X[i],Y[i] is " 664 << sumxyc <<
", but length of vectors is only " << sr <<
". This is infeasible.");
699 if (MXmissing == 0) {
702 if (MYmissing == 0) {
709 switch (whichAlloc) {
725 if (MS == Teuchos::null) {
726 #ifdef ANASAZI_ICGS_DEBUG 727 *out <<
"Allocating MS...\n";
729 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
730 OPT::Apply(*(this->_Op),S,*MS);
731 this->_OpCounter += MVT::GetNumberVecs(S);
735 if (whichAlloc == 0) {
737 for (
int i=0; i<numxy; i++) {
738 if (MX[i] == Teuchos::null) {
739 #ifdef ANASAZI_ICGS_DEBUG 740 *out <<
"Allocating MX[" << i <<
"]...\n";
743 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
745 this->_OpCounter += xyc[i];
752 MS = Teuchos::rcpFromRef(S);
756 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
766 if (isBiortho ==
false) {
767 for (
int i=0; i<numxy; i++) {
768 #ifdef ANASAZI_ICGS_DEBUG 769 *out <<
"Computing YMX[" << i <<
"] and its Cholesky factorization...\n";
773 #ifdef ANASAZI_ICGS_DEBUG 780 MagnitudeType err = ZERO;
781 for (
int jj=0; jj<YMX[i]->numCols(); jj++) {
782 err =+ SCT::magnitude(SCT::imag((*YMX[i])(jj,jj)));
783 for (
int ii=jj; ii<YMX[i]->numRows(); ii++) {
784 err += SCT::magnitude( (*YMX[i])(ii,jj) - SCT::conjugate((*YMX[i])(jj,ii)) );
787 *out <<
"Symmetry error in YMX[" << i <<
"] == " << err <<
"\n";
792 lapack.
POTRF(
'U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info);
794 "Anasazi::ICGSOrthoManager::projectGen(): Error computing Cholesky factorization of Y[i]^T * M * X[i] using POTRF: returned info " << info);
799 #ifdef ANASAZI_ICGS_DEBUG 800 std::vector<MagnitudeType> oldNorms(sc);
802 *out <<
"oldNorms = { ";
803 std::copy(oldNorms.begin(), oldNorms.end(), std::ostream_iterator<MagnitudeType>(*out,
" "));
810 for (
int i=0; i<numxy; i++) {
811 C[i]->putScalar(ZERO);
812 Ccur[i].reshape(C[i]->numRows(),C[i]->numCols());
815 for (
int iter=0; iter < numIters_; iter++) {
816 #ifdef ANASAZI_ICGS_DEBUG 817 *out <<
"beginning iteration " << iter+1 <<
"\n";
821 for (
int i=0; i<numxy; i++) {
824 if (isBiortho ==
false) {
827 lapack.
POTRS(
'U',YMX[i]->numCols(),Ccur[i].numCols(),
828 YMX[i]->values(),YMX[i]->stride(),
829 Ccur[i].values(),Ccur[i].stride(), &info);
831 "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info <<
" from lapack::POTRS." );
835 #ifdef ANASAZI_ICGS_DEBUG 836 *out <<
"Applying projector P_{X[" << i <<
"],Y[" << i <<
"]}...\n";
838 MVT::MvTimesMatAddMv( -ONE, *X[i], Ccur[i], ONE, S );
845 #ifdef ANASAZI_ICGS_DEBUG 846 *out <<
"Updating MS...\n";
848 OPT::Apply( *(this->_Op), S, *MS);
849 this->_OpCounter += sc;
851 else if (updateMS == 2) {
852 #ifdef ANASAZI_ICGS_DEBUG 853 *out <<
"Updating MS...\n";
855 MVT::MvTimesMatAddMv( -ONE, *MX[i], Ccur[i], ONE, *MS );
860 #ifdef ANASAZI_ICGS_DEBUG 861 std::vector<MagnitudeType> newNorms(sc);
863 *out <<
"newNorms = { ";
864 std::copy(newNorms.begin(), newNorms.end(), std::ostream_iterator<MagnitudeType>(*out,
" "));
869 #ifdef ANASAZI_ICGS_DEBUG 870 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
878 template<
class ScalarType,
class MV,
class OP>
906 #ifdef ANASAZI_ICGS_DEBUG 909 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
911 *out <<
"Entering Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
915 int sc = MVT::GetNumberVecs( S );
916 ptrdiff_t sr = MVT::GetGlobalLength( S );
917 int numxy = X.length();
919 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X and Y must contain the same number of multivectors.");
920 std::vector<int> xyc(numxy);
922 if (sc == 0 || sr == 0) {
923 #ifdef ANASAZI_ICGS_DEBUG 924 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
939 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." );
942 if (this->_hasOp ==
true) {
943 if (MS != Teuchos::null) {
945 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." );
947 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." );
952 ptrdiff_t sumxyc = 0;
955 for (
int i=0; i<numxy; i++) {
956 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
958 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"] length not consistent with S." );
960 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i <<
"] length not consistent with S." );
962 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"] and Y[" << i <<
"] widths not consistent." );
964 xyc[i] = MVT::GetNumberVecs( *X[i] );
966 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"],Y[" << i <<
"] have less rows than columns, and therefore cannot be full rank." );
970 if ( C[i] == Teuchos::null ) {
975 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of Q not consistent with size of C." );
979 if (MX[i] != Teuchos::null) {
980 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
981 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MX[" << i <<
"] not consistent with size of X[" << i <<
"]." );
986 if (MY[i] != Teuchos::null) {
987 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
988 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MY[" << i <<
"] not consistent with size of Y[" << i <<
"]." );
997 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): " 998 << (X[i] == Teuchos::null ?
"Y[" :
"X[") << i <<
"] was provided but " 999 << (X[i] == Teuchos::null ?
"X[" :
"Y[") << i <<
"] was not.");
1005 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): dimension of all X[i],Y[i] is " 1006 << sumxyc <<
" and requested " << sc <<
"-dimensional basis, but length of vectors is only " 1007 << sr <<
". This is infeasible.");
1042 if (MXmissing == 0) {
1045 if (MYmissing == 0) {
1052 switch (whichAlloc) {
1068 if (MS == Teuchos::null) {
1069 #ifdef ANASAZI_ICGS_DEBUG 1070 *out <<
"Allocating MS...\n";
1072 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
1073 OPT::Apply(*(this->_Op),S,*MS);
1074 this->_OpCounter += MVT::GetNumberVecs(S);
1078 if (whichAlloc == 0) {
1080 for (
int i=0; i<numxy; i++) {
1081 if (MX[i] == Teuchos::null) {
1082 #ifdef ANASAZI_ICGS_DEBUG 1083 *out <<
"Allocating MX[" << i <<
"]...\n";
1086 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
1088 this->_OpCounter += xyc[i];
1095 MS = Teuchos::rcpFromRef(S);
1099 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
1104 if ( B == Teuchos::null ) {
1110 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of S must be consistent with size of B" );
1115 projectGen(S,X,Y,isBiortho,C,MS,MX,MY);
1123 int curssize = sc - rank;
1128 #ifdef ANASAZI_ICGS_DEBUG 1129 *out <<
"Attempting to find orthonormal basis for X...\n";
1131 rank = findBasis(S,MS,*B,
false,curssize);
1133 if (oldrank != -1 && rank != oldrank) {
1139 for (
int i=0; i<sc; i++) {
1140 (*B)(i,oldrank) = oldCoeff(i,0);
1145 if (rank != oldrank) {
1153 for (
int i=0; i<sc; i++) {
1154 oldCoeff(i,0) = (*B)(i,rank);
1161 #ifdef ANASAZI_ICGS_DEBUG 1162 *out <<
"Finished computing basis.\n";
1168 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen");
1170 if (rank != oldrank) {
1178 if (numTries <= 0) {
1186 #ifdef ANASAZI_ICGS_DEBUG 1187 *out <<
"Inserting random vector in X[" << rank <<
"]. Attempt " << 10-numTries <<
".\n";
1190 std::vector<int> ind(1);
1192 curS = MVT::CloneViewNonConst(S,ind);
1193 MVT::MvRandom(*curS);
1195 #ifdef ANASAZI_ICGS_DEBUG 1196 *out <<
"Applying operator to random vector.\n";
1198 curMS = MVT::CloneViewNonConst(*MS,ind);
1199 OPT::Apply( *(this->_Op), *curS, *curMS );
1200 this->_OpCounter += MVT::GetNumberVecs(*curS);
1209 projectGen(*curS,X,Y,isBiortho,dummyC,curMS,MX,MY);
1216 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Debug error in rank variable." );
1218 #ifdef ANASAZI_ICGS_DEBUG 1219 *out <<
"Returning " << rank <<
" from Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
1230 template<
class ScalarType,
class MV,
class OP>
1234 bool completeBasis,
int howMany )
const {
1251 #ifdef ANASAZI_ICGS_DEBUG 1254 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
1256 *out <<
"Entering Anasazi::ICGSOrthoManager::findBasis(...)\n";
1259 const ScalarType ONE = SCT::one();
1260 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
1262 int xc = MVT::GetNumberVecs( X );
1264 if (howMany == -1) {
1272 "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS.");
1274 "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
1279 int xstart = xc - howMany;
1281 for (
int j = xstart; j < xc; j++) {
1290 for (
int i=j+1; i<xc; ++i) {
1295 std::vector<int> index(1);
1299 if ((this->_hasOp)) {
1301 MXj = MVT::CloneViewNonConst( *MX, index );
1309 std::vector<int> prev_idx( numX );
1313 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
1314 prevX = MVT::CloneView( X, prev_idx );
1316 prevMX = MVT::CloneView( *MX, prev_idx );
1320 bool rankDef =
true;
1325 for (
int numTrials = 0; numTrials < 10; numTrials++) {
1326 #ifdef ANASAZI_ICGS_DEBUG 1327 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
1332 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
1339 #ifdef ANASAZI_ICGS_DEBUG 1340 *out <<
"origNorm = " << origNorm[0] <<
"\n";
1351 #ifdef ANASAZI_ICGS_DEBUG 1352 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
1354 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1361 #ifdef ANASAZI_ICGS_DEBUG 1362 *out <<
"Updating MX[" << j <<
"]...\n";
1364 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1369 MagnitudeType product_norm = product.normOne();
1371 #ifdef ANASAZI_ICGS_DEBUG 1372 *out <<
"newNorm = " << newNorm[0] <<
"\n";
1373 *out <<
"prodoct_norm = " << product_norm <<
"\n";
1377 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
1378 #ifdef ANASAZI_ICGS_DEBUG 1379 if (product_norm/newNorm[0] >= tol_) {
1380 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
1383 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
1391 #ifdef ANASAZI_ICGS_DEBUG 1392 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
1394 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1395 if ((this->_hasOp)) {
1396 #ifdef ANASAZI_ICGS_DEBUG 1397 *out <<
"Updating MX[" << j <<
"]...\n";
1399 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1403 product_norm = P2.normOne();
1404 #ifdef ANASAZI_ICGS_DEBUG 1405 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
1406 *out <<
"product_norm = " << product_norm <<
"\n";
1408 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
1410 #ifdef ANASAZI_ICGS_DEBUG 1411 if (product_norm/newNorm2[0] >= tol_) {
1412 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
1414 else if (newNorm[0] < newNorm2[0]) {
1415 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
1418 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
1421 MVT::MvInit(*Xj,ZERO);
1422 if ((this->_hasOp)) {
1423 #ifdef ANASAZI_ICGS_DEBUG 1424 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
1426 MVT::MvInit(*MXj,ZERO);
1433 if (numTrials == 0) {
1434 for (
int i=0; i<numX; i++) {
1435 B(i,j) = product(i,0);
1441 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1442 #ifdef ANASAZI_ICGS_DEBUG 1443 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1447 MVT::MvScale( *Xj, ONE/newNorm[0]);
1449 #ifdef ANASAZI_ICGS_DEBUG 1450 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1453 MVT::MvScale( *MXj, ONE/newNorm[0]);
1457 if (numTrials == 0) {
1458 B(j,j) = newNorm[0];
1466 #ifdef ANASAZI_ICGS_DEBUG 1467 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1474 if (completeBasis) {
1476 #ifdef ANASAZI_ICGS_DEBUG 1477 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1479 MVT::MvRandom( *Xj );
1481 #ifdef ANASAZI_ICGS_DEBUG 1482 *out <<
"Updating M*X[" << j <<
"]...\n";
1484 OPT::Apply( *(this->_Op), *Xj, *MXj );
1485 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1496 if (rankDef ==
true) {
1498 "Anasazi::ICGSOrthoManager::findBasis(): Unable to complete basis" );
1499 #ifdef ANASAZI_ICGS_DEBUG 1500 *out <<
"Returning early, rank " << j <<
" from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1507 #ifdef ANASAZI_ICGS_DEBUG 1508 *out <<
"Returning " << xc <<
" from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1515 #endif // ANASAZI_ICSG_ORTHOMANAGER_HPP ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
int getNumIters() const
Return parameter for re-orthogonalization threshold.
Declaration of basic traits for the multivector type.
basic_FancyOStream & setShowProcRank(const bool showProcRank)
~ICGSOrthoManager()
Destructor.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
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 ...
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
ICGSOrthoManager(Teuchos::RCP< const OP > Op=Teuchos::null, int numIters=2, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType eps=0.0, typename Teuchos::ScalarTraits< ScalarType >::magnitudeType tol=0.20)
Constructor specifying the operator defining the inner product as well as the number of orthogonaliza...
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
int projectAndNormalizeGen(MV &S, Teuchos::Array< Teuchos::RCP< const MV > > X, Teuchos::Array< Teuchos::RCP< const MV > > Y, bool isBiOrtho, 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 > MS=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MX=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null)), Teuchos::Array< Teuchos::RCP< const MV > > MY=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Applies a series of generic projectors and returns an orthonormal basis for the residual data...
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...
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
ScalarType LAMCH(const char CMACH) const
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.
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 ...
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
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.
static magnitudeType magnitude(ScalarType a)
void POTRS(const char UPLO, const OrdinalType n, const OrdinalType nrhs, const ScalarType *A, const OrdinalType lda, ScalarType *B, const OrdinalType ldb, OrdinalType *info) const
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...
void projectGen(MV &S, Teuchos::Array< Teuchos::RCP< const MV > > X, Teuchos::Array< Teuchos::RCP< const MV > > Y, bool isBiOrtho, Teuchos::Array< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > C=Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > >(Teuchos::null)), Teuchos::RCP< MV > MS=Teuchos::null, Teuchos::Array< Teuchos::RCP< const MV > > MX=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null)), Teuchos::Array< Teuchos::RCP< const MV > > MY=Teuchos::tuple(Teuchos::RCP< const MV >(Teuchos::null))) const
Applies a series of generic projectors.
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 ...
void POTRF(const char UPLO, const OrdinalType n, ScalarType *A, const OrdinalType lda, OrdinalType *info) const
void setNumIters(int numIters)
Set parameter for re-orthogonalization threshold.