00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00034 #ifndef ANASAZI_BASIC_ORTHOMANAGER_HPP
00035 #define ANASAZI_BASIC_ORTHOMANAGER_HPP
00036
00044 #include "AnasaziConfigDefs.hpp"
00045 #include "AnasaziMultiVecTraits.hpp"
00046 #include "AnasaziOperatorTraits.hpp"
00047 #include "AnasaziMatOrthoManager.hpp"
00048 #include "Teuchos_TimeMonitor.hpp"
00049 #include "Teuchos_LAPACK.hpp"
00050 #include "Teuchos_BLAS.hpp"
00051 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00052 # include <Teuchos_FancyOStream.hpp>
00053 #endif
00054
00055 namespace Anasazi {
00056
00057 template<class ScalarType, class MV, class OP>
00058 class BasicOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00059
00060 private:
00061 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00062 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00063 typedef MultiVecTraits<ScalarType,MV> MVT;
00064 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00065
00066 public:
00067
00069
00070
00071 BasicOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null,
00072 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa = 1.41421356 ,
00073 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
00074 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
00075
00076
00078 ~BasicOrthoManager() {}
00080
00081
00083
00084
00085
00124 void projectMat (
00125 MV &X,
00126 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00127 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00128 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00129 Teuchos::RCP<MV> MX = Teuchos::null,
00130 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00131 ) const;
00132
00133
00172 int normalizeMat (
00173 MV &X,
00174 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00175 Teuchos::RCP<MV> MX = Teuchos::null
00176 ) const;
00177
00178
00238 int projectAndNormalizeMat (
00239 MV &X,
00240 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00241 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00242 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00243 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00244 Teuchos::RCP<MV> MX = Teuchos::null,
00245 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00246 ) const;
00247
00249
00251
00252
00257 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00258 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
00259
00264 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00265 orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const;
00266
00268
00270
00271
00273 void setKappa( typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa ) { kappa_ = kappa; }
00274
00276 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType getKappa() const { return kappa_; }
00277
00279
00280 private:
00282 MagnitudeType kappa_;
00283 MagnitudeType eps_;
00284 MagnitudeType tol_;
00285
00286
00287 int findBasis(MV &X, Teuchos::RCP<MV> MX,
00288 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
00289 bool completeBasis, int howMany = -1 ) const;
00290
00291
00292
00293
00294 Teuchos::RCP<Teuchos::Time> timerReortho_;
00295
00296 };
00297
00298
00300
00301 template<class ScalarType, class MV, class OP>
00302 BasicOrthoManager<ScalarType,MV,OP>::BasicOrthoManager( Teuchos::RCP<const OP> Op,
00303 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa,
00304 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
00305 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol ) :
00306 MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol),
00307 timerReortho_(Teuchos::TimeMonitor::getNewTimer("Anasazi::BasicOrthoManager::Re-orthogonalization"))
00308 {
00309 TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
00310 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative.");
00311 if (eps_ == 0) {
00312 Teuchos::LAPACK<int,MagnitudeType> lapack;
00313 eps_ = lapack.LAMCH('E');
00314 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.75);
00315 }
00316 TEST_FOR_EXCEPTION(
00317 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
00318 std::invalid_argument,
00319 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1].");
00320 }
00321
00322
00323
00325
00326 template<class ScalarType, class MV, class OP>
00327 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00328 BasicOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
00329 const ScalarType ONE = SCT::one();
00330 int rank = MVT::GetNumberVecs(X);
00331 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00332 innerProdMat(X,X,xTx,MX,MX);
00333 for (int i=0; i<rank; i++) {
00334 xTx(i,i) -= ONE;
00335 }
00336 return xTx.normFrobenius();
00337 }
00338
00339
00340
00342
00343 template<class ScalarType, class MV, class OP>
00344 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00345 BasicOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const {
00346 int r1 = MVT::GetNumberVecs(X1);
00347 int r2 = MVT::GetNumberVecs(X2);
00348 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
00349 innerProdMat(X1,X2,xTx,MX1,MX2);
00350 return xTx.normFrobenius();
00351 }
00352
00353
00354
00356 template<class ScalarType, class MV, class OP>
00357 void BasicOrthoManager<ScalarType, MV, OP>::projectMat(
00358 MV &X,
00359 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00360 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00361 Teuchos::RCP<MV> MX,
00362 Teuchos::Array<Teuchos::RCP<const MV> > MQ
00363 ) const {
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00380
00381 Teuchos::RCP<Teuchos::FancyOStream>
00382 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00383 out->setShowAllFrontMatter(false).setShowProcRank(true);
00384 *out << "Entering Anasazi::BasicOrthoManager::projectMat(...)\n";
00385 #endif
00386
00387 ScalarType ONE = SCT::one();
00388
00389 int xc = MVT::GetNumberVecs( X );
00390 int xr = MVT::GetVecLength( X );
00391 int nq = Q.length();
00392 std::vector<int> qcs(nq);
00393
00394 if (nq == 0 || xc == 0 || xr == 0) {
00395 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00396 *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
00397 #endif
00398 return;
00399 }
00400 int qr = MVT::GetVecLength ( *Q[0] );
00401
00402
00403
00404 C.resize(nq);
00405
00406
00407
00408 if (this->_hasOp) {
00409 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00410 *out << "Allocating MX...\n";
00411 #endif
00412 if (MX == Teuchos::null) {
00413
00414 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00415 OPT::Apply(*(this->_Op),X,*MX);
00416 this->_OpCounter += MVT::GetNumberVecs(X);
00417 }
00418 }
00419 else {
00420
00421 MX = Teuchos::rcpFromRef(X);
00422 }
00423 int mxc = MVT::GetNumberVecs( *MX );
00424 int mxr = MVT::GetVecLength( *MX );
00425
00426
00427 TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
00428 "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" );
00429
00430 TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument,
00431 "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" );
00432
00433
00434 int baslen = 0;
00435 for (int i=0; i<nq; i++) {
00436 TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument,
00437 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" );
00438 qcs[i] = MVT::GetNumberVecs( *Q[i] );
00439 TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
00440 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" );
00441 baslen += qcs[i];
00442
00443
00444 if ( C[i] == Teuchos::null ) {
00445 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00446 }
00447 else {
00448 TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument,
00449 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" );
00450 }
00451 }
00452
00453
00454
00455
00456 std::vector<ScalarType> oldDot( xc );
00457 MVT::MvDot( X, *MX, oldDot );
00458 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00459 *out << "oldDot = { ";
00460 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
00461 *out << "}\n";
00462 #endif
00463
00464 MQ.resize(nq);
00465
00466 for (int i=0; i<nq; i++) {
00467
00468 innerProdMat(*Q[i],X,*C[i],MQ[i],MX);
00469
00470 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00471 *out << "Applying projector P_Q[" << i << "]...\n";
00472 #endif
00473 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X );
00474
00475
00476
00477 if (this->_hasOp) {
00478 if (MQ[i] == Teuchos::null) {
00479 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00480 *out << "Updating MX via M*X...\n";
00481 #endif
00482 OPT::Apply( *(this->_Op), X, *MX);
00483 this->_OpCounter += MVT::GetNumberVecs(X);
00484 }
00485 else {
00486 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00487 *out << "Updating MX via M*Q...\n";
00488 #endif
00489 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX );
00490 }
00491 }
00492 }
00493
00494
00495 std::vector<ScalarType> newDot(xc);
00496 MVT::MvDot( X, *MX, newDot );
00497 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00498 *out << "newDot = { ";
00499 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out, " "));
00500 *out << "}\n";
00501 #endif
00502
00503
00504 for (int j = 0; j < xc; ++j) {
00505
00506 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) {
00507 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00508 *out << "kappa_*newDot[" <<j<< "] == " << kappa_*newDot[j] << "... another step of Gram-Schmidt.\n";
00509 #endif
00510 Teuchos::TimeMonitor lcltimer( *timerReortho_ );
00511 for (int i=0; i<nq; i++) {
00512 Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]);
00513
00514
00515 innerProdMat(*Q[i],X,C2,MQ[i],MX);
00516 *C[i] += C2;
00517 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00518 *out << "Applying projector P_Q[" << i << "]...\n";
00519 #endif
00520 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X );
00521
00522
00523 if (this->_hasOp) {
00524 if (MQ[i] == Teuchos::null) {
00525 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00526 *out << "Updating MX via M*X...\n";
00527 #endif
00528 OPT::Apply( *(this->_Op), X, *MX);
00529 this->_OpCounter += MVT::GetNumberVecs(X);
00530 }
00531 else {
00532 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00533 *out << "Updating MX via M*Q...\n";
00534 #endif
00535 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX );
00536 }
00537 }
00538 }
00539 break;
00540 }
00541 }
00542
00543 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00544 *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n";
00545 #endif
00546 }
00547
00548
00549
00551
00552 template<class ScalarType, class MV, class OP>
00553 int BasicOrthoManager<ScalarType, MV, OP>::normalizeMat(
00554 MV &X,
00555 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00556 Teuchos::RCP<MV> MX) const {
00557
00558
00559
00560 int xc = MVT::GetNumberVecs(X);
00561 int xr = MVT::GetVecLength(X);
00562
00563
00564
00565 if (this->_hasOp) {
00566 if (MX == Teuchos::null) {
00567
00568 MX = MVT::Clone(X,xc);
00569 OPT::Apply(*(this->_Op),X,*MX);
00570 this->_OpCounter += MVT::GetNumberVecs(X);
00571 }
00572 }
00573
00574
00575
00576 if ( B == Teuchos::null ) {
00577 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00578 }
00579
00580 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00581 int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX ) : xr;
00582
00583
00584 TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
00585 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" );
00586 TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
00587 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
00588 TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
00589 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
00590 TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
00591 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
00592
00593 return findBasis(X, MX, *B, true );
00594 }
00595
00596
00597
00599
00600 template<class ScalarType, class MV, class OP>
00601 int BasicOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat(
00602 MV &X,
00603 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00604 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00605 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00606 Teuchos::RCP<MV> MX,
00607 Teuchos::Array<Teuchos::RCP<const MV> > MQ
00608 ) const {
00609
00610 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00611
00612 Teuchos::RCP<Teuchos::FancyOStream>
00613 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00614 out->setShowAllFrontMatter(false).setShowProcRank(true);
00615 *out << "Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
00616 #endif
00617
00618 int nq = Q.length();
00619 int xc = MVT::GetNumberVecs( X );
00620 int xr = MVT::GetVecLength( X );
00621 int rank;
00622
00623
00624
00625
00626 if ( B == Teuchos::null ) {
00627 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00628 }
00629
00630
00631 if (this->_hasOp) {
00632 if (MX == Teuchos::null) {
00633
00634 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00635 *out << "Allocating MX...\n";
00636 #endif
00637 MX = MVT::Clone(X,MVT::GetNumberVecs(X));
00638 OPT::Apply(*(this->_Op),X,*MX);
00639 this->_OpCounter += MVT::GetNumberVecs(X);
00640 }
00641 }
00642 else {
00643
00644 MX = Teuchos::rcpFromRef(X);
00645 }
00646
00647 int mxc = MVT::GetNumberVecs( *MX );
00648 int mxr = MVT::GetVecLength( *MX );
00649
00650 TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" );
00651
00652 int numbas = 0;
00653 for (int i=0; i<nq; i++) {
00654 numbas += MVT::GetNumberVecs( *Q[i] );
00655 }
00656
00657
00658 TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
00659 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" );
00660
00661 TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument,
00662 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" );
00663
00664 TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument,
00665 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" );
00666
00667 TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument,
00668 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" );
00669
00670
00671 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00672 *out << "Orthogonalizing X against Q...\n";
00673 #endif
00674 projectMat(X,Q,C,MX,MQ);
00675
00676 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1);
00677
00678 rank = 0;
00679 int numTries = 10;
00680 int oldrank = -1;
00681 do {
00682 int curxsize = xc - rank;
00683
00684
00685
00686
00687 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00688 *out << "Attempting to find orthonormal basis for X...\n";
00689 #endif
00690 rank = findBasis(X,MX,*B,false,curxsize);
00691
00692 if (oldrank != -1 && rank != oldrank) {
00693
00694
00695
00696
00697
00698 for (int i=0; i<xc; i++) {
00699 (*B)(i,oldrank) = oldCoeff(i,0);
00700 }
00701 }
00702
00703 if (rank < xc) {
00704 if (rank != oldrank) {
00705
00706
00707
00708
00709
00710
00711
00712 for (int i=0; i<xc; i++) {
00713 oldCoeff(i,0) = (*B)(i,rank);
00714 }
00715 }
00716 }
00717
00718 if (rank == xc) {
00719
00720 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00721 *out << "Finished computing basis.\n";
00722 #endif
00723 break;
00724 }
00725 else {
00726 TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,
00727 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen");
00728
00729 if (rank != oldrank) {
00730
00731 numTries = 10;
00732
00733 oldrank = rank;
00734 }
00735 else {
00736
00737 if (numTries <= 0) {
00738 break;
00739 }
00740 }
00741
00742 numTries--;
00743
00744
00745 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00746 *out << "Randomizing X[" << rank << "]...\n";
00747 #endif
00748 Teuchos::RCP<MV> curX, curMX;
00749 std::vector<int> ind(1);
00750 ind[0] = rank;
00751 curX = MVT::CloneView(X,ind);
00752 MVT::MvRandom(*curX);
00753 if (this->_hasOp) {
00754 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00755 *out << "Applying operator to random vector.\n";
00756 #endif
00757 curMX = MVT::CloneView(*MX,ind);
00758 OPT::Apply( *(this->_Op), *curX, *curMX );
00759 this->_OpCounter += MVT::GetNumberVecs(*curX);
00760 }
00761
00762
00763
00764
00765
00766 {
00767 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
00768 projectMat(*curX,Q,dummyC,curMX,MQ);
00769 }
00770 }
00771 } while (1);
00772
00773
00774 TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error,
00775 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." );
00776
00777 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00778 *out << "Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n";
00779 #endif
00780
00781 return rank;
00782 }
00783
00784
00785
00787
00788
00789 template<class ScalarType, class MV, class OP>
00790 int BasicOrthoManager<ScalarType, MV, OP>::findBasis(
00791 MV &X, Teuchos::RCP<MV> MX,
00792 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
00793 bool completeBasis, int howMany ) const {
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00811
00812 Teuchos::RCP<Teuchos::FancyOStream>
00813 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00814 out->setShowAllFrontMatter(false).setShowProcRank(true);
00815 *out << "Entering Anasazi::BasicOrthoManager::findBasis(...)\n";
00816 #endif
00817
00818 const ScalarType ONE = SCT::one();
00819 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00820
00821 int xc = MVT::GetNumberVecs( X );
00822
00823 if (howMany == -1) {
00824 howMany = xc;
00825 }
00826
00827
00828
00829
00830 TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error,
00831 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS.");
00832 TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
00833 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" );
00834
00835
00836
00837
00838 int xstart = xc - howMany;
00839
00840 for (int j = xstart; j < xc; j++) {
00841
00842
00843 int numX = j;
00844
00845
00846
00847
00848
00849 for (int i=j+1; i<xc; ++i) {
00850 B(i,j) = ZERO;
00851 }
00852
00853
00854 std::vector<int> index(1);
00855 index[0] = j;
00856 Teuchos::RCP<MV> Xj = MVT::CloneView( X, index );
00857 Teuchos::RCP<MV> MXj;
00858 if ((this->_hasOp)) {
00859
00860 MXj = MVT::CloneView( *MX, index );
00861 }
00862 else {
00863
00864 MXj = Xj;
00865 }
00866
00867
00868 std::vector<int> prev_idx( numX );
00869 Teuchos::RCP<const MV> prevX, prevMX;
00870
00871 if (numX > 0) {
00872 for (int i=0; i<numX; ++i) prev_idx[i] = i;
00873 prevX = MVT::CloneView( X, prev_idx );
00874 if (this->_hasOp) {
00875 prevMX = MVT::CloneView( *MX, prev_idx );
00876 }
00877 }
00878
00879 bool rankDef = true;
00880
00881
00882
00883
00884 for (int numTrials = 0; numTrials < 10; numTrials++) {
00885 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00886 *out << "Trial " << numTrials << " for vector " << j << "\n";
00887 #endif
00888
00889
00890 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
00891 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
00892
00893
00894
00895
00896 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
00897 normMat(*Xj,origNorm,MXj);
00898 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00899 *out << "origNorm = " << origNorm[0] << "\n";
00900 #endif
00901
00902 if (numX > 0) {
00903
00904
00905
00906 innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj);
00907
00908
00909
00910 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00911 *out << "Orthogonalizing X[" << j << "]...\n";
00912 #endif
00913 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
00914
00915
00916 if (this->_hasOp) {
00917
00918
00919
00920 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00921 *out << "Updating MX[" << j << "]...\n";
00922 #endif
00923 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
00924 }
00925
00926
00927 normMat(*Xj,newNorm,MXj);
00928 MagnitudeType product_norm = product.normOne();
00929
00930 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00931 *out << "newNorm = " << newNorm[0] << "\n";
00932 *out << "prodoct_norm = " << product_norm << "\n";
00933 #endif
00934
00935
00936 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
00937 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00938 if (product_norm/newNorm[0] >= tol_) {
00939 *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n";
00940 }
00941 else {
00942 *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n";
00943 }
00944 #endif
00945
00946
00947 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
00948 innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj);
00949 product += P2;
00950 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00951 *out << "Orthogonalizing X[" << j << "]...\n";
00952 #endif
00953 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
00954 if ((this->_hasOp)) {
00955 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00956 *out << "Updating MX[" << j << "]...\n";
00957 #endif
00958 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
00959 }
00960
00961 normMat(*Xj,newNorm2,MXj);
00962 product_norm = P2.normOne();
00963 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00964 *out << "newNorm2 = " << newNorm2[0] << "\n";
00965 *out << "product_norm = " << product_norm << "\n";
00966 #endif
00967 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
00968
00969 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00970 if (product_norm/newNorm2[0] >= tol_) {
00971 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n";
00972 }
00973 else if (newNorm[0] < newNorm2[0]) {
00974 *out << "newNorm2 > newNorm... setting vector to zero.\n";
00975 }
00976 else {
00977 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n";
00978 }
00979 #endif
00980 MVT::MvInit(*Xj,ZERO);
00981 if ((this->_hasOp)) {
00982 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
00983 *out << "Setting MX[" << j << "] to zero as well.\n";
00984 #endif
00985 MVT::MvInit(*MXj,ZERO);
00986 }
00987 }
00988 }
00989 }
00990
00991
00992 if (numTrials == 0) {
00993 for (int i=0; i<numX; i++) {
00994 B(i,j) = product(i,0);
00995 }
00996 }
00997
00998
00999 normMat(*Xj,newNorm,MXj);
01000 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
01001 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01002 *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n";
01003 #endif
01004
01005
01006 MVT::MvScale( *Xj, ONE/newNorm[0]);
01007 if (this->_hasOp) {
01008 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01009 *out << "Normalizing M*X[" << j << "]...\n";
01010 #endif
01011
01012 MVT::MvScale( *MXj, ONE/newNorm[0]);
01013 }
01014
01015
01016 if (numTrials == 0) {
01017 B(j,j) = newNorm[0];
01018 }
01019
01020
01021 rankDef = false;
01022 break;
01023 }
01024 else {
01025 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01026 *out << "Not normalizing M*X[" << j << "]...\n";
01027 #endif
01028
01029
01030
01031 B(j,j) = ZERO;
01032
01033 if (completeBasis) {
01034
01035 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01036 *out << "Inserting random vector in X[" << j << "]...\n";
01037 #endif
01038 MVT::MvRandom( *Xj );
01039 if (this->_hasOp) {
01040 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01041 *out << "Updating M*X[" << j << "]...\n";
01042 #endif
01043 OPT::Apply( *(this->_Op), *Xj, *MXj );
01044 this->_OpCounter += MVT::GetNumberVecs(*Xj);
01045 }
01046 }
01047 else {
01048 rankDef = true;
01049 break;
01050 }
01051 }
01052 }
01053
01054
01055 if (rankDef == true) {
01056 TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
01057 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" );
01058 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01059 *out << "Returning early, rank " << j << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
01060 #endif
01061 return j;
01062 }
01063
01064 }
01065
01066 #ifdef ANASAZI_BASIC_ORTHO_DEBUG
01067 *out << "Returning " << xc << " from Anasazi::BasicOrthoManager::findBasis(...)\n";
01068 #endif
01069 return xc;
01070 }
01071
01072 }
01073
01074 #endif // ANASAZI_BASIC_ORTHOMANAGER_HPP
01075