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
00035 #ifndef ANASAZI_SVQB_ORTHOMANAGER_HPP
00036 #define ANASAZI_SVQB_ORTHOMANAGER_HPP
00037
00047 #include "AnasaziConfigDefs.hpp"
00048 #include "AnasaziMultiVecTraits.hpp"
00049 #include "AnasaziOperatorTraits.hpp"
00050 #include "AnasaziMatOrthoManager.hpp"
00051 #include "Teuchos_LAPACK.hpp"
00052
00053 namespace Anasazi {
00054
00055 template<class ScalarType, class MV, class OP>
00056 class SVQBOrthoManager : public MatOrthoManager<ScalarType,MV,OP> {
00057
00058 private:
00059 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00060 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00061 typedef Teuchos::ScalarTraits<MagnitudeType> SCTM;
00062 typedef MultiVecTraits<ScalarType,MV> MVT;
00063 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00064 std::string dbgstr;
00065
00066
00067 public:
00068
00070
00071
00072 SVQBOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, bool debug = false );
00073
00074
00076 ~SVQBOrthoManager() {};
00078
00079
00080
00082
00083
00084
00123 void projectMat (
00124 MV &X,
00125 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00126 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00127 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00128 Teuchos::RCP<MV> MX = Teuchos::null,
00129 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00130 ) const;
00131
00132
00171 int normalizeMat (
00172 MV &X,
00173 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00174 Teuchos::RCP<MV> MX = Teuchos::null
00175 ) const;
00176
00177
00237 int projectAndNormalizeMat (
00238 MV &X,
00239 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00240 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00241 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00242 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00243 Teuchos::RCP<MV> MX = Teuchos::null,
00244 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00245 ) const;
00246
00248
00250
00251
00256 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00257 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
00258
00263 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00264 orthogErrorMat(
00265 const MV &X,
00266 const MV &Y,
00267 Teuchos::RCP<const MV> MX = Teuchos::null,
00268 Teuchos::RCP<const MV> MY = Teuchos::null
00269 ) const;
00270
00272
00273 private:
00274
00275 MagnitudeType eps_;
00276 bool debug_;
00277
00278
00279 int findBasis(MV &X, Teuchos::RCP<MV> MX,
00280 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00281 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00282 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00283 bool normalize_in ) const;
00284 };
00285
00286
00288
00289 template<class ScalarType, class MV, class OP>
00290 SVQBOrthoManager<ScalarType,MV,OP>::SVQBOrthoManager( Teuchos::RCP<const OP> Op, bool debug)
00291 : MatOrthoManager<ScalarType,MV,OP>(Op), dbgstr(" *** "), debug_(debug) {
00292
00293 Teuchos::LAPACK<int,MagnitudeType> lapack;
00294 eps_ = lapack.LAMCH('E');
00295 if (debug_) {
00296 std::cout << "eps_ == " << eps_ << std::endl;
00297 }
00298 }
00299
00300
00302
00303 template<class ScalarType, class MV, class OP>
00304 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00305 SVQBOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
00306 const ScalarType ONE = SCT::one();
00307 int rank = MVT::GetNumberVecs(X);
00308 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00309 innerProdMat(X,X,xTx,MX,MX);
00310 for (int i=0; i<rank; i++) {
00311 xTx(i,i) -= ONE;
00312 }
00313 return xTx.normFrobenius();
00314 }
00315
00317
00318 template<class ScalarType, class MV, class OP>
00319 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00320 SVQBOrthoManager<ScalarType,MV,OP>::orthogErrorMat(
00321 const MV &X,
00322 const MV &Y,
00323 Teuchos::RCP<const MV> MX,
00324 Teuchos::RCP<const MV> MY
00325 ) const {
00326 int r1 = MVT::GetNumberVecs(X);
00327 int r2 = MVT::GetNumberVecs(Y);
00328 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
00329 innerProdMat(X,Y,xTx,MX,MY);
00330 return xTx.normFrobenius();
00331 }
00332
00333
00334
00336 template<class ScalarType, class MV, class OP>
00337 void SVQBOrthoManager<ScalarType, MV, OP>::projectMat(
00338 MV &X,
00339 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00340 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00341 Teuchos::RCP<MV> MX,
00342 Teuchos::Array<Teuchos::RCP<const MV> > MQ) const {
00343 (void)MQ;
00344 findBasis(X,MX,C,Teuchos::null,Q,false);
00345 }
00346
00347
00348
00350
00351 template<class ScalarType, class MV, class OP>
00352 int SVQBOrthoManager<ScalarType, MV, OP>::normalizeMat(
00353 MV &X,
00354 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00355 Teuchos::RCP<MV> MX) const {
00356 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C;
00357 Teuchos::Array<Teuchos::RCP<const MV> > Q;
00358 return findBasis(X,MX,C,B,Q,true);
00359 }
00360
00361
00362
00364
00365 template<class ScalarType, class MV, class OP>
00366 int SVQBOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat(
00367 MV &X,
00368 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00369 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00370 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00371 Teuchos::RCP<MV> MX,
00372 Teuchos::Array<Teuchos::RCP<const MV> > MQ) const {
00373 (void)MQ;
00374 return findBasis(X,MX,C,B,Q,true);
00375 }
00376
00377
00378
00379
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406 template<class ScalarType, class MV, class OP>
00407 int SVQBOrthoManager<ScalarType, MV, OP>::findBasis(
00408 MV &X, Teuchos::RCP<MV> MX,
00409 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00410 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00411 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00412 bool normalize_in) const {
00413
00414 const ScalarType ONE = SCT::one();
00415 const MagnitudeType MONE = SCTM::one();
00416 const MagnitudeType ZERO = SCTM::zero();
00417
00418 int numGS = 0,
00419 numSVQB = 0,
00420 numRand = 0;
00421
00422
00423 int xc = MVT::GetNumberVecs(X);
00424 int xr = MVT::GetVecLength( X );
00425
00426
00427 int nq = Q.length();
00428 int qr = (nq == 0) ? 0 : MVT::GetVecLength(*Q[0]);
00429 int qsize = 0;
00430 std::vector<int> qcs(nq);
00431 for (int i=0; i<nq; i++) {
00432 qcs[i] = MVT::GetNumberVecs(*Q[i]);
00433 qsize += qcs[i];
00434 }
00435
00436 if (normalize_in == true && qsize + xc > xr) {
00437
00438 TEST_FOR_EXCEPTION( true, std::invalid_argument,
00439 "Anasazi::SVQBOrthoManager::findBasis(): Orthogonalization constraints not feasible" );
00440 }
00441
00442
00443 if (normalize_in == false && (qsize == 0 || xc == 0)) {
00444
00445 return 0;
00446 }
00447 else if (normalize_in == true && (xc == 0 || xr == 0)) {
00448
00449 TEST_FOR_EXCEPTION( true, std::invalid_argument,
00450 "Anasazi::SVQBOrthoManager::findBasis(): X must be non-empty" );
00451 }
00452
00453
00454 TEST_FOR_EXCEPTION( qsize != 0 && qr != xr , std::invalid_argument,
00455 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with size of Q" );
00456
00457
00458
00459
00460
00461 C.resize(nq);
00462 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > newC(nq);
00463
00464 for (int i=0; i<nq; i++) {
00465
00466 TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument,
00467 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not mutually consistant" );
00468 TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument,
00469 "Anasazi::SVQBOrthoManager::findBasis(): Q has less rows than columns" );
00470
00471 if ( C[i] == Teuchos::null ) {
00472 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) );
00473 }
00474 else {
00475 TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc, std::invalid_argument,
00476 "Anasazi::SVQBOrthoManager::findBasis(): Size of Q not consistant with C" );
00477 }
00478
00479 C[i]->putScalar(ZERO);
00480 newC[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(*C[i]) );
00481 }
00482
00483
00485
00486
00487
00488
00489 if (normalize_in == true) {
00490 if ( B == Teuchos::null ) {
00491 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00492 }
00493 TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
00494 "Anasazi::SVQBOrthoManager::findBasis(): Size of B not consistant with X" );
00495
00496 B->putScalar(ZERO);
00497 for (int i=0; i<xc; i++) {
00498 (*B)(i,i) = MONE;
00499 }
00500 }
00501
00502
00503
00504
00505
00506
00507
00508
00509 Teuchos::RCP<MV> workX;
00510 if (normalize_in) {
00511 workX = MVT::Clone(X,xc);
00512 }
00513 if (this->_hasOp) {
00514 if (MX == Teuchos::null) {
00515
00516 MX = MVT::Clone(X,xc);
00517 OPT::Apply(*(this->_Op),X,*MX);
00518 this->_OpCounter += MVT::GetNumberVecs(X);
00519 }
00520 }
00521 else {
00522 MX = Teuchos::rcpFromRef(X);
00523 }
00524 std::vector<ScalarType> normX(xc), invnormX(xc);
00525 Teuchos::SerialDenseMatrix<int,ScalarType> XtMX(xc,xc), workU(1,1);
00526 Teuchos::LAPACK<int,ScalarType> lapack;
00527
00528
00529
00530
00531 std::vector<ScalarType> work;
00532 std::vector<MagnitudeType> lambda, lambdahi, rwork;
00533 if (normalize_in) {
00534
00535 int lwork = lapack.ILAENV(1,"hetrd","VU",xc,-1,-1,-1);
00536
00537
00538 TEST_FOR_EXCEPTION( lwork < 0, OrthoError,
00539 "Anasazi::SVQBOrthoManager::findBasis(): Error code from ILAENV" );
00540
00541 lwork = (lwork+2)*xc;
00542 work.resize(lwork);
00543
00544 lwork = (3*xc-2 > 1) ? 3*xc - 2 : 1;
00545 rwork.resize(lwork);
00546
00547 lambda.resize(xc);
00548 lambdahi.resize(xc);
00549 workU.reshape(xc,xc);
00550 }
00551
00552
00553 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00554 int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX ) : xr;
00555 TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
00556 "Anasazi::SVQBOrthoManager::findBasis(): Size of X not consistant with MX" );
00557
00558
00559 bool doGramSchmidt = true;
00560
00561 MagnitudeType tolerance = MONE/SCTM::squareroot(eps_);
00562
00563
00564 while (doGramSchmidt) {
00565
00567
00568 if (qsize > 0) {
00569
00570 numGS++;
00571
00572
00573 {
00574 std::vector<MagnitudeType> normX_mag(xc);
00575 normMat(X,normX_mag,MX);
00576 for (int i=0; i<xc; ++i) {
00577 normX[i] = normX_mag[i];
00578 invnormX[i] = (normX_mag[i] == ZERO) ? ZERO : MONE/normX_mag[i];
00579 }
00580 }
00581
00582 MVT::MvScale(X,invnormX);
00583 if (this->_hasOp) {
00584 MVT::MvScale(*MX,invnormX);
00585 }
00586
00587 if (debug_) {
00588 std::vector<MagnitudeType> nrm2(xc);
00589 std::cout << dbgstr << "max post-scale norm: (with/without MX) : ";
00590 MagnitudeType maxpsnw = ZERO, maxpsnwo = ZERO;
00591 normMat(X,nrm2,MX);
00592 for (int i=0; i<xc; i++) {
00593 maxpsnw = (nrm2[i] > maxpsnw ? nrm2[i] : maxpsnw);
00594 }
00595 this->norm(X,nrm2);
00596 for (int i=0; i<xc; i++) {
00597 maxpsnwo = (nrm2[i] > maxpsnwo ? nrm2[i] : maxpsnwo);
00598 }
00599 std::cout << "(" << maxpsnw << "," << maxpsnwo << ")" << std::endl;
00600 }
00601
00602 for (int i=0; i<nq; i++) {
00603 innerProdMat(*Q[i],X,*newC[i],Teuchos::null,MX);
00604 }
00605
00606 for (int i=0; i<nq; i++) {
00607 MVT::MvTimesMatAddMv(-ONE,*Q[i],*newC[i],ONE,X);
00608 }
00609
00610 MVT::MvScale(X,normX);
00611
00612 if (this->_hasOp) {
00613 OPT::Apply(*(this->_Op),X,*MX);
00614 this->_OpCounter += MVT::GetNumberVecs(X);
00615 }
00616
00617
00618
00619
00620
00621
00622 MagnitudeType maxNorm = ZERO;
00623 for (int j=0; j<xc; j++) {
00624 MagnitudeType sum = ZERO;
00625 for (int k=0; k<nq; k++) {
00626 for (int i=0; i<qcs[k]; i++) {
00627 sum += SCT::magnitude((*newC[k])(i,j))*SCT::magnitude((*newC[k])(i,j));
00628 }
00629 }
00630 maxNorm = (sum > maxNorm) ? sum : maxNorm;
00631 }
00632
00633
00634 if (maxNorm < 0.36) {
00635 doGramSchmidt = false;
00636 }
00637
00638
00639 for (int k=0; k<nq; k++) {
00640 for (int j=0; j<xc; j++) {
00641 for (int i=0; i<qcs[k]; i++) {
00642 (*newC[k])(i,j) *= normX[j];
00643 }
00644 }
00645 }
00646
00647 if (normalize_in) {
00648
00649 int info;
00650 for (int i=0; i<nq; i++) {
00651 info = C[i]->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,*newC[i],*B,ONE);
00652 TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
00653 }
00654 }
00655 else {
00656
00657 for (int i=0; i<nq; i++) {
00658 (*C[i]) += *newC[i];
00659 }
00660 }
00661 }
00662 else {
00663
00664 doGramSchmidt = false;
00665 }
00666
00667
00669
00670 if (normalize_in) {
00671
00672 MagnitudeType condT = tolerance;
00673
00674 while (condT >= tolerance) {
00675
00676 numSVQB++;
00677
00678
00679 innerProdMat(X,X,XtMX,MX,MX);
00680
00681
00682 std::vector<MagnitudeType> Dh(xc), Dhi(xc);
00683 for (int i=0; i<xc; i++) {
00684 Dh[i] = SCT::magnitude(SCT::squareroot(XtMX(i,i)));
00685 Dhi[i] = (Dh[i] == ZERO ? ZERO : MONE/Dh[i]);
00686 }
00687
00688 for (int i=0; i<xc; i++) {
00689 for (int j=0; j<xc; j++) {
00690 XtMX(i,j) *= Dhi[i]*Dhi[j];
00691 }
00692 }
00693
00694
00695 int info;
00696 lapack.HEEV('V', 'U', xc, XtMX.values(), XtMX.stride(), &lambda[0], &work[0], work.size(), &rwork[0], &info);
00697 std::stringstream os;
00698 os << "Anasazi::SVQBOrthoManager::findBasis(): Error code " << info << " from HEEV";
00699 TEST_FOR_EXCEPTION( info != 0, OrthoError,
00700 os.str() );
00701 if (debug_) {
00702 std::cout << dbgstr << "eigenvalues of XtMX: (";
00703 for (int i=0; i<xc-1; i++) {
00704 std::cout << lambda[i] << ",";
00705 }
00706 std::cout << lambda[xc-1] << ")" << std::endl;
00707 }
00708
00709
00710
00711 MagnitudeType maxLambda = lambda[xc-1],
00712 minLambda = lambda[0];
00713 int iZeroMax = -1;
00714 for (int i=0; i<xc; i++) {
00715 if (lambda[i] <= 10*eps_*maxLambda) {
00716 iZeroMax = i;
00717 lambda[i] = ZERO;
00718 lambdahi[i] = ZERO;
00719 }
00720
00721
00722
00723
00724
00725
00726 else {
00727 lambda[i] = SCTM::squareroot(lambda[i]);
00728 lambdahi[i] = MONE/lambda[i];
00729 }
00730 }
00731
00732
00733
00734
00735 std::vector<int> ind(xc);
00736 for (int i=0; i<xc; i++) {ind[i] = i;}
00737 MVT::SetBlock(X,ind,*workX);
00738
00739
00740 workU.assign(XtMX);
00741 for (int j=0; j<xc; j++) {
00742 for (int i=0; i<xc; i++) {
00743 workU(i,j) *= Dhi[i]*lambdahi[j];
00744 }
00745 }
00746
00747 MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,X);
00748
00749
00750
00751
00752 if (this->_hasOp) {
00753 if (maxLambda >= tolerance * minLambda) {
00754
00755 OPT::Apply(*(this->_Op),X,*MX);
00756 this->_OpCounter += MVT::GetNumberVecs(X);
00757 }
00758 else {
00759
00760
00761 MVT::SetBlock(*MX,ind,*workX);
00762
00763
00764 MVT::MvTimesMatAddMv(ONE,*workX,workU,ZERO,*MX);
00765 }
00766 }
00767
00768
00769
00770 for (int j=0; j<xc; j++) {
00771 for (int i=0; i<xc; i++) {
00772 workU(i,j) = Dh[i] * (*B)(i,j);
00773 }
00774 }
00775 info = B->multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,XtMX,workU,ZERO);
00776 TEST_FOR_EXCEPTION(info != 0, std::logic_error, "Anasazi::SVQBOrthoManager::findBasis(): Input error to SerialDenseMatrix::multiply.");
00777 for (int j=0; j<xc ;j++) {
00778 for (int i=0; i<xc; i++) {
00779 (*B)(i,j) *= lambda[i];
00780 }
00781 }
00782
00783
00784 if (iZeroMax >= 0) {
00785 if (debug_) {
00786 std::cout << dbgstr << "augmenting multivec with " << iZeroMax+1 << " random directions" << std::endl;
00787 }
00788
00789 numRand++;
00790
00791 std::vector<int> ind2(iZeroMax+1);
00792 for (int i=0; i<iZeroMax+1; i++) {
00793 ind2[i] = i;
00794 }
00795 Teuchos::RCP<MV> Xnull,MXnull;
00796 Xnull = MVT::CloneView(X,ind2);
00797 MVT::MvRandom(*Xnull);
00798 if (this->_hasOp) {
00799 MXnull = MVT::CloneView(*MX,ind2);
00800 OPT::Apply(*(this->_Op),*Xnull,*MXnull);
00801 this->_OpCounter += MVT::GetNumberVecs(*Xnull);
00802 MXnull = Teuchos::null;
00803 }
00804 Xnull = Teuchos::null;
00805 condT = tolerance;
00806 doGramSchmidt = true;
00807 break;
00808 }
00809
00810 condT = SCTM::magnitude(maxLambda / minLambda);
00811 if (debug_) {
00812 std::cout << dbgstr << "condT: " << condT << std::endl;
00813 }
00814
00815 }
00816
00817 if ((doGramSchmidt == false) && (condT > SCTM::squareroot(tolerance))) {
00818 doGramSchmidt = true;
00819 }
00820 }
00821
00822
00823 }
00824
00825 if (debug_) {
00826 std::cout << dbgstr << "(numGS,numSVQB,numRand) : "
00827 << "(" << numGS
00828 << "," << numSVQB
00829 << "," << numRand
00830 << ")" << std::endl;
00831 }
00832 return xc;
00833 }
00834
00835
00836 }
00837
00838 #endif // ANASAZI_SVQB_ORTHOMANAGER_HPP
00839