|
Anasazi
Version of the Day
|
|
Anasazi
Version of the Day
|
00001 // @HEADER 00002 // *********************************************************************** 00003 // 00004 // Anasazi: Block Eigensolvers Package 00005 // Copyright (2004) Sandia Corporation 00006 // 00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive 00008 // license for use of this work by or on behalf of the U.S. Government. 00009 // 00010 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 00025 // 00026 // *********************************************************************** 00027 // @HEADER 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 TEUCHOS_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 // ! Routine to find an orthonormal basis 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 // Constructor 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 TEUCHOS_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 TEUCHOS_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 // Compute the distance from orthonormality 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 MatOrthoManager<ScalarType,MV,OP>::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 // Compute the distance from orthogonality 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 MatOrthoManager<ScalarType,MV,OP>::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 // Find an Op-orthonormal basis for span(X), with rank numvectors(X) 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 // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X) 00475 // findBasis() requires MX 00476 00477 int xc = MVT::GetNumberVecs(X); 00478 int xr = MVT::GetVecLength(X); 00479 00480 // if Op==null, MX == X (via pointer) 00481 // Otherwise, either the user passed in MX or we will allocated and compute it 00482 if (this->_hasOp) { 00483 if (MX == Teuchos::null) { 00484 // we need to allocate space for MX 00485 MX = MVT::Clone(X,xc); 00486 OPT::Apply(*(this->_Op),X,*MX); 00487 this->_OpCounter += MVT::GetNumberVecs(X); 00488 } 00489 } 00490 00491 // if the user doesn't want to store the coefficients, 00492 // allocate some local memory for them 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 // check size of C, B 00501 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 00502 "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" ); 00503 TEUCHOS_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 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 00506 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" ); 00507 TEUCHOS_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 // Find an Op-orthonormal basis for span(X) - span(W) 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 // For the inner product defined by the operator Op or the identity (Op == 0) 00546 // -> Orthogonalize S against each Y[i], modifying it in the range of X[i] 00547 // Modify MS accordingly 00548 // 00549 // Note that when Op is 0, MS is not referenced 00550 // 00551 // Parameter variables 00552 // 00553 // S : Multivector to be transformed 00554 // 00555 // MS : Image of the block vector S by the mass matrix 00556 // 00557 // X,Y: Bases to orthogonalize against/via. 00558 // 00559 00560 #ifdef ANASAZI_ICGS_DEBUG 00561 // Get a FancyOStream from out_arg or create a new one ... 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 TEUCHOS_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 // short-circuit 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 // if we don't have enough C, expand it with null references 00587 // if we have too many, resize to throw away the latter ones 00588 // if we have exactly as many as we have X,Y this call has no effect 00589 // 00590 // same for MX, MY 00591 C.resize(numxy); 00592 MX.resize(numxy); 00593 MY.resize(numxy); 00594 00595 // check size of S w.r.t. common sense 00596 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument, 00597 "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." ); 00598 00599 // check size of MS 00600 if (this->_hasOp == true) { 00601 if (MS != Teuchos::null) { 00602 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*MS) != sr, std::invalid_argument, 00603 "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." ); 00604 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument, 00605 "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." ); 00606 } 00607 } 00608 00609 // tally up size of all X,Y and check/allocate C 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 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*X[i]) != sr, std::invalid_argument, 00616 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] length not consistent with S." ); 00617 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*Y[i]) != sr, std::invalid_argument, 00618 "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i << "] length not consistent with S." ); 00619 TEUCHOS_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 TEUCHOS_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 // check size of C[i] 00628 if ( C[i] == Teuchos::null ) { 00629 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) ); 00630 } 00631 else { 00632 TEUCHOS_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 // check sizes of MX[i], MY[i] if present 00636 // if not present, count their absence 00637 if (MX[i] != Teuchos::null) { 00638 TEUCHOS_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 TEUCHOS_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 // if one is null and the other is not... the user may have made a mistake 00654 TEUCHOS_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 // is this operation even feasible? 00662 TEUCHOS_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 /****** DO NO MODIFY *MS IF _hasOp == false 00668 * if _hasOp == false, we don't need MS, MX or MY 00669 * 00670 * if _hasOp == true, we need MS (for S M-norms) and 00671 * therefore, we must also update MS, regardless of whether 00672 * it gets returned to the user (i.e., whether the user provided it) 00673 * we may need to allocate and compute MX or MY 00674 * 00675 * let BXY denote: 00676 * "X" - we have all M*X[i] 00677 * "Y" - we have all M*Y[i] 00678 * "B" - we have biorthogonality (for all X[i],Y[i]) 00679 * Encode these as values 00680 * X = 1 00681 * Y = 2 00682 * B = 4 00683 * We have 8 possibilities, 0-7 00684 * 00685 * We must allocate storage and computer the following, lest 00686 * innerProdMat do it for us: 00687 * none (0) - allocate MX, for inv(<X,Y>) and M*S 00688 * 00689 * for the following, we will compute M*S manually before returning 00690 * B(4) BY(6) Y(2) --> updateMS = 1 00691 * for the following, we will update M*S as we go, using M*X 00692 * XY(3) X(1) none(0) BXY(7) BX(5) --> updateMS = 2 00693 * these choices favor applications of M over allocation of memory 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 // produce MS 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 // allocate the rest 00735 if (whichAlloc == 0) { 00736 // allocate and compute missing MX 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 // Op == I --> MS == S 00752 MS = Teuchos::rcpFromRef(S); 00753 updateMS = 0; 00754 } 00755 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error, 00756 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic."); 00757 00758 00760 // Perform the Gram-Schmidt transformation for a block of vectors 00762 00763 // Compute Cholesky factorizations for the Y'*M*X 00764 // YMX stores the YMX (initially) and their Cholesky factorizations (utlimately) 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 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Y[i],*X[i],*YMX[i],MY[i],MX[i]); 00773 #ifdef ANASAZI_ICGS_DEBUG 00774 // YMX should be symmetric positive definite 00775 // the cholesky will check the positive definiteness, but it looks only at the upper half 00776 // we will check the symmetry by removing the upper half from the lower half, which should 00777 // result in zeros 00778 // also, diagonal of YMX should be real; consider the complex part as error 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 // take the LU factorization 00791 int info; 00792 lapack.POTRF('U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info); 00793 TEUCHOS_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 // Compute the initial Op-norms 00799 #ifdef ANASAZI_ICGS_DEBUG 00800 std::vector<MagnitudeType> oldNorms(sc); 00801 MatOrthoManager<ScalarType,MV,OP>::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 // clear the C[i] and allocate Ccur 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 // Define the product Y[i]'*Op*S 00821 for (int i=0; i<numxy; i++) { 00822 // Compute Y[i]'*M*S 00823 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Y[i],S,Ccur[i],MY[i],MS); 00824 if (isBiortho == false) { 00825 // C[i] = inv(YMX[i])*(Y[i]'*M*S) 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 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, 00831 "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info << " from lapack::POTRS." ); 00832 } 00833 00834 // Multiply by X[i] and subtract the result in S 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 // Accumulate coeffs into previous step 00841 *C[i] += Ccur[i]; 00842 00843 // Update MS as required 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 // Compute new Op-norms 00860 #ifdef ANASAZI_ICGS_DEBUG 00861 std::vector<MagnitudeType> newNorms(sc); 00862 MatOrthoManager<ScalarType,MV,OP>::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 // Find an Op-orthonormal basis for span(S) - span(Y) 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 // For the inner product defined by the operator Op or the identity (Op == 0) 00891 // -> Orthogonalize S against each Y[i], modifying it in the range of X[i] 00892 // Modify MS accordingly 00893 // Then construct a M-orthonormal basis for the remaining part 00894 // 00895 // Note that when Op is 0, MS is not referenced 00896 // 00897 // Parameter variables 00898 // 00899 // S : Multivector to be transformed 00900 // 00901 // MS : Image of the block vector S by the mass matrix 00902 // 00903 // X,Y: Bases to orthogonalize against/via. 00904 // 00905 00906 #ifdef ANASAZI_ICGS_DEBUG 00907 // Get a FancyOStream from out_arg or create a new one ... 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 TEUCHOS_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 // short-circuit 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 // if we don't have enough C, expand it with null references 00929 // if we have too many, resize to throw away the latter ones 00930 // if we have exactly as many as we have X,Y this call has no effect 00931 // 00932 // same for MX, MY 00933 C.resize(numxy); 00934 MX.resize(numxy); 00935 MY.resize(numxy); 00936 00937 // check size of S w.r.t. common sense 00938 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument, 00939 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." ); 00940 00941 // check size of MS 00942 if (this->_hasOp == true) { 00943 if (MS != Teuchos::null) { 00944 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*MS) != sr, std::invalid_argument, 00945 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." ); 00946 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument, 00947 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." ); 00948 } 00949 } 00950 00951 // tally up size of all X,Y and check/allocate C 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 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*X[i]) != sr, std::invalid_argument, 00958 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] length not consistent with S." ); 00959 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*Y[i]) != sr, std::invalid_argument, 00960 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i << "] length not consistent with S." ); 00961 TEUCHOS_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 TEUCHOS_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 // check size of C[i] 00970 if ( C[i] == Teuchos::null ) { 00971 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) ); 00972 } 00973 else { 00974 TEUCHOS_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 // check sizes of MX[i], MY[i] if present 00978 // if not present, count their absence 00979 if (MX[i] != Teuchos::null) { 00980 TEUCHOS_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 TEUCHOS_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 // if one is null and the other is not... the user may have made a mistake 00996 TEUCHOS_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 // is this operation even feasible? 01004 TEUCHOS_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 /****** DO NO MODIFY *MS IF _hasOp == false 01011 * if _hasOp == false, we don't need MS, MX or MY 01012 * 01013 * if _hasOp == true, we need MS (for S M-norms and normalization) and 01014 * therefore, we must also update MS, regardless of whether 01015 * it gets returned to the user (i.e., whether the user provided it) 01016 * we may need to allocate and compute MX or MY 01017 * 01018 * let BXY denote: 01019 * "X" - we have all M*X[i] 01020 * "Y" - we have all M*Y[i] 01021 * "B" - we have biorthogonality (for all X[i],Y[i]) 01022 * Encode these as values 01023 * X = 1 01024 * Y = 2 01025 * B = 4 01026 * We have 8 possibilities, 0-7 01027 * 01028 * We must allocate storage and computer the following, lest 01029 * innerProdMat do it for us: 01030 * none (0) - allocate MX, for inv(<X,Y>) and M*S 01031 * 01032 * for the following, we will compute M*S manually before returning 01033 * B(4) BY(6) Y(2) --> updateMS = 1 01034 * for the following, we will update M*S as we go, using M*X 01035 * XY(3) X(1) none(0) BXY(7) BX(5) --> updateMS = 2 01036 * these choices favor applications of M over allocation of memory 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 // produce MS 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 // allocate the rest 01078 if (whichAlloc == 0) { 01079 // allocate and compute missing MX 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 // Op == I --> MS == S 01095 MS = Teuchos::rcpFromRef(S); 01096 updateMS = 0; 01097 } 01098 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error, 01099 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic."); 01100 01101 01102 // if the user doesn't want to store the coefficients, 01103 // allocate some local memory for them 01104 if ( B == Teuchos::null ) { 01105 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(sc,sc) ); 01106 } 01107 else { 01108 // check size of B 01109 TEUCHOS_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 // orthogonalize all of S against X,Y 01115 projectGen(S,X,Y,isBiortho,C,MS,MX,MY); 01116 01117 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(sc,1); 01118 // start working 01119 rank = 0; 01120 int numTries = 10; // each vector in X gets 10 random chances to escape degeneracy 01121 int oldrank = -1; 01122 do { 01123 int curssize = sc - rank; 01124 01125 // orthonormalize S, but quit if it is rank deficient 01126 // we can't let findBasis generated random vectors to complete the basis, 01127 // because it doesn't know about Q; we will do this ourselves below 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 // we had previously stopped before, after operating on vector oldrank 01135 // we saved its coefficients, augmented it with a random vector, and 01136 // then called findBasis() again, which proceeded to add vector oldrank 01137 // to the basis. 01138 // now, restore the saved coefficients into B 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 // we quit on this vector and will augment it with random below 01147 // this is the first time that we have quit on this vector 01148 // therefor, (*B)(:,rank) contains the actual coefficients of the 01149 // input vectors with respect to the previous vectors in the basis 01150 // save these values, as (*B)(:,rank) will be overwritten by our next 01151 // call to findBasis() 01152 // we will restore it after we are done working on this vector 01153 for (int i=0; i<sc; i++) { 01154 oldCoeff(i,0) = (*B)(i,rank); 01155 } 01156 } 01157 } 01158 01159 if (rank == sc) { 01160 // we are done 01161 #ifdef ANASAZI_ICGS_DEBUG 01162 *out << "Finished computing basis.\n"; 01163 #endif 01164 break; 01165 } 01166 else { 01167 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError, 01168 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen"); 01169 01170 if (rank != oldrank) { 01171 // we added a basis vector from random info; reset the chance counter 01172 numTries = 10; 01173 // store old rank 01174 oldrank = rank; 01175 } 01176 else { 01177 // has this vector run out of chances to escape degeneracy? 01178 if (numTries <= 0) { 01179 break; 01180 } 01181 } 01182 // use one of this vector's chances 01183 numTries--; 01184 01185 // randomize troubled direction 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::CloneViewNonConst(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::CloneViewNonConst(*MS,ind); 01199 OPT::Apply( *(this->_Op), *curS, *curMS ); 01200 this->_OpCounter += MVT::GetNumberVecs(*curS); 01201 } 01202 01203 // orthogonalize against X,Y 01204 // if !this->_hasOp, the curMS will be ignored. 01205 // we don't care about these coefficients 01206 // on the contrary, we need to preserve the previous coeffs 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 // this should never raise an exception; but our post-conditions oblige us to check 01215 TEUCHOS_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 // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 01229 // the rank is numvectors(X) 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 // For the inner product defined by the operator Op or the identity (Op == 0) 01237 // -> Orthonormalize X 01238 // Modify MX accordingly 01239 // 01240 // Note that when Op is 0, MX is not referenced 01241 // 01242 // Parameter variables 01243 // 01244 // X : Vectors to be orthonormalized 01245 // 01246 // MX : Image of the multivector X under the operator Op 01247 // 01248 // Op : Pointer to the operator for the inner product 01249 // 01250 01251 #ifdef ANASAZI_ICGS_DEBUG 01252 // Get a FancyOStream from out_arg or create a new one ... 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 * If _hasOp == false, we will not reference MX below * 01270 *******************************************************/ 01271 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error, 01272 "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS."); 01273 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error, 01274 "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" ); 01275 01276 /* xstart is which column we are starting the process with, based on howMany 01277 * columns before xstart are assumed to be Op-orthonormal already 01278 */ 01279 int xstart = xc - howMany; 01280 01281 for (int j = xstart; j < xc; j++) { 01282 01283 // numX represents the number of currently orthonormal columns of X 01284 int numX = j; 01285 // j represents the index of the current column of X 01286 // these are different interpretations of the same value 01287 01288 // 01289 // set the lower triangular part of B to zero 01290 for (int i=j+1; i<xc; ++i) { 01291 B(i,j) = ZERO; 01292 } 01293 01294 // Get a view of the vector currently being worked on. 01295 std::vector<int> index(1); 01296 index[0] = j; 01297 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index ); 01298 Teuchos::RCP<MV> MXj; 01299 if ((this->_hasOp)) { 01300 // MXj is a view of the current vector in MX 01301 MXj = MVT::CloneViewNonConst( *MX, index ); 01302 } 01303 else { 01304 // MXj is a pointer to Xj, and MUST NOT be modified 01305 MXj = Xj; 01306 } 01307 01308 // Get a view of the previous vectors. 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 /* numTrials>0 will denote that the current vector was randomized for the purpose 01322 * of finding a basis vector, and that the coefficients of that vector should 01323 * not be stored in B 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 // Make storage for these Gram-Schmidt iterations. 01331 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1); 01332 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1); 01333 01334 // 01335 // Save old MXj vector and compute Op-norm 01336 // 01337 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 01338 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,origNorm,MXj); 01339 #ifdef ANASAZI_ICGS_DEBUG 01340 *out << "origNorm = " << origNorm[0] << "\n"; 01341 #endif 01342 01343 if (numX > 0) { 01344 // Apply the first step of Gram-Schmidt 01345 01346 // product <- prevX^T MXj 01347 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj); 01348 01349 // Xj <- Xj - prevX prevX^T MXj 01350 // = Xj - prevX product 01351 #ifdef ANASAZI_ICGS_DEBUG 01352 *out << "Orthogonalizing X[" << j << "]...\n"; 01353 #endif 01354 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj ); 01355 01356 // Update MXj 01357 if (this->_hasOp) { 01358 // MXj <- Op*Xj_new 01359 // = Op*(Xj_old - prevX prevX^T MXj) 01360 // = MXj - prevMX product 01361 #ifdef ANASAZI_ICGS_DEBUG 01362 *out << "Updating MX[" << j << "]...\n"; 01363 #endif 01364 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj ); 01365 } 01366 01367 // Compute new Op-norm 01368 MatOrthoManager<ScalarType,MV,OP>::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 // Check if a correction is needed. 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 // Apply the second step of Gram-Schmidt 01387 // This is the same as above 01388 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1); 01389 MatOrthoManager<ScalarType,MV,OP>::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 // Compute new Op-norms 01402 MatOrthoManager<ScalarType,MV,OP>::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 // we don't do another GS, we just set it to zero. 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 } // if (numX > 0) do GS 01431 01432 // save the coefficients, if we are working on the original vector and not a randomly generated one 01433 if (numTrials == 0) { 01434 for (int i=0; i<numX; i++) { 01435 B(i,j) = product(i,0); 01436 } 01437 } 01438 01439 // Check if Xj has any directional information left after the orthogonalization. 01440 MatOrthoManager<ScalarType,MV,OP>::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 // Normalize Xj. 01446 // Xj <- Xj / norm 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 // Update MXj. 01453 MVT::MvScale( *MXj, ONE/newNorm[0]); 01454 } 01455 01456 // save it, if it corresponds to the original vector and not a randomly generated one 01457 if (numTrials == 0) { 01458 B(j,j) = newNorm[0]; 01459 } 01460 01461 // We are not rank deficient in this vector. Move on to the next vector in X. 01462 rankDef = false; 01463 break; 01464 } 01465 else { 01466 #ifdef ANASAZI_ICGS_DEBUG 01467 *out << "Not normalizing M*X[" << j << "]...\n"; 01468 #endif 01469 // There was nothing left in Xj after orthogonalizing against previous columns in X. 01470 // X is rank deficient. 01471 // reflect this in the coefficients 01472 B(j,j) = ZERO; 01473 01474 if (completeBasis) { 01475 // Fill it with random information and keep going. 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 } // for (numTrials = 0; numTrials < 10; ++numTrials) 01494 01495 // if rankDef == true, then quit and notify user of rank obtained 01496 if (rankDef == true) { 01497 TEUCHOS_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 } // for (j = 0; j < xc; ++j) 01506 01507 #ifdef ANASAZI_ICGS_DEBUG 01508 *out << "Returning " << xc << " from Anasazi::ICGSOrthoManager::findBasis(...)\n"; 01509 #endif 01510 return xc; 01511 } 01512 01513 } // namespace Anasazi 01514 01515 #endif // ANASAZI_ICSG_ORTHOMANAGER_HPP 01516
1.7.6.1
