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_ICSG_ORTHOMANAGER_HPP
00035 #define ANASAZI_ICSG_ORTHOMANAGER_HPP
00036
00044 #include "AnasaziConfigDefs.hpp"
00045 #include "AnasaziMultiVecTraits.hpp"
00046 #include "AnasaziOperatorTraits.hpp"
00047 #include "AnasaziGenOrthoManager.hpp"
00048 #include "Teuchos_TimeMonitor.hpp"
00049 #include "Teuchos_LAPACK.hpp"
00050 #include "Teuchos_BLAS.hpp"
00051 #ifdef ANASAZI_ICGS_DEBUG
00052 # include <Teuchos_FancyOStream.hpp>
00053 #endif
00054
00055 namespace Anasazi {
00056
00057 template<class ScalarType, class MV, class OP>
00058 class ICGSOrthoManager : public GenOrthoManager<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 ICGSOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, int numIters = 2,
00072 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0,
00073 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 );
00074
00075
00077 ~ICGSOrthoManager() {}
00079
00080
00082
00083
00165 void projectGen(
00166 MV &S,
00167 Teuchos::Array<Teuchos::RCP<const MV> > X,
00168 Teuchos::Array<Teuchos::RCP<const MV> > Y,
00169 bool isBiOrtho,
00170 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00171 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00172 Teuchos::RCP<MV> MS = Teuchos::null,
00173 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
00174 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00175 ) const;
00176
00177
00269 int projectAndNormalizeGen (
00270 MV &S,
00271 Teuchos::Array<Teuchos::RCP<const MV> > X,
00272 Teuchos::Array<Teuchos::RCP<const MV> > Y,
00273 bool isBiOrtho,
00274 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00275 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00276 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00277 Teuchos::RCP<MV> MS = Teuchos::null,
00278 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)),
00279 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00280 ) const;
00281
00282
00284
00285
00287
00288
00289
00301 void projectMat (
00302 MV &X,
00303 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00304 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00305 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00306 Teuchos::RCP<MV> MX = Teuchos::null,
00307 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00308 ) const;
00309
00310
00319 int normalizeMat (
00320 MV &X,
00321 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00322 Teuchos::RCP<MV> MX = Teuchos::null
00323 ) const;
00324
00325
00334 int projectAndNormalizeMat (
00335 MV &X,
00336 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00337 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C
00338 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)),
00339 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null,
00340 Teuchos::RCP<MV> MX = Teuchos::null,
00341 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null))
00342 ) const;
00343
00345
00347
00348
00353 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00354 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const;
00355
00360 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00361 orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const;
00362
00364
00366
00367
00369 void setNumIters( int numIters ) {
00370 numIters_ = numIters;
00371 TEST_FOR_EXCEPTION(numIters_ < 1,std::invalid_argument,
00372 "Anasazi::ICGSOrthoManager::setNumIters(): input must be >= 1.");
00373 }
00374
00376 int getNumIters() const { return numIters_; }
00377
00379
00380 private:
00381 MagnitudeType eps_;
00382 MagnitudeType tol_;
00383
00385 int numIters_;
00386
00387
00388 int findBasis(MV &X, Teuchos::RCP<MV> MX,
00389 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
00390 bool completeBasis, int howMany = -1) const;
00391 };
00392
00393
00394
00396
00397 template<class ScalarType, class MV, class OP>
00398 ICGSOrthoManager<ScalarType,MV,OP>::ICGSOrthoManager( Teuchos::RCP<const OP> Op,
00399 int numIters,
00400 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps,
00401 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol) :
00402 GenOrthoManager<ScalarType,MV,OP>(Op), eps_(eps), tol_(tol)
00403 {
00404 setNumIters(numIters);
00405 TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument,
00406 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative.");
00407 if (eps_ == 0) {
00408 Teuchos::LAPACK<int,MagnitudeType> lapack;
00409 eps_ = lapack.LAMCH('E');
00410 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.50);
00411 }
00412 TEST_FOR_EXCEPTION(
00413 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()),
00414 std::invalid_argument,
00415 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1].");
00416 }
00417
00418
00419
00421
00422 template<class ScalarType, class MV, class OP>
00423 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00424 ICGSOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const {
00425 const ScalarType ONE = SCT::one();
00426 int rank = MVT::GetNumberVecs(X);
00427 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank);
00428 innerProdMat(X,X,xTx,MX,MX);
00429 for (int i=0; i<rank; i++) {
00430 xTx(i,i) -= ONE;
00431 }
00432 return xTx.normFrobenius();
00433 }
00434
00435
00436
00438
00439 template<class ScalarType, class MV, class OP>
00440 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType
00441 ICGSOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const {
00442 int r1 = MVT::GetNumberVecs(X1);
00443 int r2 = MVT::GetNumberVecs(X2);
00444 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2);
00445 innerProdMat(X1,X2,xTx,MX1,MX2);
00446 return xTx.normFrobenius();
00447 }
00448
00449
00450
00452 template<class ScalarType, class MV, class OP>
00453 void ICGSOrthoManager<ScalarType, MV, OP>::projectMat(
00454 MV &X,
00455 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00456 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00457 Teuchos::RCP<MV> MX,
00458 Teuchos::Array<Teuchos::RCP<const MV> > MQ
00459 ) const
00460 {
00461 projectGen(X,Q,Q,true,C,MX,MQ,MQ);
00462 }
00463
00464
00465
00467
00468 template<class ScalarType, class MV, class OP>
00469 int ICGSOrthoManager<ScalarType, MV, OP>::normalizeMat(
00470 MV &X,
00471 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00472 Teuchos::RCP<MV> MX) const
00473 {
00474
00475
00476
00477 int xc = MVT::GetNumberVecs(X);
00478 int xr = MVT::GetVecLength(X);
00479
00480
00481
00482 if (this->_hasOp) {
00483 if (MX == Teuchos::null) {
00484
00485 MX = MVT::Clone(X,xc);
00486 OPT::Apply(*(this->_Op),X,*MX);
00487 this->_OpCounter += MVT::GetNumberVecs(X);
00488 }
00489 }
00490
00491
00492
00493 if ( B == Teuchos::null ) {
00494 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) );
00495 }
00496
00497 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc;
00498 int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX ) : xr;
00499
00500
00501 TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument,
00502 "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" );
00503 TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument,
00504 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of B" );
00505 TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument,
00506 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" );
00507 TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument,
00508 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not feasible for normalization" );
00509
00510 return findBasis(X, MX, *B, true );
00511 }
00512
00513
00514
00516
00517 template<class ScalarType, class MV, class OP>
00518 int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat(
00519 MV &X,
00520 Teuchos::Array<Teuchos::RCP<const MV> > Q,
00521 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00522 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00523 Teuchos::RCP<MV> MX,
00524 Teuchos::Array<Teuchos::RCP<const MV> > MQ
00525 ) const
00526 {
00527 return projectAndNormalizeGen(X,Q,Q,true,C,B,MX,MQ,MQ);
00528 }
00529
00530
00531
00533 template<class ScalarType, class MV, class OP>
00534 void ICGSOrthoManager<ScalarType, MV, OP>::projectGen(
00535 MV &S,
00536 Teuchos::Array<Teuchos::RCP<const MV> > X,
00537 Teuchos::Array<Teuchos::RCP<const MV> > Y,
00538 bool isBiortho,
00539 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00540 Teuchos::RCP<MV> MS,
00541 Teuchos::Array<Teuchos::RCP<const MV> > MX,
00542 Teuchos::Array<Teuchos::RCP<const MV> > MY
00543 ) const
00544 {
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560 #ifdef ANASAZI_ICGS_DEBUG
00561
00562 Teuchos::RCP<Teuchos::FancyOStream>
00563 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00564 out->setShowAllFrontMatter(false).setShowProcRank(true);
00565 *out << "Entering Anasazi::ICGSOrthoManager::projectGen(...)\n";
00566 #endif
00567
00568 const ScalarType ONE = SCT::one();
00569 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
00570 Teuchos::LAPACK<int,ScalarType> lapack;
00571 Teuchos::BLAS<int,ScalarType> blas;
00572
00573 int sc = MVT::GetNumberVecs( S );
00574 int sr = MVT::GetVecLength( S );
00575 int numxy = X.length();
00576 TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument,
00577 "Anasazi::ICGSOrthoManager::projectGen(): X and Y must contain the same number of multivectors.");
00578 std::vector<int> xyc(numxy);
00579
00580 if (numxy == 0 || sc == 0 || sr == 0) {
00581 #ifdef ANASAZI_ICGS_DEBUG
00582 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
00583 #endif
00584 return;
00585 }
00586
00587
00588
00589
00590
00591 C.resize(numxy);
00592 MX.resize(numxy);
00593 MY.resize(numxy);
00594
00595
00596 TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument,
00597 "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." );
00598
00599
00600 if (this->_hasOp == true) {
00601 if (MS != Teuchos::null) {
00602 TEST_FOR_EXCEPTION( MVT::GetVecLength(*MS) != sr, std::invalid_argument,
00603 "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." );
00604 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
00605 "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." );
00606 }
00607 }
00608
00609
00610 int sumxyc = 0;
00611 int MYmissing = 0;
00612 int MXmissing = 0;
00613 for (int i=0; i<numxy; i++) {
00614 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
00615 TEST_FOR_EXCEPTION( MVT::GetVecLength(*X[i]) != sr, std::invalid_argument,
00616 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] length not consistent with S." );
00617 TEST_FOR_EXCEPTION( MVT::GetVecLength(*Y[i]) != sr, std::invalid_argument,
00618 "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i << "] length not consistent with S." );
00619 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
00620 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] and Y[" << i << "] widths not consistent." );
00621
00622 xyc[i] = MVT::GetNumberVecs( *X[i] );
00623 TEST_FOR_EXCEPTION( sr < xyc[i], std::invalid_argument,
00624 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "],Y[" << i << "] have less rows than columns, and therefore cannot be full rank." );
00625 sumxyc += xyc[i];
00626
00627
00628 if ( C[i] == Teuchos::null ) {
00629 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
00630 }
00631 else {
00632 TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument,
00633 "Anasazi::ICGSOrthoManager::projectGen(): Size of Q not consistent with size of C." );
00634 }
00635
00636
00637 if (MX[i] != Teuchos::null) {
00638 TEST_FOR_EXCEPTION( MVT::GetVecLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
00639 "Anasazi::ICGSOrthoManager::projectGen(): Size of MX[" << i << "] not consistent with size of X[" << i << "]." );
00640 }
00641 else {
00642 MXmissing += xyc[i];
00643 }
00644 if (MY[i] != Teuchos::null) {
00645 TEST_FOR_EXCEPTION( MVT::GetVecLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
00646 "Anasazi::ICGSOrthoManager::projectGen(): Size of MY[" << i << "] not consistent with size of Y[" << i << "]." );
00647 }
00648 else {
00649 MYmissing += xyc[i];
00650 }
00651 }
00652 else {
00653
00654 TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
00655 "Anasazi::ICGSOrthoManager::projectGen(): "
00656 << (X[i] == Teuchos::null ? "Y[" : "X[") << i << "] was provided but "
00657 << (X[i] == Teuchos::null ? "X[" : "Y[") << i << "] was not.");
00658 }
00659 }
00660
00661
00662 TEST_FOR_EXCEPTION(sumxyc > sr, std::invalid_argument,
00663 "Anasazi::ICGSOrthoManager::projectGen(): dimension of all X[i],Y[i] is "
00664 << sumxyc << ", but length of vectors is only " << sr << ". This is infeasible.");
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 int updateMS = -1;
00697 if (this->_hasOp) {
00698 int whichAlloc = 0;
00699 if (MXmissing == 0) {
00700 whichAlloc += 1;
00701 }
00702 if (MYmissing == 0) {
00703 whichAlloc += 2;
00704 }
00705 if (isBiortho) {
00706 whichAlloc += 4;
00707 }
00708
00709 switch (whichAlloc) {
00710 case 2:
00711 case 4:
00712 case 6:
00713 updateMS = 1;
00714 break;
00715 case 0:
00716 case 1:
00717 case 3:
00718 case 5:
00719 case 7:
00720 updateMS = 2;
00721 break;
00722 }
00723
00724
00725 if (MS == Teuchos::null) {
00726 #ifdef ANASAZI_ICGS_DEBUG
00727 *out << "Allocating MS...\n";
00728 #endif
00729 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
00730 OPT::Apply(*(this->_Op),S,*MS);
00731 this->_OpCounter += MVT::GetNumberVecs(S);
00732 }
00733
00734
00735 if (whichAlloc == 0) {
00736
00737 for (int i=0; i<numxy; i++) {
00738 if (MX[i] == Teuchos::null) {
00739 #ifdef ANASAZI_ICGS_DEBUG
00740 *out << "Allocating MX[" << i << "]...\n";
00741 #endif
00742 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
00743 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
00744 MX[i] = tmpMX;
00745 this->_OpCounter += xyc[i];
00746 }
00747 }
00748 }
00749 }
00750 else {
00751
00752 MS = Teuchos::rcpFromRef(S);
00753 updateMS = 0;
00754 }
00755 TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
00756 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
00757
00758
00760
00762
00763
00764
00765 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > YMX(numxy);
00766 if (isBiortho == false) {
00767 for (int i=0; i<numxy; i++) {
00768 #ifdef ANASAZI_ICGS_DEBUG
00769 *out << "Computing YMX[" << i << "] and its Cholesky factorization...\n";
00770 #endif
00771 YMX[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],xyc[i]) );
00772 innerProdMat(*Y[i],*X[i],*YMX[i],MY[i],MX[i]);
00773 #ifdef ANASAZI_ICGS_DEBUG
00774
00775
00776
00777
00778
00779 {
00780 MagnitudeType err = ZERO;
00781 for (int jj=0; jj<YMX[i]->numCols(); jj++) {
00782 err =+ SCT::magnitude(SCT::imag((*YMX[i])(jj,jj)));
00783 for (int ii=jj; ii<YMX[i]->numRows(); ii++) {
00784 err += SCT::magnitude( (*YMX[i])(ii,jj) - SCT::conjugate((*YMX[i])(jj,ii)) );
00785 }
00786 }
00787 *out << "Symmetry error in YMX[" << i << "] == " << err << "\n";
00788 }
00789 #endif
00790
00791 int info;
00792 lapack.POTRF('U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info);
00793 TEST_FOR_EXCEPTION(info != 0,std::logic_error,
00794 "Anasazi::ICGSOrthoManager::projectGen(): Error computing Cholesky factorization of Y[i]^T * M * X[i] using POTRF: returned info " << info);
00795 }
00796 }
00797
00798
00799 #ifdef ANASAZI_ICGS_DEBUG
00800 std::vector<MagnitudeType> oldNorms(sc);
00801 normMat(S,oldNorms,MS);
00802 *out << "oldNorms = { ";
00803 std::copy(oldNorms.begin(), oldNorms.end(), std::ostream_iterator<MagnitudeType>(*out, " "));
00804 *out << "}\n";
00805 #endif
00806
00807
00808
00809 Teuchos::Array<Teuchos::SerialDenseMatrix<int,ScalarType> > Ccur(numxy);
00810 for (int i=0; i<numxy; i++) {
00811 C[i]->putScalar(ZERO);
00812 Ccur[i].reshape(C[i]->numRows(),C[i]->numCols());
00813 }
00814
00815 for (int iter=0; iter < numIters_; iter++) {
00816 #ifdef ANASAZI_ICGS_DEBUG
00817 *out << "beginning iteration " << iter+1 << "\n";
00818 #endif
00819
00820
00821 for (int i=0; i<numxy; i++) {
00822
00823 innerProdMat(*Y[i],S,Ccur[i],MY[i],MS);
00824 if (isBiortho == false) {
00825
00826 int info;
00827 lapack.POTRS('U',YMX[i]->numCols(),Ccur[i].numCols(),
00828 YMX[i]->values(),YMX[i]->stride(),
00829 Ccur[i].values(),Ccur[i].stride(), &info);
00830 TEST_FOR_EXCEPTION(info != 0, std::logic_error,
00831 "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info << " from lapack::POTRS." );
00832 }
00833
00834
00835 #ifdef ANASAZI_ICGS_DEBUG
00836 *out << "Applying projector P_{X[" << i << "],Y[" << i << "]}...\n";
00837 #endif
00838 MVT::MvTimesMatAddMv( -ONE, *X[i], Ccur[i], ONE, S );
00839
00840
00841 *C[i] += Ccur[i];
00842
00843
00844 if (updateMS == 1) {
00845 #ifdef ANASAZI_ICGS_DEBUG
00846 *out << "Updating MS...\n";
00847 #endif
00848 OPT::Apply( *(this->_Op), S, *MS);
00849 this->_OpCounter += sc;
00850 }
00851 else if (updateMS == 2) {
00852 #ifdef ANASAZI_ICGS_DEBUG
00853 *out << "Updating MS...\n";
00854 #endif
00855 MVT::MvTimesMatAddMv( -ONE, *MX[i], Ccur[i], ONE, *MS );
00856 }
00857 }
00858
00859
00860 #ifdef ANASAZI_ICGS_DEBUG
00861 std::vector<MagnitudeType> newNorms(sc);
00862 normMat(S,newNorms,MS);
00863 *out << "newNorms = { ";
00864 std::copy(newNorms.begin(), newNorms.end(), std::ostream_iterator<MagnitudeType>(*out, " "));
00865 *out << "}\n";
00866 #endif
00867 }
00868
00869 #ifdef ANASAZI_ICGS_DEBUG
00870 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
00871 #endif
00872 }
00873
00874
00875
00877
00878 template<class ScalarType, class MV, class OP>
00879 int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalizeGen(
00880 MV &S,
00881 Teuchos::Array<Teuchos::RCP<const MV> > X,
00882 Teuchos::Array<Teuchos::RCP<const MV> > Y,
00883 bool isBiortho,
00884 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C,
00885 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B,
00886 Teuchos::RCP<MV> MS,
00887 Teuchos::Array<Teuchos::RCP<const MV> > MX,
00888 Teuchos::Array<Teuchos::RCP<const MV> > MY
00889 ) const {
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906 #ifdef ANASAZI_ICGS_DEBUG
00907
00908 Teuchos::RCP<Teuchos::FancyOStream>
00909 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
00910 out->setShowAllFrontMatter(false).setShowProcRank(true);
00911 *out << "Entering Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
00912 #endif
00913
00914 int rank;
00915 int sc = MVT::GetNumberVecs( S );
00916 int sr = MVT::GetVecLength( S );
00917 int numxy = X.length();
00918 TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument,
00919 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X and Y must contain the same number of multivectors.");
00920 std::vector<int> xyc(numxy);
00921
00922 if (sc == 0 || sr == 0) {
00923 #ifdef ANASAZI_ICGS_DEBUG
00924 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n";
00925 #endif
00926 return 0;
00927 }
00928
00929
00930
00931
00932
00933 C.resize(numxy);
00934 MX.resize(numxy);
00935 MY.resize(numxy);
00936
00937
00938 TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument,
00939 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." );
00940
00941
00942 if (this->_hasOp == true) {
00943 if (MS != Teuchos::null) {
00944 TEST_FOR_EXCEPTION( MVT::GetVecLength(*MS) != sr, std::invalid_argument,
00945 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." );
00946 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument,
00947 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." );
00948 }
00949 }
00950
00951
00952 int sumxyc = 0;
00953 int MYmissing = 0;
00954 int MXmissing = 0;
00955 for (int i=0; i<numxy; i++) {
00956 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) {
00957 TEST_FOR_EXCEPTION( MVT::GetVecLength(*X[i]) != sr, std::invalid_argument,
00958 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] length not consistent with S." );
00959 TEST_FOR_EXCEPTION( MVT::GetVecLength(*Y[i]) != sr, std::invalid_argument,
00960 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i << "] length not consistent with S." );
00961 TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument,
00962 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] and Y[" << i << "] widths not consistent." );
00963
00964 xyc[i] = MVT::GetNumberVecs( *X[i] );
00965 TEST_FOR_EXCEPTION( sr < xyc[i], std::invalid_argument,
00966 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "],Y[" << i << "] have less rows than columns, and therefore cannot be full rank." );
00967 sumxyc += xyc[i];
00968
00969
00970 if ( C[i] == Teuchos::null ) {
00971 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) );
00972 }
00973 else {
00974 TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument,
00975 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of Q not consistent with size of C." );
00976 }
00977
00978
00979 if (MX[i] != Teuchos::null) {
00980 TEST_FOR_EXCEPTION( MVT::GetVecLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument,
00981 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MX[" << i << "] not consistent with size of X[" << i << "]." );
00982 }
00983 else {
00984 MXmissing += xyc[i];
00985 }
00986 if (MY[i] != Teuchos::null) {
00987 TEST_FOR_EXCEPTION( MVT::GetVecLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument,
00988 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MY[" << i << "] not consistent with size of Y[" << i << "]." );
00989 }
00990 else {
00991 MYmissing += xyc[i];
00992 }
00993 }
00994 else {
00995
00996 TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument,
00997 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): "
00998 << (X[i] == Teuchos::null ? "Y[" : "X[") << i << "] was provided but "
00999 << (X[i] == Teuchos::null ? "X[" : "Y[") << i << "] was not.");
01000 }
01001 }
01002
01003
01004 TEST_FOR_EXCEPTION(sumxyc + sc > sr, std::invalid_argument,
01005 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): dimension of all X[i],Y[i] is "
01006 << sumxyc << " and requested " << sc << "-dimensional basis, but length of vectors is only "
01007 << sr << ". This is infeasible.");
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020
01021
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039 int updateMS = -1;
01040 if (this->_hasOp) {
01041 int whichAlloc = 0;
01042 if (MXmissing == 0) {
01043 whichAlloc += 1;
01044 }
01045 if (MYmissing == 0) {
01046 whichAlloc += 2;
01047 }
01048 if (isBiortho) {
01049 whichAlloc += 4;
01050 }
01051
01052 switch (whichAlloc) {
01053 case 2:
01054 case 4:
01055 case 6:
01056 updateMS = 1;
01057 break;
01058 case 0:
01059 case 1:
01060 case 3:
01061 case 5:
01062 case 7:
01063 updateMS = 2;
01064 break;
01065 }
01066
01067
01068 if (MS == Teuchos::null) {
01069 #ifdef ANASAZI_ICGS_DEBUG
01070 *out << "Allocating MS...\n";
01071 #endif
01072 MS = MVT::Clone(S,MVT::GetNumberVecs(S));
01073 OPT::Apply(*(this->_Op),S,*MS);
01074 this->_OpCounter += MVT::GetNumberVecs(S);
01075 }
01076
01077
01078 if (whichAlloc == 0) {
01079
01080 for (int i=0; i<numxy; i++) {
01081 if (MX[i] == Teuchos::null) {
01082 #ifdef ANASAZI_ICGS_DEBUG
01083 *out << "Allocating MX[" << i << "]...\n";
01084 #endif
01085 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]);
01086 OPT::Apply(*(this->_Op),*X[i],*tmpMX);
01087 MX[i] = tmpMX;
01088 this->_OpCounter += xyc[i];
01089 }
01090 }
01091 }
01092 }
01093 else {
01094
01095 MS = Teuchos::rcpFromRef(S);
01096 updateMS = 0;
01097 }
01098 TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error,
01099 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic.");
01100
01101
01102
01103
01104 if ( B == Teuchos::null ) {
01105 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(sc,sc) );
01106 }
01107 else {
01108
01109 TEST_FOR_EXCEPTION( B->numRows() != sc || B->numCols() != sc, std::invalid_argument,
01110 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of S must be consistent with size of B" );
01111 }
01112
01113
01114
01115 projectGen(S,X,Y,isBiortho,C,MS,MX,MY);
01116
01117 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(sc,1);
01118
01119 rank = 0;
01120 int numTries = 10;
01121 int oldrank = -1;
01122 do {
01123 int curssize = sc - rank;
01124
01125
01126
01127
01128 #ifdef ANASAZI_ICGS_DEBUG
01129 *out << "Attempting to find orthonormal basis for X...\n";
01130 #endif
01131 rank = findBasis(S,MS,*B,false,curssize);
01132
01133 if (oldrank != -1 && rank != oldrank) {
01134
01135
01136
01137
01138
01139 for (int i=0; i<sc; i++) {
01140 (*B)(i,oldrank) = oldCoeff(i,0);
01141 }
01142 }
01143
01144 if (rank < sc) {
01145 if (rank != oldrank) {
01146
01147
01148
01149
01150
01151
01152
01153 for (int i=0; i<sc; i++) {
01154 oldCoeff(i,0) = (*B)(i,rank);
01155 }
01156 }
01157 }
01158
01159 if (rank == sc) {
01160
01161 #ifdef ANASAZI_ICGS_DEBUG
01162 *out << "Finished computing basis.\n";
01163 #endif
01164 break;
01165 }
01166 else {
01167 TEST_FOR_EXCEPTION( rank < oldrank, OrthoError,
01168 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen");
01169
01170 if (rank != oldrank) {
01171
01172 numTries = 10;
01173
01174 oldrank = rank;
01175 }
01176 else {
01177
01178 if (numTries <= 0) {
01179 break;
01180 }
01181 }
01182
01183 numTries--;
01184
01185
01186 #ifdef ANASAZI_ICGS_DEBUG
01187 *out << "Inserting random vector in X[" << rank << "]. Attempt " << 10-numTries << ".\n";
01188 #endif
01189 Teuchos::RCP<MV> curS, curMS;
01190 std::vector<int> ind(1);
01191 ind[0] = rank;
01192 curS = MVT::CloneView(S,ind);
01193 MVT::MvRandom(*curS);
01194 if (this->_hasOp) {
01195 #ifdef ANASAZI_ICGS_DEBUG
01196 *out << "Applying operator to random vector.\n";
01197 #endif
01198 curMS = MVT::CloneView(*MS,ind);
01199 OPT::Apply( *(this->_Op), *curS, *curMS );
01200 this->_OpCounter += MVT::GetNumberVecs(*curS);
01201 }
01202
01203
01204
01205
01206
01207 {
01208 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0);
01209 projectGen(*curS,X,Y,isBiortho,dummyC,curMS,MX,MY);
01210 }
01211 }
01212 } while (1);
01213
01214
01215 TEST_FOR_EXCEPTION( rank > sc || rank < 0, std::logic_error,
01216 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Debug error in rank variable." );
01217
01218 #ifdef ANASAZI_ICGS_DEBUG
01219 *out << "Returning " << rank << " from Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n";
01220 #endif
01221
01222 return rank;
01223 }
01224
01225
01226
01228
01229
01230 template<class ScalarType, class MV, class OP>
01231 int ICGSOrthoManager<ScalarType, MV, OP>::findBasis(
01232 MV &X, Teuchos::RCP<MV> MX,
01233 Teuchos::SerialDenseMatrix<int,ScalarType> &B,
01234 bool completeBasis, int howMany ) const {
01235
01236
01237
01238
01239
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251 #ifdef ANASAZI_ICGS_DEBUG
01252
01253 Teuchos::RCP<Teuchos::FancyOStream>
01254 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
01255 out->setShowAllFrontMatter(false).setShowProcRank(true);
01256 *out << "Entering Anasazi::ICGSOrthoManager::findBasis(...)\n";
01257 #endif
01258
01259 const ScalarType ONE = SCT::one();
01260 const MagnitudeType ZERO = SCT::magnitude(SCT::zero());
01261
01262 int xc = MVT::GetNumberVecs( X );
01263
01264 if (howMany == -1) {
01265 howMany = xc;
01266 }
01267
01268
01269
01270
01271 TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error,
01272 "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS.");
01273 TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error,
01274 "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" );
01275
01276
01277
01278
01279 int xstart = xc - howMany;
01280
01281 for (int j = xstart; j < xc; j++) {
01282
01283
01284 int numX = j;
01285
01286
01287
01288
01289
01290 for (int i=j+1; i<xc; ++i) {
01291 B(i,j) = ZERO;
01292 }
01293
01294
01295 std::vector<int> index(1);
01296 index[0] = j;
01297 Teuchos::RCP<MV> Xj = MVT::CloneView( X, index );
01298 Teuchos::RCP<MV> MXj;
01299 if ((this->_hasOp)) {
01300
01301 MXj = MVT::CloneView( *MX, index );
01302 }
01303 else {
01304
01305 MXj = Xj;
01306 }
01307
01308
01309 std::vector<int> prev_idx( numX );
01310 Teuchos::RCP<const MV> prevX, prevMX;
01311
01312 if (numX > 0) {
01313 for (int i=0; i<numX; ++i) prev_idx[i] = i;
01314 prevX = MVT::CloneView( X, prev_idx );
01315 if (this->_hasOp) {
01316 prevMX = MVT::CloneView( *MX, prev_idx );
01317 }
01318 }
01319
01320 bool rankDef = true;
01321
01322
01323
01324
01325 for (int numTrials = 0; numTrials < 10; numTrials++) {
01326 #ifdef ANASAZI_ICGS_DEBUG
01327 *out << "Trial " << numTrials << " for vector " << j << "\n";
01328 #endif
01329
01330
01331 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1);
01332 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1);
01333
01334
01335
01336
01337 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj );
01338 normMat(*Xj,origNorm,MXj);
01339 #ifdef ANASAZI_ICGS_DEBUG
01340 *out << "origNorm = " << origNorm[0] << "\n";
01341 #endif
01342
01343 if (numX > 0) {
01344
01345
01346
01347 innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj);
01348
01349
01350
01351 #ifdef ANASAZI_ICGS_DEBUG
01352 *out << "Orthogonalizing X[" << j << "]...\n";
01353 #endif
01354 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj );
01355
01356
01357 if (this->_hasOp) {
01358
01359
01360
01361 #ifdef ANASAZI_ICGS_DEBUG
01362 *out << "Updating MX[" << j << "]...\n";
01363 #endif
01364 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj );
01365 }
01366
01367
01368 normMat(*Xj,newNorm,MXj);
01369 MagnitudeType product_norm = product.normOne();
01370
01371 #ifdef ANASAZI_ICGS_DEBUG
01372 *out << "newNorm = " << newNorm[0] << "\n";
01373 *out << "prodoct_norm = " << product_norm << "\n";
01374 #endif
01375
01376
01377 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) {
01378 #ifdef ANASAZI_ICGS_DEBUG
01379 if (product_norm/newNorm[0] >= tol_) {
01380 *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n";
01381 }
01382 else {
01383 *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n";
01384 }
01385 #endif
01386
01387
01388 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1);
01389 innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj);
01390 product += P2;
01391 #ifdef ANASAZI_ICGS_DEBUG
01392 *out << "Orthogonalizing X[" << j << "]...\n";
01393 #endif
01394 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj );
01395 if ((this->_hasOp)) {
01396 #ifdef ANASAZI_ICGS_DEBUG
01397 *out << "Updating MX[" << j << "]...\n";
01398 #endif
01399 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj );
01400 }
01401
01402 normMat(*Xj,newNorm2,MXj);
01403 product_norm = P2.normOne();
01404 #ifdef ANASAZI_ICGS_DEBUG
01405 *out << "newNorm2 = " << newNorm2[0] << "\n";
01406 *out << "product_norm = " << product_norm << "\n";
01407 #endif
01408 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) {
01409
01410 #ifdef ANASAZI_ICGS_DEBUG
01411 if (product_norm/newNorm2[0] >= tol_) {
01412 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n";
01413 }
01414 else if (newNorm[0] < newNorm2[0]) {
01415 *out << "newNorm2 > newNorm... setting vector to zero.\n";
01416 }
01417 else {
01418 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n";
01419 }
01420 #endif
01421 MVT::MvInit(*Xj,ZERO);
01422 if ((this->_hasOp)) {
01423 #ifdef ANASAZI_ICGS_DEBUG
01424 *out << "Setting MX[" << j << "] to zero as well.\n";
01425 #endif
01426 MVT::MvInit(*MXj,ZERO);
01427 }
01428 }
01429 }
01430 }
01431
01432
01433 if (numTrials == 0) {
01434 for (int i=0; i<numX; i++) {
01435 B(i,j) = product(i,0);
01436 }
01437 }
01438
01439
01440 normMat(*Xj,newNorm,MXj);
01441 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) {
01442 #ifdef ANASAZI_ICGS_DEBUG
01443 *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n";
01444 #endif
01445
01446
01447 MVT::MvScale( *Xj, ONE/newNorm[0]);
01448 if (this->_hasOp) {
01449 #ifdef ANASAZI_ICGS_DEBUG
01450 *out << "Normalizing M*X[" << j << "]...\n";
01451 #endif
01452
01453 MVT::MvScale( *MXj, ONE/newNorm[0]);
01454 }
01455
01456
01457 if (numTrials == 0) {
01458 B(j,j) = newNorm[0];
01459 }
01460
01461
01462 rankDef = false;
01463 break;
01464 }
01465 else {
01466 #ifdef ANASAZI_ICGS_DEBUG
01467 *out << "Not normalizing M*X[" << j << "]...\n";
01468 #endif
01469
01470
01471
01472 B(j,j) = ZERO;
01473
01474 if (completeBasis) {
01475
01476 #ifdef ANASAZI_ICGS_DEBUG
01477 *out << "Inserting random vector in X[" << j << "]...\n";
01478 #endif
01479 MVT::MvRandom( *Xj );
01480 if (this->_hasOp) {
01481 #ifdef ANASAZI_ICGS_DEBUG
01482 *out << "Updating M*X[" << j << "]...\n";
01483 #endif
01484 OPT::Apply( *(this->_Op), *Xj, *MXj );
01485 this->_OpCounter += MVT::GetNumberVecs(*Xj);
01486 }
01487 }
01488 else {
01489 rankDef = true;
01490 break;
01491 }
01492 }
01493 }
01494
01495
01496 if (rankDef == true) {
01497 TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError,
01498 "Anasazi::ICGSOrthoManager::findBasis(): Unable to complete basis" );
01499 #ifdef ANASAZI_ICGS_DEBUG
01500 *out << "Returning early, rank " << j << " from Anasazi::ICGSOrthoManager::findBasis(...)\n";
01501 #endif
01502 return j;
01503 }
01504
01505 }
01506
01507 #ifdef ANASAZI_ICGS_DEBUG
01508 *out << "Returning " << xc << " from Anasazi::ICGSOrthoManager::findBasis(...)\n";
01509 #endif
01510 return xc;
01511 }
01512
01513 }
01514
01515 #endif // ANASAZI_ICSG_ORTHOMANAGER_HPP
01516