15 #ifndef ANASAZI_ICSG_ORTHOMANAGER_HPP
16 #define ANASAZI_ICSG_ORTHOMANAGER_HPP
32 #ifdef ANASAZI_ICGS_DEBUG
33 # include <Teuchos_FancyOStream.hpp>
38 template<
class ScalarType,
class MV,
class OP>
351 numIters_ = numIters;
353 "Anasazi::ICGSOrthoManager::setNumIters(): input must be >= 1.");
371 bool completeBasis,
int howMany = -1)
const;
378 template<
class ScalarType,
class MV,
class OP>
387 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative.");
390 eps_ = lapack.
LAMCH(
'E');
395 std::invalid_argument,
396 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1].");
403 template<
class ScalarType,
class MV,
class OP>
406 const ScalarType ONE = SCT::one();
407 int rank = MVT::GetNumberVecs(X);
410 for (
int i=0; i<rank; i++) {
420 template<
class ScalarType,
class MV,
class OP>
423 int r1 = MVT::GetNumberVecs(X1);
424 int r2 = MVT::GetNumberVecs(X2);
433 template<
class ScalarType,
class MV,
class OP>
442 projectGen(X,Q,Q,
true,C,MX,MQ,MQ);
449 template<
class ScalarType,
class MV,
class OP>
458 int xc = MVT::GetNumberVecs(X);
459 ptrdiff_t xr = MVT::GetGlobalLength(X);
464 if (MX == Teuchos::null) {
466 MX = MVT::Clone(X,xc);
467 OPT::Apply(*(this->_Op),X,*MX);
468 this->_OpCounter += MVT::GetNumberVecs(X);
474 if ( B == Teuchos::null ) {
478 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
479 ptrdiff_t mxr = (this->_hasOp) ? MVT::GetGlobalLength( *MX ) : xr;
483 "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" );
485 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
487 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
488 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument,
489 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
491 return findBasis(X, MX, *B,
true );
498 template<
class ScalarType,
class MV,
class OP>
508 return projectAndNormalizeGen(X,Q,Q,
true,C,B,MX,MQ,MQ);
514 template<
class ScalarType,
class MV,
class OP>
541 #ifdef ANASAZI_ICGS_DEBUG
544 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
546 *out <<
"Entering Anasazi::ICGSOrthoManager::projectGen(...)\n";
549 const ScalarType ONE = SCT::one();
550 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
554 int sc = MVT::GetNumberVecs( S );
555 ptrdiff_t sr = MVT::GetGlobalLength( S );
556 int numxy = X.length();
558 "Anasazi::ICGSOrthoManager::projectGen(): X and Y must contain the same number of multivectors.");
559 std::vector<int> xyc(numxy);
561 if (numxy == 0 || sc == 0 || sr == 0) {
562 #ifdef ANASAZI_ICGS_DEBUG
563 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
578 "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." );
581 if (this->_hasOp ==
true) {
582 if (MS != Teuchos::null) {
584 "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." );
586 "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." );
591 ptrdiff_t sumxyc = 0;
594 for (
int i=0; i<numxy; i++) {
595 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
597 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"] length not consistent with S." );
599 "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i <<
"] length not consistent with S." );
601 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"] and Y[" << i <<
"] widths not consistent." );
603 xyc[i] = MVT::GetNumberVecs( *X[i] );
605 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i <<
"],Y[" << i <<
"] have less rows than columns, and therefore cannot be full rank." );
609 if ( C[i] == Teuchos::null ) {
614 "Anasazi::ICGSOrthoManager::projectGen(): Size of Q not consistent with size of C." );
618 if (MX[i] != Teuchos::null) {
619 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
620 "Anasazi::ICGSOrthoManager::projectGen(): Size of MX[" << i <<
"] not consistent with size of X[" << i <<
"]." );
625 if (MY[i] != Teuchos::null) {
626 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
627 "Anasazi::ICGSOrthoManager::projectGen(): Size of MY[" << i <<
"] not consistent with size of Y[" << i <<
"]." );
636 "Anasazi::ICGSOrthoManager::projectGen(): "
637 << (X[i] == Teuchos::null ?
"Y[" :
"X[") << i <<
"] was provided but "
638 << (X[i] == Teuchos::null ?
"X[" :
"Y[") << i <<
"] was not.");
644 "Anasazi::ICGSOrthoManager::projectGen(): dimension of all X[i],Y[i] is "
645 << sumxyc <<
", but length of vectors is only " << sr <<
". This is infeasible.");
680 if (MXmissing == 0) {
683 if (MYmissing == 0) {
690 switch (whichAlloc) {
706 if (MS == Teuchos::null) {
707 #ifdef ANASAZI_ICGS_DEBUG
708 *out <<
"Allocating MS...\n";
710 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
711 OPT::Apply(*(this->_Op),S,*MS);
712 this->_OpCounter += MVT::GetNumberVecs(S);
716 if (whichAlloc == 0) {
718 for (
int i=0; i<numxy; i++) {
719 if (MX[i] == Teuchos::null) {
720 #ifdef ANASAZI_ICGS_DEBUG
721 *out <<
"Allocating MX[" << i <<
"]...\n";
724 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
726 this->_OpCounter += xyc[i];
733 MS = Teuchos::rcpFromRef(S);
737 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
747 if (isBiortho ==
false) {
748 for (
int i=0; i<numxy; i++) {
749 #ifdef ANASAZI_ICGS_DEBUG
750 *out <<
"Computing YMX[" << i <<
"] and its Cholesky factorization...\n";
754 #ifdef ANASAZI_ICGS_DEBUG
761 MagnitudeType err = ZERO;
762 for (
int jj=0; jj<YMX[i]->numCols(); jj++) {
763 err =+ SCT::magnitude(SCT::imag((*YMX[i])(jj,jj)));
764 for (
int ii=jj; ii<YMX[i]->numRows(); ii++) {
765 err += SCT::magnitude( (*YMX[i])(ii,jj) - SCT::conjugate((*YMX[i])(jj,ii)) );
768 *out <<
"Symmetry error in YMX[" << i <<
"] == " << err <<
"\n";
773 lapack.
POTRF(
'U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info);
775 "Anasazi::ICGSOrthoManager::projectGen(): Error computing Cholesky factorization of Y[i]^T * M * X[i] using POTRF: returned info " << info);
780 #ifdef ANASAZI_ICGS_DEBUG
781 std::vector<MagnitudeType> oldNorms(sc);
783 *out <<
"oldNorms = { ";
784 std::copy(oldNorms.begin(), oldNorms.end(), std::ostream_iterator<MagnitudeType>(*out,
" "));
791 for (
int i=0; i<numxy; i++) {
792 C[i]->putScalar(ZERO);
793 Ccur[i].reshape(C[i]->numRows(),C[i]->numCols());
796 for (
int iter=0; iter < numIters_; iter++) {
797 #ifdef ANASAZI_ICGS_DEBUG
798 *out <<
"beginning iteration " << iter+1 <<
"\n";
802 for (
int i=0; i<numxy; i++) {
805 if (isBiortho ==
false) {
808 lapack.
POTRS(
'U',YMX[i]->numCols(),Ccur[i].numCols(),
809 YMX[i]->values(),YMX[i]->stride(),
810 Ccur[i].values(),Ccur[i].stride(), &info);
812 "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info <<
" from lapack::POTRS." );
816 #ifdef ANASAZI_ICGS_DEBUG
817 *out <<
"Applying projector P_{X[" << i <<
"],Y[" << i <<
"]}...\n";
819 MVT::MvTimesMatAddMv( -ONE, *X[i], Ccur[i], ONE, S );
826 #ifdef ANASAZI_ICGS_DEBUG
827 *out <<
"Updating MS...\n";
829 OPT::Apply( *(this->_Op), S, *MS);
830 this->_OpCounter += sc;
832 else if (updateMS == 2) {
833 #ifdef ANASAZI_ICGS_DEBUG
834 *out <<
"Updating MS...\n";
836 MVT::MvTimesMatAddMv( -ONE, *MX[i], Ccur[i], ONE, *MS );
841 #ifdef ANASAZI_ICGS_DEBUG
842 std::vector<MagnitudeType> newNorms(sc);
844 *out <<
"newNorms = { ";
845 std::copy(newNorms.begin(), newNorms.end(), std::ostream_iterator<MagnitudeType>(*out,
" "));
850 #ifdef ANASAZI_ICGS_DEBUG
851 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
859 template<
class ScalarType,
class MV,
class OP>
887 #ifdef ANASAZI_ICGS_DEBUG
890 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
892 *out <<
"Entering Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
896 int sc = MVT::GetNumberVecs( S );
897 ptrdiff_t sr = MVT::GetGlobalLength( S );
898 int numxy = X.length();
900 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X and Y must contain the same number of multivectors.");
901 std::vector<int> xyc(numxy);
903 if (sc == 0 || sr == 0) {
904 #ifdef ANASAZI_ICGS_DEBUG
905 *out <<
"Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
920 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." );
923 if (this->_hasOp ==
true) {
924 if (MS != Teuchos::null) {
926 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." );
928 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." );
933 ptrdiff_t sumxyc = 0;
936 for (
int i=0; i<numxy; i++) {
937 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
939 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"] length not consistent with S." );
941 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i <<
"] length not consistent with S." );
943 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"] and Y[" << i <<
"] widths not consistent." );
945 xyc[i] = MVT::GetNumberVecs( *X[i] );
947 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i <<
"],Y[" << i <<
"] have less rows than columns, and therefore cannot be full rank." );
951 if ( C[i] == Teuchos::null ) {
956 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of Q not consistent with size of C." );
960 if (MX[i] != Teuchos::null) {
961 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
962 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MX[" << i <<
"] not consistent with size of X[" << i <<
"]." );
967 if (MY[i] != Teuchos::null) {
968 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
969 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MY[" << i <<
"] not consistent with size of Y[" << i <<
"]." );
978 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): "
979 << (X[i] == Teuchos::null ?
"Y[" :
"X[") << i <<
"] was provided but "
980 << (X[i] == Teuchos::null ?
"X[" :
"Y[") << i <<
"] was not.");
986 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): dimension of all X[i],Y[i] is "
987 << sumxyc <<
" and requested " << sc <<
"-dimensional basis, but length of vectors is only "
988 << sr <<
". This is infeasible.");
1023 if (MXmissing == 0) {
1026 if (MYmissing == 0) {
1033 switch (whichAlloc) {
1049 if (MS == Teuchos::null) {
1050 #ifdef ANASAZI_ICGS_DEBUG
1051 *out <<
"Allocating MS...\n";
1053 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
1054 OPT::Apply(*(this->_Op),S,*MS);
1055 this->_OpCounter += MVT::GetNumberVecs(S);
1059 if (whichAlloc == 0) {
1061 for (
int i=0; i<numxy; i++) {
1062 if (MX[i] == Teuchos::null) {
1063 #ifdef ANASAZI_ICGS_DEBUG
1064 *out <<
"Allocating MX[" << i <<
"]...\n";
1067 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
1069 this->_OpCounter += xyc[i];
1076 MS = Teuchos::rcpFromRef(S);
1080 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
1085 if ( B == Teuchos::null ) {
1091 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of S must be consistent with size of B" );
1096 projectGen(S,X,Y,isBiortho,C,MS,MX,MY);
1104 int curssize = sc - rank;
1109 #ifdef ANASAZI_ICGS_DEBUG
1110 *out <<
"Attempting to find orthonormal basis for X...\n";
1112 rank = findBasis(S,MS,*B,
false,curssize);
1114 if (oldrank != -1 && rank != oldrank) {
1120 for (
int i=0; i<sc; i++) {
1121 (*B)(i,oldrank) = oldCoeff(i,0);
1126 if (rank != oldrank) {
1134 for (
int i=0; i<sc; i++) {
1135 oldCoeff(i,0) = (*B)(i,rank);
1142 #ifdef ANASAZI_ICGS_DEBUG
1143 *out <<
"Finished computing basis.\n";
1149 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen");
1151 if (rank != oldrank) {
1159 if (numTries <= 0) {
1167 #ifdef ANASAZI_ICGS_DEBUG
1168 *out <<
"Inserting random vector in X[" << rank <<
"]. Attempt " << 10-numTries <<
".\n";
1171 std::vector<int> ind(1);
1173 curS = MVT::CloneViewNonConst(S,ind);
1174 MVT::MvRandom(*curS);
1176 #ifdef ANASAZI_ICGS_DEBUG
1177 *out <<
"Applying operator to random vector.\n";
1179 curMS = MVT::CloneViewNonConst(*MS,ind);
1180 OPT::Apply( *(this->_Op), *curS, *curMS );
1181 this->_OpCounter += MVT::GetNumberVecs(*curS);
1190 projectGen(*curS,X,Y,isBiortho,dummyC,curMS,MX,MY);
1197 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Debug error in rank variable." );
1199 #ifdef ANASAZI_ICGS_DEBUG
1200 *out <<
"Returning " << rank <<
" from Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
1211 template<
class ScalarType,
class MV,
class OP>
1215 bool completeBasis,
int howMany )
const {
1232 #ifdef ANASAZI_ICGS_DEBUG
1235 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
1237 *out <<
"Entering Anasazi::ICGSOrthoManager::findBasis(...)\n";
1240 const ScalarType ONE = SCT::one();
1241 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
1243 int xc = MVT::GetNumberVecs( X );
1245 if (howMany == -1) {
1253 "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS.");
1255 "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
1260 int xstart = xc - howMany;
1262 for (
int j = xstart; j < xc; j++) {
1271 for (
int i=j+1; i<xc; ++i) {
1276 std::vector<int> index(1);
1280 if ((this->_hasOp)) {
1282 MXj = MVT::CloneViewNonConst( *MX, index );
1290 std::vector<int> prev_idx( numX );
1294 for (
int i=0; i<numX; ++i) prev_idx[i] = i;
1295 prevX = MVT::CloneView( X, prev_idx );
1297 prevMX = MVT::CloneView( *MX, prev_idx );
1301 bool rankDef =
true;
1306 for (
int numTrials = 0; numTrials < 10; numTrials++) {
1307 #ifdef ANASAZI_ICGS_DEBUG
1308 *out <<
"Trial " << numTrials <<
" for vector " << j <<
"\n";
1313 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
1320 #ifdef ANASAZI_ICGS_DEBUG
1321 *out <<
"origNorm = " << origNorm[0] <<
"\n";
1332 #ifdef ANASAZI_ICGS_DEBUG
1333 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
1335 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
1342 #ifdef ANASAZI_ICGS_DEBUG
1343 *out <<
"Updating MX[" << j <<
"]...\n";
1345 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
1350 MagnitudeType product_norm = product.normOne();
1352 #ifdef ANASAZI_ICGS_DEBUG
1353 *out <<
"newNorm = " << newNorm[0] <<
"\n";
1354 *out <<
"prodoct_norm = " << product_norm <<
"\n";
1358 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
1359 #ifdef ANASAZI_ICGS_DEBUG
1360 if (product_norm/newNorm[0] >= tol_) {
1361 *out <<
"product_norm/newNorm == " << product_norm/newNorm[0] <<
"... another step of Gram-Schmidt.\n";
1364 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... another step of Gram-Schmidt.\n";
1372 #ifdef ANASAZI_ICGS_DEBUG
1373 *out <<
"Orthogonalizing X[" << j <<
"]...\n";
1375 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
1376 if ((this->_hasOp)) {
1377 #ifdef ANASAZI_ICGS_DEBUG
1378 *out <<
"Updating MX[" << j <<
"]...\n";
1380 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
1384 product_norm = P2.normOne();
1385 #ifdef ANASAZI_ICGS_DEBUG
1386 *out <<
"newNorm2 = " << newNorm2[0] <<
"\n";
1387 *out <<
"product_norm = " << product_norm <<
"\n";
1389 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
1391 #ifdef ANASAZI_ICGS_DEBUG
1392 if (product_norm/newNorm2[0] >= tol_) {
1393 *out <<
"product_norm/newNorm2 == " << product_norm/newNorm2[0] <<
"... setting vector to zero.\n";
1395 else if (newNorm[0] < newNorm2[0]) {
1396 *out <<
"newNorm2 > newNorm... setting vector to zero.\n";
1399 *out <<
"eps*origNorm == " << eps_*origNorm[0] <<
"... setting vector to zero.\n";
1402 MVT::MvInit(*Xj,ZERO);
1403 if ((this->_hasOp)) {
1404 #ifdef ANASAZI_ICGS_DEBUG
1405 *out <<
"Setting MX[" << j <<
"] to zero as well.\n";
1407 MVT::MvInit(*MXj,ZERO);
1414 if (numTrials == 0) {
1415 for (
int i=0; i<numX; i++) {
1416 B(i,j) = product(i,0);
1422 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
1423 #ifdef ANASAZI_ICGS_DEBUG
1424 *out <<
"Normalizing X[" << j <<
"], norm(X[" << j <<
"]) = " << newNorm[0] <<
"\n";
1428 MVT::MvScale( *Xj, ONE/newNorm[0]);
1430 #ifdef ANASAZI_ICGS_DEBUG
1431 *out <<
"Normalizing M*X[" << j <<
"]...\n";
1434 MVT::MvScale( *MXj, ONE/newNorm[0]);
1438 if (numTrials == 0) {
1439 B(j,j) = newNorm[0];
1447 #ifdef ANASAZI_ICGS_DEBUG
1448 *out <<
"Not normalizing M*X[" << j <<
"]...\n";
1455 if (completeBasis) {
1457 #ifdef ANASAZI_ICGS_DEBUG
1458 *out <<
"Inserting random vector in X[" << j <<
"]...\n";
1460 MVT::MvRandom( *Xj );
1462 #ifdef ANASAZI_ICGS_DEBUG
1463 *out <<
"Updating M*X[" << j <<
"]...\n";
1465 OPT::Apply( *(this->_Op), *Xj, *MXj );
1466 this->_OpCounter += MVT::GetNumberVecs(*Xj);
1477 if (rankDef ==
true) {
1479 "Anasazi::ICGSOrthoManager::findBasis(): Unable to complete basis" );
1480 #ifdef ANASAZI_ICGS_DEBUG
1481 *out <<
"Returning early, rank " << j <<
" from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1488 #ifdef ANASAZI_ICGS_DEBUG
1489 *out <<
"Returning " << xc <<
" from Anasazi::ICGSOrthoManager::findBasis(...)\n";
1496 #endif // ANASAZI_ICSG_ORTHOMANAGER_HPP
Declaration of basic traits for the multivector type.
basic_FancyOStream & setShowProcRank(const bool showProcRank)
~ICGSOrthoManager()
Destructor.
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...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Virtual base class which defines basic traits for the operator type.
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
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.
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 ...
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.
ScalarTraits< ScalarType >::magnitudeType normFrobenius() const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void POTRF(const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
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 ...
Traits class which defines basic operations on multivectors.
Virtual base class which defines basic traits for the operator type.
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
Anasazi header file which uses auto-configuration information to include necessary C++ headers...
Exception thrown to signal error in an orthogonalization manager method.
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 ...
static magnitudeType magnitude(ScalarTypea)
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...
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 ...
int getNumIters() const
Return parameter for re-orthogonalization threshold.
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.
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...
void setNumIters(int numIters)
Set parameter for re-orthogonalization threshold.