|
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., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 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 MultiVecTraitsExt<ScalarType,MV> MVText; 00065 typedef OperatorTraits<ScalarType,MV,OP> OPT; 00066 00067 public: 00068 00070 00071 00072 ICGSOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, int numIters = 2, 00073 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0, 00074 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 ); 00075 00076 00078 ~ICGSOrthoManager() {} 00080 00081 00083 00084 00166 void projectGen( 00167 MV &S, 00168 Teuchos::Array<Teuchos::RCP<const MV> > X, 00169 Teuchos::Array<Teuchos::RCP<const MV> > Y, 00170 bool isBiOrtho, 00171 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C 00172 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)), 00173 Teuchos::RCP<MV> MS = Teuchos::null, 00174 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)), 00175 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)) 00176 ) const; 00177 00178 00270 int projectAndNormalizeGen ( 00271 MV &S, 00272 Teuchos::Array<Teuchos::RCP<const MV> > X, 00273 Teuchos::Array<Teuchos::RCP<const MV> > Y, 00274 bool isBiOrtho, 00275 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C 00276 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)), 00277 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null, 00278 Teuchos::RCP<MV> MS = Teuchos::null, 00279 Teuchos::Array<Teuchos::RCP<const MV> > MX = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)), 00280 Teuchos::Array<Teuchos::RCP<const MV> > MY = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)) 00281 ) const; 00282 00283 00285 00286 00288 00289 00290 00302 void projectMat ( 00303 MV &X, 00304 Teuchos::Array<Teuchos::RCP<const MV> > Q, 00305 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C 00306 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)), 00307 Teuchos::RCP<MV> MX = Teuchos::null, 00308 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)) 00309 ) const; 00310 00311 00320 int normalizeMat ( 00321 MV &X, 00322 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null, 00323 Teuchos::RCP<MV> MX = Teuchos::null 00324 ) const; 00325 00326 00335 int projectAndNormalizeMat ( 00336 MV &X, 00337 Teuchos::Array<Teuchos::RCP<const MV> > Q, 00338 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C 00339 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)), 00340 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null, 00341 Teuchos::RCP<MV> MX = Teuchos::null, 00342 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)) 00343 ) const; 00344 00346 00348 00349 00354 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 00355 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const; 00356 00361 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 00362 orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const; 00363 00365 00367 00368 00370 void setNumIters( int numIters ) { 00371 numIters_ = numIters; 00372 TEUCHOS_TEST_FOR_EXCEPTION(numIters_ < 1,std::invalid_argument, 00373 "Anasazi::ICGSOrthoManager::setNumIters(): input must be >= 1."); 00374 } 00375 00377 int getNumIters() const { return numIters_; } 00378 00380 00381 private: 00382 MagnitudeType eps_; 00383 MagnitudeType tol_; 00384 00386 int numIters_; 00387 00388 // ! Routine to find an orthonormal basis 00389 int findBasis(MV &X, Teuchos::RCP<MV> MX, 00390 Teuchos::SerialDenseMatrix<int,ScalarType> &B, 00391 bool completeBasis, int howMany = -1) const; 00392 }; 00393 00394 00395 00397 // Constructor 00398 template<class ScalarType, class MV, class OP> 00399 ICGSOrthoManager<ScalarType,MV,OP>::ICGSOrthoManager( Teuchos::RCP<const OP> Op, 00400 int numIters, 00401 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps, 00402 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol) : 00403 GenOrthoManager<ScalarType,MV,OP>(Op), eps_(eps), tol_(tol) 00404 { 00405 setNumIters(numIters); 00406 TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument, 00407 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"eps\" must be non-negative."); 00408 if (eps_ == 0) { 00409 Teuchos::LAPACK<int,MagnitudeType> lapack; 00410 eps_ = lapack.LAMCH('E'); 00411 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.50); 00412 } 00413 TEUCHOS_TEST_FOR_EXCEPTION( 00414 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()), 00415 std::invalid_argument, 00416 "Anasazi::ICGSOrthoManager::ICGSOrthoManager(): argument \"tol\" must be in [0,1]."); 00417 } 00418 00419 00420 00422 // Compute the distance from orthonormality 00423 template<class ScalarType, class MV, class OP> 00424 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 00425 ICGSOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const { 00426 const ScalarType ONE = SCT::one(); 00427 int rank = MVT::GetNumberVecs(X); 00428 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank); 00429 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X,X,xTx,MX,MX); 00430 for (int i=0; i<rank; i++) { 00431 xTx(i,i) -= ONE; 00432 } 00433 return xTx.normFrobenius(); 00434 } 00435 00436 00437 00439 // Compute the distance from orthogonality 00440 template<class ScalarType, class MV, class OP> 00441 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 00442 ICGSOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const { 00443 int r1 = MVT::GetNumberVecs(X1); 00444 int r2 = MVT::GetNumberVecs(X2); 00445 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2); 00446 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X1,X2,xTx,MX1,MX2); 00447 return xTx.normFrobenius(); 00448 } 00449 00450 00451 00453 template<class ScalarType, class MV, class OP> 00454 void ICGSOrthoManager<ScalarType, MV, OP>::projectMat( 00455 MV &X, 00456 Teuchos::Array<Teuchos::RCP<const MV> > Q, 00457 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 00458 Teuchos::RCP<MV> MX, 00459 Teuchos::Array<Teuchos::RCP<const MV> > MQ 00460 ) const 00461 { 00462 projectGen(X,Q,Q,true,C,MX,MQ,MQ); 00463 } 00464 00465 00466 00468 // Find an Op-orthonormal basis for span(X), with rank numvectors(X) 00469 template<class ScalarType, class MV, class OP> 00470 int ICGSOrthoManager<ScalarType, MV, OP>::normalizeMat( 00471 MV &X, 00472 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 00473 Teuchos::RCP<MV> MX) const 00474 { 00475 // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X) 00476 // findBasis() requires MX 00477 00478 int xc = MVT::GetNumberVecs(X); 00479 ptrdiff_t xr = MVText::GetGlobalLength(X); 00480 00481 // if Op==null, MX == X (via pointer) 00482 // Otherwise, either the user passed in MX or we will allocated and compute it 00483 if (this->_hasOp) { 00484 if (MX == Teuchos::null) { 00485 // we need to allocate space for MX 00486 MX = MVT::Clone(X,xc); 00487 OPT::Apply(*(this->_Op),X,*MX); 00488 this->_OpCounter += MVT::GetNumberVecs(X); 00489 } 00490 } 00491 00492 // if the user doesn't want to store the coefficients, 00493 // allocate some local memory for them 00494 if ( B == Teuchos::null ) { 00495 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) ); 00496 } 00497 00498 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc; 00499 ptrdiff_t mxr = (this->_hasOp) ? MVText::GetGlobalLength( *MX ) : xr; 00500 00501 // check size of C, B 00502 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 00503 "Anasazi::ICGSOrthoManager::normalizeMat(): X must be non-empty" ); 00504 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 00505 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of B" ); 00506 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 00507 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not consistent with size of MX" ); 00508 TEUCHOS_TEST_FOR_EXCEPTION( static_cast<ptrdiff_t>(xc) > xr, std::invalid_argument, 00509 "Anasazi::ICGSOrthoManager::normalizeMat(): Size of X not feasible for normalization" ); 00510 00511 return findBasis(X, MX, *B, true ); 00512 } 00513 00514 00515 00517 // Find an Op-orthonormal basis for span(X) - span(W) 00518 template<class ScalarType, class MV, class OP> 00519 int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat( 00520 MV &X, 00521 Teuchos::Array<Teuchos::RCP<const MV> > Q, 00522 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 00523 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 00524 Teuchos::RCP<MV> MX, 00525 Teuchos::Array<Teuchos::RCP<const MV> > MQ 00526 ) const 00527 { 00528 return projectAndNormalizeGen(X,Q,Q,true,C,B,MX,MQ,MQ); 00529 } 00530 00531 00532 00534 template<class ScalarType, class MV, class OP> 00535 void ICGSOrthoManager<ScalarType, MV, OP>::projectGen( 00536 MV &S, 00537 Teuchos::Array<Teuchos::RCP<const MV> > X, 00538 Teuchos::Array<Teuchos::RCP<const MV> > Y, 00539 bool isBiortho, 00540 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 00541 Teuchos::RCP<MV> MS, 00542 Teuchos::Array<Teuchos::RCP<const MV> > MX, 00543 Teuchos::Array<Teuchos::RCP<const MV> > MY 00544 ) const 00545 { 00546 // For the inner product defined by the operator Op or the identity (Op == 0) 00547 // -> Orthogonalize S against each Y[i], modifying it in the range of X[i] 00548 // Modify MS accordingly 00549 // 00550 // Note that when Op is 0, MS is not referenced 00551 // 00552 // Parameter variables 00553 // 00554 // S : Multivector to be transformed 00555 // 00556 // MS : Image of the block vector S by the mass matrix 00557 // 00558 // X,Y: Bases to orthogonalize against/via. 00559 // 00560 00561 #ifdef ANASAZI_ICGS_DEBUG 00562 // Get a FancyOStream from out_arg or create a new one ... 00563 Teuchos::RCP<Teuchos::FancyOStream> 00564 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout)); 00565 out->setShowAllFrontMatter(false).setShowProcRank(true); 00566 *out << "Entering Anasazi::ICGSOrthoManager::projectGen(...)\n"; 00567 #endif 00568 00569 const ScalarType ONE = SCT::one(); 00570 const MagnitudeType ZERO = SCT::magnitude(SCT::zero()); 00571 Teuchos::LAPACK<int,ScalarType> lapack; 00572 Teuchos::BLAS<int,ScalarType> blas; 00573 00574 int sc = MVT::GetNumberVecs( S ); 00575 ptrdiff_t sr = MVText::GetGlobalLength( S ); 00576 int numxy = X.length(); 00577 TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument, 00578 "Anasazi::ICGSOrthoManager::projectGen(): X and Y must contain the same number of multivectors."); 00579 std::vector<int> xyc(numxy); 00580 // short-circuit 00581 if (numxy == 0 || sc == 0 || sr == 0) { 00582 #ifdef ANASAZI_ICGS_DEBUG 00583 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n"; 00584 #endif 00585 return; 00586 } 00587 // if we don't have enough C, expand it with null references 00588 // if we have too many, resize to throw away the latter ones 00589 // if we have exactly as many as we have X,Y this call has no effect 00590 // 00591 // same for MX, MY 00592 C.resize(numxy); 00593 MX.resize(numxy); 00594 MY.resize(numxy); 00595 00596 // check size of S w.r.t. common sense 00597 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument, 00598 "Anasazi::ICGSOrthoManager::projectGen(): MVT returned negative dimensions for S." ); 00599 00600 // check size of MS 00601 if (this->_hasOp == true) { 00602 if (MS != Teuchos::null) { 00603 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MS) != sr, std::invalid_argument, 00604 "Anasazi::ICGSOrthoManager::projectGen(): MS length not consistent with S." ); 00605 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument, 00606 "Anasazi::ICGSOrthoManager::projectGen(): MS width not consistent with S." ); 00607 } 00608 } 00609 00610 // tally up size of all X,Y and check/allocate C 00611 ptrdiff_t sumxyc = 0; 00612 int MYmissing = 0; 00613 int MXmissing = 0; 00614 for (int i=0; i<numxy; i++) { 00615 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) { 00616 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*X[i]) != sr, std::invalid_argument, 00617 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] length not consistent with S." ); 00618 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*Y[i]) != sr, std::invalid_argument, 00619 "Anasazi::ICGSOrthoManager::projectGen(): Y[" << i << "] length not consistent with S." ); 00620 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument, 00621 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "] and Y[" << i << "] widths not consistent." ); 00622 00623 xyc[i] = MVT::GetNumberVecs( *X[i] ); 00624 TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument, 00625 "Anasazi::ICGSOrthoManager::projectGen(): X[" << i << "],Y[" << i << "] have less rows than columns, and therefore cannot be full rank." ); 00626 sumxyc += xyc[i]; 00627 00628 // check size of C[i] 00629 if ( C[i] == Teuchos::null ) { 00630 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) ); 00631 } 00632 else { 00633 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument, 00634 "Anasazi::ICGSOrthoManager::projectGen(): Size of Q not consistent with size of C." ); 00635 } 00636 // check sizes of MX[i], MY[i] if present 00637 // if not present, count their absence 00638 if (MX[i] != Teuchos::null) { 00639 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument, 00640 "Anasazi::ICGSOrthoManager::projectGen(): Size of MX[" << i << "] not consistent with size of X[" << i << "]." ); 00641 } 00642 else { 00643 MXmissing += xyc[i]; 00644 } 00645 if (MY[i] != Teuchos::null) { 00646 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument, 00647 "Anasazi::ICGSOrthoManager::projectGen(): Size of MY[" << i << "] not consistent with size of Y[" << i << "]." ); 00648 } 00649 else { 00650 MYmissing += xyc[i]; 00651 } 00652 } 00653 else { 00654 // if one is null and the other is not... the user may have made a mistake 00655 TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument, 00656 "Anasazi::ICGSOrthoManager::projectGen(): " 00657 << (X[i] == Teuchos::null ? "Y[" : "X[") << i << "] was provided but " 00658 << (X[i] == Teuchos::null ? "X[" : "Y[") << i << "] was not."); 00659 } 00660 } 00661 00662 // is this operation even feasible? 00663 TEUCHOS_TEST_FOR_EXCEPTION(sumxyc > sr, std::invalid_argument, 00664 "Anasazi::ICGSOrthoManager::projectGen(): dimension of all X[i],Y[i] is " 00665 << sumxyc << ", but length of vectors is only " << sr << ". This is infeasible."); 00666 00667 00668 /****** DO NO MODIFY *MS IF _hasOp == false 00669 * if _hasOp == false, we don't need MS, MX or MY 00670 * 00671 * if _hasOp == true, we need MS (for S M-norms) and 00672 * therefore, we must also update MS, regardless of whether 00673 * it gets returned to the user (i.e., whether the user provided it) 00674 * we may need to allocate and compute MX or MY 00675 * 00676 * let BXY denote: 00677 * "X" - we have all M*X[i] 00678 * "Y" - we have all M*Y[i] 00679 * "B" - we have biorthogonality (for all X[i],Y[i]) 00680 * Encode these as values 00681 * X = 1 00682 * Y = 2 00683 * B = 4 00684 * We have 8 possibilities, 0-7 00685 * 00686 * We must allocate storage and computer the following, lest 00687 * innerProdMat do it for us: 00688 * none (0) - allocate MX, for inv(<X,Y>) and M*S 00689 * 00690 * for the following, we will compute M*S manually before returning 00691 * B(4) BY(6) Y(2) --> updateMS = 1 00692 * for the following, we will update M*S as we go, using M*X 00693 * XY(3) X(1) none(0) BXY(7) BX(5) --> updateMS = 2 00694 * these choices favor applications of M over allocation of memory 00695 * 00696 */ 00697 int updateMS = -1; 00698 if (this->_hasOp) { 00699 int whichAlloc = 0; 00700 if (MXmissing == 0) { 00701 whichAlloc += 1; 00702 } 00703 if (MYmissing == 0) { 00704 whichAlloc += 2; 00705 } 00706 if (isBiortho) { 00707 whichAlloc += 4; 00708 } 00709 00710 switch (whichAlloc) { 00711 case 2: 00712 case 4: 00713 case 6: 00714 updateMS = 1; 00715 break; 00716 case 0: 00717 case 1: 00718 case 3: 00719 case 5: 00720 case 7: 00721 updateMS = 2; 00722 break; 00723 } 00724 00725 // produce MS 00726 if (MS == Teuchos::null) { 00727 #ifdef ANASAZI_ICGS_DEBUG 00728 *out << "Allocating MS...\n"; 00729 #endif 00730 MS = MVT::Clone(S,MVT::GetNumberVecs(S)); 00731 OPT::Apply(*(this->_Op),S,*MS); 00732 this->_OpCounter += MVT::GetNumberVecs(S); 00733 } 00734 00735 // allocate the rest 00736 if (whichAlloc == 0) { 00737 // allocate and compute missing MX 00738 for (int i=0; i<numxy; i++) { 00739 if (MX[i] == Teuchos::null) { 00740 #ifdef ANASAZI_ICGS_DEBUG 00741 *out << "Allocating MX[" << i << "]...\n"; 00742 #endif 00743 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]); 00744 OPT::Apply(*(this->_Op),*X[i],*tmpMX); 00745 MX[i] = tmpMX; 00746 this->_OpCounter += xyc[i]; 00747 } 00748 } 00749 } 00750 } 00751 else { 00752 // Op == I --> MS == S 00753 MS = Teuchos::rcpFromRef(S); 00754 updateMS = 0; 00755 } 00756 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error, 00757 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic."); 00758 00759 00761 // Perform the Gram-Schmidt transformation for a block of vectors 00763 00764 // Compute Cholesky factorizations for the Y'*M*X 00765 // YMX stores the YMX (initially) and their Cholesky factorizations (utlimately) 00766 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > YMX(numxy); 00767 if (isBiortho == false) { 00768 for (int i=0; i<numxy; i++) { 00769 #ifdef ANASAZI_ICGS_DEBUG 00770 *out << "Computing YMX[" << i << "] and its Cholesky factorization...\n"; 00771 #endif 00772 YMX[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],xyc[i]) ); 00773 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Y[i],*X[i],*YMX[i],MY[i],MX[i]); 00774 #ifdef ANASAZI_ICGS_DEBUG 00775 // YMX should be symmetric positive definite 00776 // the cholesky will check the positive definiteness, but it looks only at the upper half 00777 // we will check the symmetry by removing the upper half from the lower half, which should 00778 // result in zeros 00779 // also, diagonal of YMX should be real; consider the complex part as error 00780 { 00781 MagnitudeType err = ZERO; 00782 for (int jj=0; jj<YMX[i]->numCols(); jj++) { 00783 err =+ SCT::magnitude(SCT::imag((*YMX[i])(jj,jj))); 00784 for (int ii=jj; ii<YMX[i]->numRows(); ii++) { 00785 err += SCT::magnitude( (*YMX[i])(ii,jj) - SCT::conjugate((*YMX[i])(jj,ii)) ); 00786 } 00787 } 00788 *out << "Symmetry error in YMX[" << i << "] == " << err << "\n"; 00789 } 00790 #endif 00791 // take the LU factorization 00792 int info; 00793 lapack.POTRF('U',YMX[i]->numRows(),YMX[i]->values(),YMX[i]->stride(),&info); 00794 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error, 00795 "Anasazi::ICGSOrthoManager::projectGen(): Error computing Cholesky factorization of Y[i]^T * M * X[i] using POTRF: returned info " << info); 00796 } 00797 } 00798 00799 // Compute the initial Op-norms 00800 #ifdef ANASAZI_ICGS_DEBUG 00801 std::vector<MagnitudeType> oldNorms(sc); 00802 MatOrthoManager<ScalarType,MV,OP>::normMat(S,oldNorms,MS); 00803 *out << "oldNorms = { "; 00804 std::copy(oldNorms.begin(), oldNorms.end(), std::ostream_iterator<MagnitudeType>(*out, " ")); 00805 *out << "}\n"; 00806 #endif 00807 00808 00809 // clear the C[i] and allocate Ccur 00810 Teuchos::Array<Teuchos::SerialDenseMatrix<int,ScalarType> > Ccur(numxy); 00811 for (int i=0; i<numxy; i++) { 00812 C[i]->putScalar(ZERO); 00813 Ccur[i].reshape(C[i]->numRows(),C[i]->numCols()); 00814 } 00815 00816 for (int iter=0; iter < numIters_; iter++) { 00817 #ifdef ANASAZI_ICGS_DEBUG 00818 *out << "beginning iteration " << iter+1 << "\n"; 00819 #endif 00820 00821 // Define the product Y[i]'*Op*S 00822 for (int i=0; i<numxy; i++) { 00823 // Compute Y[i]'*M*S 00824 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Y[i],S,Ccur[i],MY[i],MS); 00825 if (isBiortho == false) { 00826 // C[i] = inv(YMX[i])*(Y[i]'*M*S) 00827 int info; 00828 lapack.POTRS('U',YMX[i]->numCols(),Ccur[i].numCols(), 00829 YMX[i]->values(),YMX[i]->stride(), 00830 Ccur[i].values(),Ccur[i].stride(), &info); 00831 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, std::logic_error, 00832 "Anasazi::ICGSOrthoManager::projectGen(): Error code " << info << " from lapack::POTRS." ); 00833 } 00834 00835 // Multiply by X[i] and subtract the result in S 00836 #ifdef ANASAZI_ICGS_DEBUG 00837 *out << "Applying projector P_{X[" << i << "],Y[" << i << "]}...\n"; 00838 #endif 00839 MVT::MvTimesMatAddMv( -ONE, *X[i], Ccur[i], ONE, S ); 00840 00841 // Accumulate coeffs into previous step 00842 *C[i] += Ccur[i]; 00843 00844 // Update MS as required 00845 if (updateMS == 1) { 00846 #ifdef ANASAZI_ICGS_DEBUG 00847 *out << "Updating MS...\n"; 00848 #endif 00849 OPT::Apply( *(this->_Op), S, *MS); 00850 this->_OpCounter += sc; 00851 } 00852 else if (updateMS == 2) { 00853 #ifdef ANASAZI_ICGS_DEBUG 00854 *out << "Updating MS...\n"; 00855 #endif 00856 MVT::MvTimesMatAddMv( -ONE, *MX[i], Ccur[i], ONE, *MS ); 00857 } 00858 } 00859 00860 // Compute new Op-norms 00861 #ifdef ANASAZI_ICGS_DEBUG 00862 std::vector<MagnitudeType> newNorms(sc); 00863 MatOrthoManager<ScalarType,MV,OP>::normMat(S,newNorms,MS); 00864 *out << "newNorms = { "; 00865 std::copy(newNorms.begin(), newNorms.end(), std::ostream_iterator<MagnitudeType>(*out, " ")); 00866 *out << "}\n"; 00867 #endif 00868 } 00869 00870 #ifdef ANASAZI_ICGS_DEBUG 00871 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n"; 00872 #endif 00873 } 00874 00875 00876 00878 // Find an Op-orthonormal basis for span(S) - span(Y) 00879 template<class ScalarType, class MV, class OP> 00880 int ICGSOrthoManager<ScalarType, MV, OP>::projectAndNormalizeGen( 00881 MV &S, 00882 Teuchos::Array<Teuchos::RCP<const MV> > X, 00883 Teuchos::Array<Teuchos::RCP<const MV> > Y, 00884 bool isBiortho, 00885 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 00886 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 00887 Teuchos::RCP<MV> MS, 00888 Teuchos::Array<Teuchos::RCP<const MV> > MX, 00889 Teuchos::Array<Teuchos::RCP<const MV> > MY 00890 ) const { 00891 // For the inner product defined by the operator Op or the identity (Op == 0) 00892 // -> Orthogonalize S against each Y[i], modifying it in the range of X[i] 00893 // Modify MS accordingly 00894 // Then construct a M-orthonormal basis for the remaining part 00895 // 00896 // Note that when Op is 0, MS is not referenced 00897 // 00898 // Parameter variables 00899 // 00900 // S : Multivector to be transformed 00901 // 00902 // MS : Image of the block vector S by the mass matrix 00903 // 00904 // X,Y: Bases to orthogonalize against/via. 00905 // 00906 00907 #ifdef ANASAZI_ICGS_DEBUG 00908 // Get a FancyOStream from out_arg or create a new one ... 00909 Teuchos::RCP<Teuchos::FancyOStream> 00910 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout)); 00911 out->setShowAllFrontMatter(false).setShowProcRank(true); 00912 *out << "Entering Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n"; 00913 #endif 00914 00915 int rank; 00916 int sc = MVT::GetNumberVecs( S ); 00917 ptrdiff_t sr = MVText::GetGlobalLength( S ); 00918 int numxy = X.length(); 00919 TEUCHOS_TEST_FOR_EXCEPTION(X.length() != Y.length(),std::invalid_argument, 00920 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X and Y must contain the same number of multivectors."); 00921 std::vector<int> xyc(numxy); 00922 // short-circuit 00923 if (sc == 0 || sr == 0) { 00924 #ifdef ANASAZI_ICGS_DEBUG 00925 *out << "Leaving Anasazi::ICGSOrthoManager::projectGen(...)\n"; 00926 #endif 00927 return 0; 00928 } 00929 // if we don't have enough C, expand it with null references 00930 // if we have too many, resize to throw away the latter ones 00931 // if we have exactly as many as we have X,Y this call has no effect 00932 // 00933 // same for MX, MY 00934 C.resize(numxy); 00935 MX.resize(numxy); 00936 MY.resize(numxy); 00937 00938 // check size of S w.r.t. common sense 00939 TEUCHOS_TEST_FOR_EXCEPTION( sc<0 || sr<0, std::invalid_argument, 00940 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MVT returned negative dimensions for S." ); 00941 00942 // check size of MS 00943 if (this->_hasOp == true) { 00944 if (MS != Teuchos::null) { 00945 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MS) != sr, std::invalid_argument, 00946 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS length not consistent with S." ); 00947 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*MS) != sc, std::invalid_argument, 00948 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): MS width not consistent with S." ); 00949 } 00950 } 00951 00952 // tally up size of all X,Y and check/allocate C 00953 ptrdiff_t sumxyc = 0; 00954 int MYmissing = 0; 00955 int MXmissing = 0; 00956 for (int i=0; i<numxy; i++) { 00957 if (X[i] != Teuchos::null && Y[i] != Teuchos::null) { 00958 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*X[i]) != sr, std::invalid_argument, 00959 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] length not consistent with S." ); 00960 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*Y[i]) != sr, std::invalid_argument, 00961 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Y[" << i << "] length not consistent with S." ); 00962 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*X[i]) != MVT::GetNumberVecs(*Y[i]), std::invalid_argument, 00963 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "] and Y[" << i << "] widths not consistent." ); 00964 00965 xyc[i] = MVT::GetNumberVecs( *X[i] ); 00966 TEUCHOS_TEST_FOR_EXCEPTION( sr < static_cast<ptrdiff_t>(xyc[i]), std::invalid_argument, 00967 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): X[" << i << "],Y[" << i << "] have less rows than columns, and therefore cannot be full rank." ); 00968 sumxyc += xyc[i]; 00969 00970 // check size of C[i] 00971 if ( C[i] == Teuchos::null ) { 00972 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xyc[i],sc) ); 00973 } 00974 else { 00975 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != xyc[i] || C[i]->numCols() != sc , std::invalid_argument, 00976 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of Q not consistent with size of C." ); 00977 } 00978 // check sizes of MX[i], MY[i] if present 00979 // if not present, count their absence 00980 if (MX[i] != Teuchos::null) { 00981 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MX[i]) != sr || MVT::GetNumberVecs(*MX[i]) != xyc[i], std::invalid_argument, 00982 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MX[" << i << "] not consistent with size of X[" << i << "]." ); 00983 } 00984 else { 00985 MXmissing += xyc[i]; 00986 } 00987 if (MY[i] != Teuchos::null) { 00988 TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*MY[i]) != sr || MVT::GetNumberVecs(*MY[i]) != xyc[i], std::invalid_argument, 00989 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of MY[" << i << "] not consistent with size of Y[" << i << "]." ); 00990 } 00991 else { 00992 MYmissing += xyc[i]; 00993 } 00994 } 00995 else { 00996 // if one is null and the other is not... the user may have made a mistake 00997 TEUCHOS_TEST_FOR_EXCEPTION(X[i] != Teuchos::null || Y[i] != Teuchos::null, std::invalid_argument, 00998 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): " 00999 << (X[i] == Teuchos::null ? "Y[" : "X[") << i << "] was provided but " 01000 << (X[i] == Teuchos::null ? "X[" : "Y[") << i << "] was not."); 01001 } 01002 } 01003 01004 // is this operation even feasible? 01005 TEUCHOS_TEST_FOR_EXCEPTION(sumxyc + sc > sr, std::invalid_argument, 01006 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): dimension of all X[i],Y[i] is " 01007 << sumxyc << " and requested " << sc << "-dimensional basis, but length of vectors is only " 01008 << sr << ". This is infeasible."); 01009 01010 01011 /****** DO NO MODIFY *MS IF _hasOp == false 01012 * if _hasOp == false, we don't need MS, MX or MY 01013 * 01014 * if _hasOp == true, we need MS (for S M-norms and normalization) and 01015 * therefore, we must also update MS, regardless of whether 01016 * it gets returned to the user (i.e., whether the user provided it) 01017 * we may need to allocate and compute MX or MY 01018 * 01019 * let BXY denote: 01020 * "X" - we have all M*X[i] 01021 * "Y" - we have all M*Y[i] 01022 * "B" - we have biorthogonality (for all X[i],Y[i]) 01023 * Encode these as values 01024 * X = 1 01025 * Y = 2 01026 * B = 4 01027 * We have 8 possibilities, 0-7 01028 * 01029 * We must allocate storage and computer the following, lest 01030 * innerProdMat do it for us: 01031 * none (0) - allocate MX, for inv(<X,Y>) and M*S 01032 * 01033 * for the following, we will compute M*S manually before returning 01034 * B(4) BY(6) Y(2) --> updateMS = 1 01035 * for the following, we will update M*S as we go, using M*X 01036 * XY(3) X(1) none(0) BXY(7) BX(5) --> updateMS = 2 01037 * these choices favor applications of M over allocation of memory 01038 * 01039 */ 01040 int updateMS = -1; 01041 if (this->_hasOp) { 01042 int whichAlloc = 0; 01043 if (MXmissing == 0) { 01044 whichAlloc += 1; 01045 } 01046 if (MYmissing == 0) { 01047 whichAlloc += 2; 01048 } 01049 if (isBiortho) { 01050 whichAlloc += 4; 01051 } 01052 01053 switch (whichAlloc) { 01054 case 2: 01055 case 4: 01056 case 6: 01057 updateMS = 1; 01058 break; 01059 case 0: 01060 case 1: 01061 case 3: 01062 case 5: 01063 case 7: 01064 updateMS = 2; 01065 break; 01066 } 01067 01068 // produce MS 01069 if (MS == Teuchos::null) { 01070 #ifdef ANASAZI_ICGS_DEBUG 01071 *out << "Allocating MS...\n"; 01072 #endif 01073 MS = MVT::Clone(S,MVT::GetNumberVecs(S)); 01074 OPT::Apply(*(this->_Op),S,*MS); 01075 this->_OpCounter += MVT::GetNumberVecs(S); 01076 } 01077 01078 // allocate the rest 01079 if (whichAlloc == 0) { 01080 // allocate and compute missing MX 01081 for (int i=0; i<numxy; i++) { 01082 if (MX[i] == Teuchos::null) { 01083 #ifdef ANASAZI_ICGS_DEBUG 01084 *out << "Allocating MX[" << i << "]...\n"; 01085 #endif 01086 Teuchos::RCP<MV> tmpMX = MVT::Clone(*X[i],xyc[i]); 01087 OPT::Apply(*(this->_Op),*X[i],*tmpMX); 01088 MX[i] = tmpMX; 01089 this->_OpCounter += xyc[i]; 01090 } 01091 } 01092 } 01093 } 01094 else { 01095 // Op == I --> MS == S 01096 MS = Teuchos::rcpFromRef(S); 01097 updateMS = 0; 01098 } 01099 TEUCHOS_TEST_FOR_EXCEPTION(updateMS == -1,std::logic_error, 01100 "Anasazi::ICGSOrthoManager::projectGen(): Error in updateMS logic."); 01101 01102 01103 // if the user doesn't want to store the coefficients, 01104 // allocate some local memory for them 01105 if ( B == Teuchos::null ) { 01106 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(sc,sc) ); 01107 } 01108 else { 01109 // check size of B 01110 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != sc || B->numCols() != sc, std::invalid_argument, 01111 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Size of S must be consistent with size of B" ); 01112 } 01113 01114 01115 // orthogonalize all of S against X,Y 01116 projectGen(S,X,Y,isBiortho,C,MS,MX,MY); 01117 01118 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(sc,1); 01119 // start working 01120 rank = 0; 01121 int numTries = 10; // each vector in X gets 10 random chances to escape degeneracy 01122 int oldrank = -1; 01123 do { 01124 int curssize = sc - rank; 01125 01126 // orthonormalize S, but quit if it is rank deficient 01127 // we can't let findBasis generated random vectors to complete the basis, 01128 // because it doesn't know about Q; we will do this ourselves below 01129 #ifdef ANASAZI_ICGS_DEBUG 01130 *out << "Attempting to find orthonormal basis for X...\n"; 01131 #endif 01132 rank = findBasis(S,MS,*B,false,curssize); 01133 01134 if (oldrank != -1 && rank != oldrank) { 01135 // we had previously stopped before, after operating on vector oldrank 01136 // we saved its coefficients, augmented it with a random vector, and 01137 // then called findBasis() again, which proceeded to add vector oldrank 01138 // to the basis. 01139 // now, restore the saved coefficients into B 01140 for (int i=0; i<sc; i++) { 01141 (*B)(i,oldrank) = oldCoeff(i,0); 01142 } 01143 } 01144 01145 if (rank < sc) { 01146 if (rank != oldrank) { 01147 // we quit on this vector and will augment it with random below 01148 // this is the first time that we have quit on this vector 01149 // therefor, (*B)(:,rank) contains the actual coefficients of the 01150 // input vectors with respect to the previous vectors in the basis 01151 // save these values, as (*B)(:,rank) will be overwritten by our next 01152 // call to findBasis() 01153 // we will restore it after we are done working on this vector 01154 for (int i=0; i<sc; i++) { 01155 oldCoeff(i,0) = (*B)(i,rank); 01156 } 01157 } 01158 } 01159 01160 if (rank == sc) { 01161 // we are done 01162 #ifdef ANASAZI_ICGS_DEBUG 01163 *out << "Finished computing basis.\n"; 01164 #endif 01165 break; 01166 } 01167 else { 01168 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError, 01169 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): basis lost rank; this shouldn't happen"); 01170 01171 if (rank != oldrank) { 01172 // we added a basis vector from random info; reset the chance counter 01173 numTries = 10; 01174 // store old rank 01175 oldrank = rank; 01176 } 01177 else { 01178 // has this vector run out of chances to escape degeneracy? 01179 if (numTries <= 0) { 01180 break; 01181 } 01182 } 01183 // use one of this vector's chances 01184 numTries--; 01185 01186 // randomize troubled direction 01187 #ifdef ANASAZI_ICGS_DEBUG 01188 *out << "Inserting random vector in X[" << rank << "]. Attempt " << 10-numTries << ".\n"; 01189 #endif 01190 Teuchos::RCP<MV> curS, curMS; 01191 std::vector<int> ind(1); 01192 ind[0] = rank; 01193 curS = MVT::CloneViewNonConst(S,ind); 01194 MVT::MvRandom(*curS); 01195 if (this->_hasOp) { 01196 #ifdef ANASAZI_ICGS_DEBUG 01197 *out << "Applying operator to random vector.\n"; 01198 #endif 01199 curMS = MVT::CloneViewNonConst(*MS,ind); 01200 OPT::Apply( *(this->_Op), *curS, *curMS ); 01201 this->_OpCounter += MVT::GetNumberVecs(*curS); 01202 } 01203 01204 // orthogonalize against X,Y 01205 // if !this->_hasOp, the curMS will be ignored. 01206 // we don't care about these coefficients 01207 // on the contrary, we need to preserve the previous coeffs 01208 { 01209 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0); 01210 projectGen(*curS,X,Y,isBiortho,dummyC,curMS,MX,MY); 01211 } 01212 } 01213 } while (1); 01214 01215 // this should never raise an exception; but our post-conditions oblige us to check 01216 TEUCHOS_TEST_FOR_EXCEPTION( rank > sc || rank < 0, std::logic_error, 01217 "Anasazi::ICGSOrthoManager::projectAndNormalizeGen(): Debug error in rank variable." ); 01218 01219 #ifdef ANASAZI_ICGS_DEBUG 01220 *out << "Returning " << rank << " from Anasazi::ICGSOrthoManager::projectAndNormalizeGen(...)\n"; 01221 #endif 01222 01223 return rank; 01224 } 01225 01226 01227 01229 // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 01230 // the rank is numvectors(X) 01231 template<class ScalarType, class MV, class OP> 01232 int ICGSOrthoManager<ScalarType, MV, OP>::findBasis( 01233 MV &X, Teuchos::RCP<MV> MX, 01234 Teuchos::SerialDenseMatrix<int,ScalarType> &B, 01235 bool completeBasis, int howMany ) const { 01236 01237 // For the inner product defined by the operator Op or the identity (Op == 0) 01238 // -> Orthonormalize X 01239 // Modify MX accordingly 01240 // 01241 // Note that when Op is 0, MX is not referenced 01242 // 01243 // Parameter variables 01244 // 01245 // X : Vectors to be orthonormalized 01246 // 01247 // MX : Image of the multivector X under the operator Op 01248 // 01249 // Op : Pointer to the operator for the inner product 01250 // 01251 01252 #ifdef ANASAZI_ICGS_DEBUG 01253 // Get a FancyOStream from out_arg or create a new one ... 01254 Teuchos::RCP<Teuchos::FancyOStream> 01255 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout)); 01256 out->setShowAllFrontMatter(false).setShowProcRank(true); 01257 *out << "Entering Anasazi::ICGSOrthoManager::findBasis(...)\n"; 01258 #endif 01259 01260 const ScalarType ONE = SCT::one(); 01261 const MagnitudeType ZERO = SCT::magnitude(SCT::zero()); 01262 01263 int xc = MVT::GetNumberVecs( X ); 01264 01265 if (howMany == -1) { 01266 howMany = xc; 01267 } 01268 01269 /******************************************************* 01270 * If _hasOp == false, we will not reference MX below * 01271 *******************************************************/ 01272 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error, 01273 "Anasazi::ICGSOrthoManager::findBasis(): calling routine did not specify MS."); 01274 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error, 01275 "Anasazi::ICGSOrthoManager::findBasis(): Invalid howMany parameter" ); 01276 01277 /* xstart is which column we are starting the process with, based on howMany 01278 * columns before xstart are assumed to be Op-orthonormal already 01279 */ 01280 int xstart = xc - howMany; 01281 01282 for (int j = xstart; j < xc; j++) { 01283 01284 // numX represents the number of currently orthonormal columns of X 01285 int numX = j; 01286 // j represents the index of the current column of X 01287 // these are different interpretations of the same value 01288 01289 // 01290 // set the lower triangular part of B to zero 01291 for (int i=j+1; i<xc; ++i) { 01292 B(i,j) = ZERO; 01293 } 01294 01295 // Get a view of the vector currently being worked on. 01296 std::vector<int> index(1); 01297 index[0] = j; 01298 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index ); 01299 Teuchos::RCP<MV> MXj; 01300 if ((this->_hasOp)) { 01301 // MXj is a view of the current vector in MX 01302 MXj = MVT::CloneViewNonConst( *MX, index ); 01303 } 01304 else { 01305 // MXj is a pointer to Xj, and MUST NOT be modified 01306 MXj = Xj; 01307 } 01308 01309 // Get a view of the previous vectors. 01310 std::vector<int> prev_idx( numX ); 01311 Teuchos::RCP<const MV> prevX, prevMX; 01312 01313 if (numX > 0) { 01314 for (int i=0; i<numX; ++i) prev_idx[i] = i; 01315 prevX = MVT::CloneView( X, prev_idx ); 01316 if (this->_hasOp) { 01317 prevMX = MVT::CloneView( *MX, prev_idx ); 01318 } 01319 } 01320 01321 bool rankDef = true; 01322 /* numTrials>0 will denote that the current vector was randomized for the purpose 01323 * of finding a basis vector, and that the coefficients of that vector should 01324 * not be stored in B 01325 */ 01326 for (int numTrials = 0; numTrials < 10; numTrials++) { 01327 #ifdef ANASAZI_ICGS_DEBUG 01328 *out << "Trial " << numTrials << " for vector " << j << "\n"; 01329 #endif 01330 01331 // Make storage for these Gram-Schmidt iterations. 01332 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1); 01333 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1); 01334 01335 // 01336 // Save old MXj vector and compute Op-norm 01337 // 01338 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 01339 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,origNorm,MXj); 01340 #ifdef ANASAZI_ICGS_DEBUG 01341 *out << "origNorm = " << origNorm[0] << "\n"; 01342 #endif 01343 01344 if (numX > 0) { 01345 // Apply the first step of Gram-Schmidt 01346 01347 // product <- prevX^T MXj 01348 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj); 01349 01350 // Xj <- Xj - prevX prevX^T MXj 01351 // = Xj - prevX product 01352 #ifdef ANASAZI_ICGS_DEBUG 01353 *out << "Orthogonalizing X[" << j << "]...\n"; 01354 #endif 01355 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj ); 01356 01357 // Update MXj 01358 if (this->_hasOp) { 01359 // MXj <- Op*Xj_new 01360 // = Op*(Xj_old - prevX prevX^T MXj) 01361 // = MXj - prevMX product 01362 #ifdef ANASAZI_ICGS_DEBUG 01363 *out << "Updating MX[" << j << "]...\n"; 01364 #endif 01365 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj ); 01366 } 01367 01368 // Compute new Op-norm 01369 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj); 01370 MagnitudeType product_norm = product.normOne(); 01371 01372 #ifdef ANASAZI_ICGS_DEBUG 01373 *out << "newNorm = " << newNorm[0] << "\n"; 01374 *out << "prodoct_norm = " << product_norm << "\n"; 01375 #endif 01376 01377 // Check if a correction is needed. 01378 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) { 01379 #ifdef ANASAZI_ICGS_DEBUG 01380 if (product_norm/newNorm[0] >= tol_) { 01381 *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n"; 01382 } 01383 else { 01384 *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n"; 01385 } 01386 #endif 01387 // Apply the second step of Gram-Schmidt 01388 // This is the same as above 01389 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1); 01390 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj); 01391 product += P2; 01392 #ifdef ANASAZI_ICGS_DEBUG 01393 *out << "Orthogonalizing X[" << j << "]...\n"; 01394 #endif 01395 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj ); 01396 if ((this->_hasOp)) { 01397 #ifdef ANASAZI_ICGS_DEBUG 01398 *out << "Updating MX[" << j << "]...\n"; 01399 #endif 01400 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj ); 01401 } 01402 // Compute new Op-norms 01403 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm2,MXj); 01404 product_norm = P2.normOne(); 01405 #ifdef ANASAZI_ICGS_DEBUG 01406 *out << "newNorm2 = " << newNorm2[0] << "\n"; 01407 *out << "product_norm = " << product_norm << "\n"; 01408 #endif 01409 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) { 01410 // we don't do another GS, we just set it to zero. 01411 #ifdef ANASAZI_ICGS_DEBUG 01412 if (product_norm/newNorm2[0] >= tol_) { 01413 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n"; 01414 } 01415 else if (newNorm[0] < newNorm2[0]) { 01416 *out << "newNorm2 > newNorm... setting vector to zero.\n"; 01417 } 01418 else { 01419 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n"; 01420 } 01421 #endif 01422 MVT::MvInit(*Xj,ZERO); 01423 if ((this->_hasOp)) { 01424 #ifdef ANASAZI_ICGS_DEBUG 01425 *out << "Setting MX[" << j << "] to zero as well.\n"; 01426 #endif 01427 MVT::MvInit(*MXj,ZERO); 01428 } 01429 } 01430 } 01431 } // if (numX > 0) do GS 01432 01433 // save the coefficients, if we are working on the original vector and not a randomly generated one 01434 if (numTrials == 0) { 01435 for (int i=0; i<numX; i++) { 01436 B(i,j) = product(i,0); 01437 } 01438 } 01439 01440 // Check if Xj has any directional information left after the orthogonalization. 01441 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj); 01442 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) { 01443 #ifdef ANASAZI_ICGS_DEBUG 01444 *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n"; 01445 #endif 01446 // Normalize Xj. 01447 // Xj <- Xj / norm 01448 MVT::MvScale( *Xj, ONE/newNorm[0]); 01449 if (this->_hasOp) { 01450 #ifdef ANASAZI_ICGS_DEBUG 01451 *out << "Normalizing M*X[" << j << "]...\n"; 01452 #endif 01453 // Update MXj. 01454 MVT::MvScale( *MXj, ONE/newNorm[0]); 01455 } 01456 01457 // save it, if it corresponds to the original vector and not a randomly generated one 01458 if (numTrials == 0) { 01459 B(j,j) = newNorm[0]; 01460 } 01461 01462 // We are not rank deficient in this vector. Move on to the next vector in X. 01463 rankDef = false; 01464 break; 01465 } 01466 else { 01467 #ifdef ANASAZI_ICGS_DEBUG 01468 *out << "Not normalizing M*X[" << j << "]...\n"; 01469 #endif 01470 // There was nothing left in Xj after orthogonalizing against previous columns in X. 01471 // X is rank deficient. 01472 // reflect this in the coefficients 01473 B(j,j) = ZERO; 01474 01475 if (completeBasis) { 01476 // Fill it with random information and keep going. 01477 #ifdef ANASAZI_ICGS_DEBUG 01478 *out << "Inserting random vector in X[" << j << "]...\n"; 01479 #endif 01480 MVT::MvRandom( *Xj ); 01481 if (this->_hasOp) { 01482 #ifdef ANASAZI_ICGS_DEBUG 01483 *out << "Updating M*X[" << j << "]...\n"; 01484 #endif 01485 OPT::Apply( *(this->_Op), *Xj, *MXj ); 01486 this->_OpCounter += MVT::GetNumberVecs(*Xj); 01487 } 01488 } 01489 else { 01490 rankDef = true; 01491 break; 01492 } 01493 } 01494 } // for (numTrials = 0; numTrials < 10; ++numTrials) 01495 01496 // if rankDef == true, then quit and notify user of rank obtained 01497 if (rankDef == true) { 01498 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError, 01499 "Anasazi::ICGSOrthoManager::findBasis(): Unable to complete basis" ); 01500 #ifdef ANASAZI_ICGS_DEBUG 01501 *out << "Returning early, rank " << j << " from Anasazi::ICGSOrthoManager::findBasis(...)\n"; 01502 #endif 01503 return j; 01504 } 01505 01506 } // for (j = 0; j < xc; ++j) 01507 01508 #ifdef ANASAZI_ICGS_DEBUG 01509 *out << "Returning " << xc << " from Anasazi::ICGSOrthoManager::findBasis(...)\n"; 01510 #endif 01511 return xc; 01512 } 01513 01514 } // namespace Anasazi 01515 01516 #endif // ANASAZI_ICSG_ORTHOMANAGER_HPP 01517
1.7.6.1
