|
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_BASIC_ORTHOMANAGER_HPP 00035 #define ANASAZI_BASIC_ORTHOMANAGER_HPP 00036 00044 #include "AnasaziConfigDefs.hpp" 00045 #include "AnasaziMultiVecTraits.hpp" 00046 #include "AnasaziOperatorTraits.hpp" 00047 #include "AnasaziMatOrthoManager.hpp" 00048 #include "Teuchos_TimeMonitor.hpp" 00049 #include "Teuchos_LAPACK.hpp" 00050 #include "Teuchos_BLAS.hpp" 00051 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00052 # include <Teuchos_FancyOStream.hpp> 00053 #endif 00054 00055 namespace Anasazi { 00056 00057 template<class ScalarType, class MV, class OP> 00058 class BasicOrthoManager : public MatOrthoManager<ScalarType,MV,OP> { 00059 00060 private: 00061 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType; 00062 typedef Teuchos::ScalarTraits<ScalarType> SCT; 00063 typedef MultiVecTraits<ScalarType,MV> MVT; 00064 typedef OperatorTraits<ScalarType,MV,OP> OPT; 00065 00066 public: 00067 00069 00070 00071 BasicOrthoManager( Teuchos::RCP<const OP> Op = Teuchos::null, 00072 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa = 1.41421356 /* sqrt(2) */, 00073 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps = 0.0, 00074 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol = 0.20 ); 00075 00076 00078 ~BasicOrthoManager() {} 00080 00081 00083 00084 00085 00124 void projectMat ( 00125 MV &X, 00126 Teuchos::Array<Teuchos::RCP<const MV> > Q, 00127 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C 00128 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)), 00129 Teuchos::RCP<MV> MX = Teuchos::null, 00130 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)) 00131 ) const; 00132 00133 00172 int normalizeMat ( 00173 MV &X, 00174 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null, 00175 Teuchos::RCP<MV> MX = Teuchos::null 00176 ) const; 00177 00178 00238 int projectAndNormalizeMat ( 00239 MV &X, 00240 Teuchos::Array<Teuchos::RCP<const MV> > Q, 00241 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C 00242 = Teuchos::tuple(Teuchos::RCP< Teuchos::SerialDenseMatrix<int,ScalarType> >(Teuchos::null)), 00243 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B = Teuchos::null, 00244 Teuchos::RCP<MV> MX = Teuchos::null, 00245 Teuchos::Array<Teuchos::RCP<const MV> > MQ = Teuchos::tuple(Teuchos::RCP<const MV>(Teuchos::null)) 00246 ) const; 00247 00249 00251 00252 00257 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 00258 orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX = Teuchos::null) const; 00259 00264 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 00265 orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const; 00266 00268 00270 00271 00273 void setKappa( typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa ) { kappa_ = kappa; } 00274 00276 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType getKappa() const { return kappa_; } 00277 00279 00280 private: 00282 MagnitudeType kappa_; 00283 MagnitudeType eps_; 00284 MagnitudeType tol_; 00285 00286 // ! Routine to find an orthonormal basis 00287 int findBasis(MV &X, Teuchos::RCP<MV> MX, 00288 Teuchos::SerialDenseMatrix<int,ScalarType> &B, 00289 bool completeBasis, int howMany = -1 ) const; 00290 00291 // 00292 // Internal timers 00293 // 00294 Teuchos::RCP<Teuchos::Time> timerReortho_; 00295 00296 }; 00297 00298 00300 // Constructor 00301 template<class ScalarType, class MV, class OP> 00302 BasicOrthoManager<ScalarType,MV,OP>::BasicOrthoManager( Teuchos::RCP<const OP> Op, 00303 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType kappa, 00304 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType eps, 00305 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType tol ) : 00306 MatOrthoManager<ScalarType,MV,OP>(Op), kappa_(kappa), eps_(eps), tol_(tol) 00307 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00308 , timerReortho_(Teuchos::TimeMonitor::getNewTimer("Anasazi::BasicOrthoManager::Re-orthogonalization")) 00309 #endif 00310 { 00311 TEUCHOS_TEST_FOR_EXCEPTION(eps_ < SCT::magnitude(SCT::zero()),std::invalid_argument, 00312 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"eps\" must be non-negative."); 00313 if (eps_ == 0) { 00314 Teuchos::LAPACK<int,MagnitudeType> lapack; 00315 eps_ = lapack.LAMCH('E'); 00316 eps_ = Teuchos::ScalarTraits<MagnitudeType>::pow(eps_,.75); 00317 } 00318 TEUCHOS_TEST_FOR_EXCEPTION( 00319 tol_ < SCT::magnitude(SCT::zero()) || tol_ > SCT::magnitude(SCT::one()), 00320 std::invalid_argument, 00321 "Anasazi::BasicOrthoManager::BasicOrthoManager(): argument \"tol\" must be in [0,1]."); 00322 } 00323 00324 00325 00327 // Compute the distance from orthonormality 00328 template<class ScalarType, class MV, class OP> 00329 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 00330 BasicOrthoManager<ScalarType,MV,OP>::orthonormErrorMat(const MV &X, Teuchos::RCP<const MV> MX) const { 00331 const ScalarType ONE = SCT::one(); 00332 int rank = MVT::GetNumberVecs(X); 00333 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(rank,rank); 00334 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X,X,xTx,MX,MX); 00335 for (int i=0; i<rank; i++) { 00336 xTx(i,i) -= ONE; 00337 } 00338 return xTx.normFrobenius(); 00339 } 00340 00341 00342 00344 // Compute the distance from orthogonality 00345 template<class ScalarType, class MV, class OP> 00346 typename Teuchos::ScalarTraits<ScalarType>::magnitudeType 00347 BasicOrthoManager<ScalarType,MV,OP>::orthogErrorMat(const MV &X1, const MV &X2, Teuchos::RCP<const MV> MX1, Teuchos::RCP<const MV> MX2) const { 00348 int r1 = MVT::GetNumberVecs(X1); 00349 int r2 = MVT::GetNumberVecs(X2); 00350 Teuchos::SerialDenseMatrix<int,ScalarType> xTx(r1,r2); 00351 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(X1,X2,xTx,MX1,MX2); 00352 return xTx.normFrobenius(); 00353 } 00354 00355 00356 00358 template<class ScalarType, class MV, class OP> 00359 void BasicOrthoManager<ScalarType, MV, OP>::projectMat( 00360 MV &X, 00361 Teuchos::Array<Teuchos::RCP<const MV> > Q, 00362 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 00363 Teuchos::RCP<MV> MX, 00364 Teuchos::Array<Teuchos::RCP<const MV> > MQ 00365 ) const { 00366 // For the inner product defined by the operator Op or the identity (Op == 0) 00367 // -> Orthogonalize X against each Q[i] 00368 // Modify MX accordingly 00369 // 00370 // Note that when Op is 0, MX is not referenced 00371 // 00372 // Parameter variables 00373 // 00374 // X : Vectors to be transformed 00375 // 00376 // MX : Image of the block vector X by the mass matrix 00377 // 00378 // Q : Bases to orthogonalize against. These are assumed orthonormal, mutually and independently. 00379 // 00380 00381 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00382 // Get a FancyOStream from out_arg or create a new one ... 00383 Teuchos::RCP<Teuchos::FancyOStream> 00384 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout)); 00385 out->setShowAllFrontMatter(false).setShowProcRank(true); 00386 *out << "Entering Anasazi::BasicOrthoManager::projectMat(...)\n"; 00387 #endif 00388 00389 ScalarType ONE = SCT::one(); 00390 00391 int xc = MVT::GetNumberVecs( X ); 00392 int xr = MVT::GetVecLength( X ); 00393 int nq = Q.length(); 00394 std::vector<int> qcs(nq); 00395 // short-circuit 00396 if (nq == 0 || xc == 0 || xr == 0) { 00397 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00398 *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n"; 00399 #endif 00400 return; 00401 } 00402 int qr = MVT::GetVecLength ( *Q[0] ); 00403 // if we don't have enough C, expand it with null references 00404 // if we have too many, resize to throw away the latter ones 00405 // if we have exactly as many as we have Q, this call has no effect 00406 C.resize(nq); 00407 00408 00409 /****** DO NO MODIFY *MX IF _hasOp == false ******/ 00410 if (this->_hasOp) { 00411 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00412 *out << "Allocating MX...\n"; 00413 #endif 00414 if (MX == Teuchos::null) { 00415 // we need to allocate space for MX 00416 MX = MVT::Clone(X,MVT::GetNumberVecs(X)); 00417 OPT::Apply(*(this->_Op),X,*MX); 00418 this->_OpCounter += MVT::GetNumberVecs(X); 00419 } 00420 } 00421 else { 00422 // Op == I --> MX = X (ignore it if the user passed it in) 00423 MX = Teuchos::rcpFromRef(X); 00424 } 00425 int mxc = MVT::GetNumberVecs( *MX ); 00426 int mxr = MVT::GetVecLength( *MX ); 00427 00428 // check size of X and Q w.r.t. common sense 00429 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 00430 "Anasazi::BasicOrthoManager::projectMat(): MVT returned negative dimensions for X,MX" ); 00431 // check size of X w.r.t. MX and Q 00432 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr || xr!=qr, std::invalid_argument, 00433 "Anasazi::BasicOrthoManager::projectMat(): Size of X not consistent with MX,Q" ); 00434 00435 // tally up size of all Q and check/allocate C 00436 int baslen = 0; 00437 for (int i=0; i<nq; i++) { 00438 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength( *Q[i] ) != qr, std::invalid_argument, 00439 "Anasazi::BasicOrthoManager::projectMat(): Q lengths not mutually consistent" ); 00440 qcs[i] = MVT::GetNumberVecs( *Q[i] ); 00441 TEUCHOS_TEST_FOR_EXCEPTION( qr < qcs[i], std::invalid_argument, 00442 "Anasazi::BasicOrthoManager::projectMat(): Q has less rows than columns" ); 00443 baslen += qcs[i]; 00444 00445 // check size of C[i] 00446 if ( C[i] == Teuchos::null ) { 00447 C[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(qcs[i],xc) ); 00448 } 00449 else { 00450 TEUCHOS_TEST_FOR_EXCEPTION( C[i]->numRows() != qcs[i] || C[i]->numCols() != xc , std::invalid_argument, 00451 "Anasazi::BasicOrthoManager::projectMat(): Size of Q not consistent with size of C" ); 00452 } 00453 } 00454 00455 // Perform the Gram-Schmidt transformation for a block of vectors 00456 00457 // Compute the initial Op-norms 00458 std::vector<ScalarType> oldDot( xc ); 00459 MVT::MvDot( X, *MX, oldDot ); 00460 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00461 *out << "oldDot = { "; 00462 std::copy(oldDot.begin(), oldDot.end(), std::ostream_iterator<ScalarType>(*out, " ")); 00463 *out << "}\n"; 00464 #endif 00465 00466 MQ.resize(nq); 00467 // Define the product Q^T * (Op*X) 00468 for (int i=0; i<nq; i++) { 00469 // Multiply Q' with MX 00470 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,*C[i],MQ[i],MX); 00471 // Multiply by Q and subtract the result in X 00472 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00473 *out << "Applying projector P_Q[" << i << "]...\n"; 00474 #endif 00475 MVT::MvTimesMatAddMv( -ONE, *Q[i], *C[i], ONE, X ); 00476 00477 // Update MX, with the least number of applications of Op as possible 00478 // Update MX. If we have MQ, use it. Otherwise, just multiply by Op 00479 if (this->_hasOp) { 00480 if (MQ[i] == Teuchos::null) { 00481 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00482 *out << "Updating MX via M*X...\n"; 00483 #endif 00484 OPT::Apply( *(this->_Op), X, *MX); 00485 this->_OpCounter += MVT::GetNumberVecs(X); 00486 } 00487 else { 00488 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00489 *out << "Updating MX via M*Q...\n"; 00490 #endif 00491 MVT::MvTimesMatAddMv( -ONE, *MQ[i], *C[i], ONE, *MX ); 00492 } 00493 } 00494 } 00495 00496 // Compute new Op-norms 00497 std::vector<ScalarType> newDot(xc); 00498 MVT::MvDot( X, *MX, newDot ); 00499 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00500 *out << "newDot = { "; 00501 std::copy(newDot.begin(), newDot.end(), std::ostream_iterator<ScalarType>(*out, " ")); 00502 *out << "}\n"; 00503 #endif 00504 00505 // determine (individually) whether to do another step of classical Gram-Schmidt 00506 for (int j = 0; j < xc; ++j) { 00507 00508 if ( SCT::magnitude(kappa_*newDot[j]) < SCT::magnitude(oldDot[j]) ) { 00509 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00510 *out << "kappa_*newDot[" <<j<< "] == " << kappa_*newDot[j] << "... another step of Gram-Schmidt.\n"; 00511 #endif 00512 #ifdef ANASAZI_TEUCHOS_TIME_MONITOR 00513 Teuchos::TimeMonitor lcltimer( *timerReortho_ ); 00514 #endif 00515 for (int i=0; i<nq; i++) { 00516 Teuchos::SerialDenseMatrix<int,ScalarType> C2(*C[i]); 00517 00518 // Apply another step of classical Gram-Schmidt 00519 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*Q[i],X,C2,MQ[i],MX); 00520 *C[i] += C2; 00521 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00522 *out << "Applying projector P_Q[" << i << "]...\n"; 00523 #endif 00524 MVT::MvTimesMatAddMv( -ONE, *Q[i], C2, ONE, X ); 00525 00526 // Update MX as above 00527 if (this->_hasOp) { 00528 if (MQ[i] == Teuchos::null) { 00529 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00530 *out << "Updating MX via M*X...\n"; 00531 #endif 00532 OPT::Apply( *(this->_Op), X, *MX); 00533 this->_OpCounter += MVT::GetNumberVecs(X); 00534 } 00535 else { 00536 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00537 *out << "Updating MX via M*Q...\n"; 00538 #endif 00539 MVT::MvTimesMatAddMv( -ONE, *MQ[i], C2, ONE, *MX ); 00540 } 00541 } 00542 } 00543 break; 00544 } // if (kappa_*newDot[j] < oldDot[j]) 00545 } // for (int j = 0; j < xc; ++j) 00546 00547 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00548 *out << "Leaving Anasazi::BasicOrthoManager::projectMat(...)\n"; 00549 #endif 00550 } 00551 00552 00553 00555 // Find an Op-orthonormal basis for span(X), with rank numvectors(X) 00556 template<class ScalarType, class MV, class OP> 00557 int BasicOrthoManager<ScalarType, MV, OP>::normalizeMat( 00558 MV &X, 00559 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 00560 Teuchos::RCP<MV> MX) const { 00561 // call findBasis(), with the instruction to try to generate a basis of rank numvecs(X) 00562 // findBasis() requires MX 00563 00564 int xc = MVT::GetNumberVecs(X); 00565 int xr = MVT::GetVecLength(X); 00566 00567 // if Op==null, MX == X (via pointer) 00568 // Otherwise, either the user passed in MX or we will allocated and compute it 00569 if (this->_hasOp) { 00570 if (MX == Teuchos::null) { 00571 // we need to allocate space for MX 00572 MX = MVT::Clone(X,xc); 00573 OPT::Apply(*(this->_Op),X,*MX); 00574 this->_OpCounter += MVT::GetNumberVecs(X); 00575 } 00576 } 00577 00578 // if the user doesn't want to store the coefficients, 00579 // allocate some local memory for them 00580 if ( B == Teuchos::null ) { 00581 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) ); 00582 } 00583 00584 int mxc = (this->_hasOp) ? MVT::GetNumberVecs( *MX ) : xc; 00585 int mxr = (this->_hasOp) ? MVT::GetVecLength( *MX ) : xr; 00586 00587 // check size of C, B 00588 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, 00589 "Anasazi::BasicOrthoManager::normalizeMat(): X must be non-empty" ); 00590 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 00591 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of B" ); 00592 TEUCHOS_TEST_FOR_EXCEPTION( xc != mxc || xr != mxr, std::invalid_argument, 00593 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not consistent with size of MX" ); 00594 TEUCHOS_TEST_FOR_EXCEPTION( xc > xr, std::invalid_argument, 00595 "Anasazi::BasicOrthoManager::normalizeMat(): Size of X not feasible for normalization" ); 00596 00597 return findBasis(X, MX, *B, true ); 00598 } 00599 00600 00601 00603 // Find an Op-orthonormal basis for span(X) - span(W) 00604 template<class ScalarType, class MV, class OP> 00605 int BasicOrthoManager<ScalarType, MV, OP>::projectAndNormalizeMat( 00606 MV &X, 00607 Teuchos::Array<Teuchos::RCP<const MV> > Q, 00608 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > C, 00609 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B, 00610 Teuchos::RCP<MV> MX, 00611 Teuchos::Array<Teuchos::RCP<const MV> > MQ 00612 ) const { 00613 00614 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00615 // Get a FancyOStream from out_arg or create a new one ... 00616 Teuchos::RCP<Teuchos::FancyOStream> 00617 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout)); 00618 out->setShowAllFrontMatter(false).setShowProcRank(true); 00619 *out << "Entering Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n"; 00620 #endif 00621 00622 int nq = Q.length(); 00623 int xc = MVT::GetNumberVecs( X ); 00624 int xr = MVT::GetVecLength( X ); 00625 int rank; 00626 00627 /* if the user doesn't want to store the coefficients, 00628 * allocate some local memory for them 00629 */ 00630 if ( B == Teuchos::null ) { 00631 B = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(xc,xc) ); 00632 } 00633 00634 /****** DO NO MODIFY *MX IF _hasOp == false ******/ 00635 if (this->_hasOp) { 00636 if (MX == Teuchos::null) { 00637 // we need to allocate space for MX 00638 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00639 *out << "Allocating MX...\n"; 00640 #endif 00641 MX = MVT::Clone(X,MVT::GetNumberVecs(X)); 00642 OPT::Apply(*(this->_Op),X,*MX); 00643 this->_OpCounter += MVT::GetNumberVecs(X); 00644 } 00645 } 00646 else { 00647 // Op == I --> MX = X (ignore it if the user passed it in) 00648 MX = Teuchos::rcpFromRef(X); 00649 } 00650 00651 int mxc = MVT::GetNumberVecs( *MX ); 00652 int mxr = MVT::GetVecLength( *MX ); 00653 00654 TEUCHOS_TEST_FOR_EXCEPTION( xc == 0 || xr == 0, std::invalid_argument, "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): X must be non-empty" ); 00655 00656 int numbas = 0; 00657 for (int i=0; i<nq; i++) { 00658 numbas += MVT::GetNumberVecs( *Q[i] ); 00659 } 00660 00661 // check size of B 00662 TEUCHOS_TEST_FOR_EXCEPTION( B->numRows() != xc || B->numCols() != xc, std::invalid_argument, 00663 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of B" ); 00664 // check size of X and MX 00665 TEUCHOS_TEST_FOR_EXCEPTION( xc<0 || xr<0 || mxc<0 || mxr<0, std::invalid_argument, 00666 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): MVT returned negative dimensions for X,MX" ); 00667 // check size of X w.r.t. MX 00668 TEUCHOS_TEST_FOR_EXCEPTION( xc!=mxc || xr!=mxr, std::invalid_argument, 00669 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Size of X must be consistent with size of MX" ); 00670 // check feasibility 00671 TEUCHOS_TEST_FOR_EXCEPTION( numbas+xc > xr, std::invalid_argument, 00672 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Orthogonality constraints not feasible" ); 00673 00674 // orthogonalize all of X against Q 00675 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00676 *out << "Orthogonalizing X against Q...\n"; 00677 #endif 00678 projectMat(X,Q,C,MX,MQ); 00679 00680 Teuchos::SerialDenseMatrix<int,ScalarType> oldCoeff(xc,1); 00681 // start working 00682 rank = 0; 00683 int numTries = 10; // each vector in X gets 10 random chances to escape degeneracy 00684 int oldrank = -1; 00685 do { 00686 int curxsize = xc - rank; 00687 00688 // orthonormalize X, but quit if it is rank deficient 00689 // we can't let findBasis generated random vectors to complete the basis, 00690 // because it doesn't know about Q; we will do this ourselves below 00691 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00692 *out << "Attempting to find orthonormal basis for X...\n"; 00693 #endif 00694 rank = findBasis(X,MX,*B,false,curxsize); 00695 00696 if (oldrank != -1 && rank != oldrank) { 00697 // we had previously stopped before, after operating on vector oldrank 00698 // we saved its coefficients, augmented it with a random vector, and 00699 // then called findBasis() again, which proceeded to add vector oldrank 00700 // to the basis. 00701 // now, restore the saved coefficients into B 00702 for (int i=0; i<xc; i++) { 00703 (*B)(i,oldrank) = oldCoeff(i,0); 00704 } 00705 } 00706 00707 if (rank < xc) { 00708 if (rank != oldrank) { 00709 // we quit on this vector and will augment it with random below 00710 // this is the first time that we have quit on this vector 00711 // therefor, (*B)(:,rank) contains the actual coefficients of the 00712 // input vectors with respect to the previous vectors in the basis 00713 // save these values, as (*B)(:,rank) will be overwritten by our next 00714 // call to findBasis() 00715 // we will restore it after we are done working on this vector 00716 for (int i=0; i<xc; i++) { 00717 oldCoeff(i,0) = (*B)(i,rank); 00718 } 00719 } 00720 } 00721 00722 if (rank == xc) { 00723 // we are done 00724 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00725 *out << "Finished computing basis.\n"; 00726 #endif 00727 break; 00728 } 00729 else { 00730 TEUCHOS_TEST_FOR_EXCEPTION( rank < oldrank, OrthoError, 00731 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): basis lost rank; this shouldn't happen"); 00732 00733 if (rank != oldrank) { 00734 // we added a vector to the basis; reset the chance counter 00735 numTries = 10; 00736 // store old rank 00737 oldrank = rank; 00738 } 00739 else { 00740 // has this vector run out of chances to escape degeneracy? 00741 if (numTries <= 0) { 00742 break; 00743 } 00744 } 00745 // use one of this vector's chances 00746 numTries--; 00747 00748 // randomize troubled direction 00749 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00750 *out << "Randomizing X[" << rank << "]...\n"; 00751 #endif 00752 Teuchos::RCP<MV> curX, curMX; 00753 std::vector<int> ind(1); 00754 ind[0] = rank; 00755 curX = MVT::CloneViewNonConst(X,ind); 00756 MVT::MvRandom(*curX); 00757 if (this->_hasOp) { 00758 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00759 *out << "Applying operator to random vector.\n"; 00760 #endif 00761 curMX = MVT::CloneViewNonConst(*MX,ind); 00762 OPT::Apply( *(this->_Op), *curX, *curMX ); 00763 this->_OpCounter += MVT::GetNumberVecs(*curX); 00764 } 00765 00766 // orthogonalize against Q 00767 // if !this->_hasOp, the curMX will be ignored. 00768 // we don't care about these coefficients 00769 // on the contrary, we need to preserve the previous coeffs 00770 { 00771 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC(0); 00772 projectMat(*curX,Q,dummyC,curMX,MQ); 00773 } 00774 } 00775 } while (1); 00776 00777 // this should never raise an exception; but our post-conditions oblige us to check 00778 TEUCHOS_TEST_FOR_EXCEPTION( rank > xc || rank < 0, std::logic_error, 00779 "Anasazi::BasicOrthoManager::projectAndNormalizeMat(): Debug error in rank variable." ); 00780 00781 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00782 *out << "Leaving Anasazi::BasicOrthoManager::projectAndNormalizeMat(...)\n"; 00783 #endif 00784 00785 return rank; 00786 } 00787 00788 00789 00791 // Find an Op-orthonormal basis for span(X), with the option of extending the subspace so that 00792 // the rank is numvectors(X) 00793 template<class ScalarType, class MV, class OP> 00794 int BasicOrthoManager<ScalarType, MV, OP>::findBasis( 00795 MV &X, Teuchos::RCP<MV> MX, 00796 Teuchos::SerialDenseMatrix<int,ScalarType> &B, 00797 bool completeBasis, int howMany ) const { 00798 00799 // For the inner product defined by the operator Op or the identity (Op == 0) 00800 // -> Orthonormalize X 00801 // Modify MX accordingly 00802 // 00803 // Note that when Op is 0, MX is not referenced 00804 // 00805 // Parameter variables 00806 // 00807 // X : Vectors to be orthonormalized 00808 // 00809 // MX : Image of the multivector X under the operator Op 00810 // 00811 // Op : Pointer to the operator for the inner product 00812 // 00813 00814 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00815 // Get a FancyOStream from out_arg or create a new one ... 00816 Teuchos::RCP<Teuchos::FancyOStream> 00817 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout)); 00818 out->setShowAllFrontMatter(false).setShowProcRank(true); 00819 *out << "Entering Anasazi::BasicOrthoManager::findBasis(...)\n"; 00820 #endif 00821 00822 const ScalarType ONE = SCT::one(); 00823 const MagnitudeType ZERO = SCT::magnitude(SCT::zero()); 00824 00825 int xc = MVT::GetNumberVecs( X ); 00826 00827 if (howMany == -1) { 00828 howMany = xc; 00829 } 00830 00831 /******************************************************* 00832 * If _hasOp == false, we will not reference MX below * 00833 *******************************************************/ 00834 TEUCHOS_TEST_FOR_EXCEPTION(this->_hasOp == true && MX == Teuchos::null, std::logic_error, 00835 "Anasazi::BasicOrthoManager::findBasis(): calling routine did not specify MS."); 00836 TEUCHOS_TEST_FOR_EXCEPTION( howMany < 0 || howMany > xc, std::logic_error, 00837 "Anasazi::BasicOrthoManager::findBasis(): Invalid howMany parameter" ); 00838 00839 /* xstart is which column we are starting the process with, based on howMany 00840 * columns before xstart are assumed to be Op-orthonormal already 00841 */ 00842 int xstart = xc - howMany; 00843 00844 for (int j = xstart; j < xc; j++) { 00845 00846 // numX represents the number of currently orthonormal columns of X 00847 int numX = j; 00848 // j represents the index of the current column of X 00849 // these are different interpretations of the same value 00850 00851 // 00852 // set the lower triangular part of B to zero 00853 for (int i=j+1; i<xc; ++i) { 00854 B(i,j) = ZERO; 00855 } 00856 00857 // Get a view of the vector currently being worked on. 00858 std::vector<int> index(1); 00859 index[0] = j; 00860 Teuchos::RCP<MV> Xj = MVT::CloneViewNonConst( X, index ); 00861 Teuchos::RCP<MV> MXj; 00862 if ((this->_hasOp)) { 00863 // MXj is a view of the current vector in MX 00864 MXj = MVT::CloneViewNonConst( *MX, index ); 00865 } 00866 else { 00867 // MXj is a pointer to Xj, and MUST NOT be modified 00868 MXj = Xj; 00869 } 00870 00871 // Get a view of the previous vectors. 00872 std::vector<int> prev_idx( numX ); 00873 Teuchos::RCP<const MV> prevX, prevMX; 00874 00875 if (numX > 0) { 00876 for (int i=0; i<numX; ++i) prev_idx[i] = i; 00877 prevX = MVT::CloneViewNonConst( X, prev_idx ); 00878 if (this->_hasOp) { 00879 prevMX = MVT::CloneViewNonConst( *MX, prev_idx ); 00880 } 00881 } 00882 00883 bool rankDef = true; 00884 /* numTrials>0 will denote that the current vector was randomized for the purpose 00885 * of finding a basis vector, and that the coefficients of that vector should 00886 * not be stored in B 00887 */ 00888 for (int numTrials = 0; numTrials < 10; numTrials++) { 00889 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00890 *out << "Trial " << numTrials << " for vector " << j << "\n"; 00891 #endif 00892 00893 // Make storage for these Gram-Schmidt iterations. 00894 Teuchos::SerialDenseMatrix<int,ScalarType> product(numX, 1); 00895 std::vector<MagnitudeType> origNorm(1), newNorm(1), newNorm2(1); 00896 00897 // 00898 // Save old MXj vector and compute Op-norm 00899 // 00900 Teuchos::RCP<MV> oldMXj = MVT::CloneCopy( *MXj ); 00901 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,origNorm,MXj); 00902 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00903 *out << "origNorm = " << origNorm[0] << "\n"; 00904 #endif 00905 00906 if (numX > 0) { 00907 // Apply the first step of Gram-Schmidt 00908 00909 // product <- prevX^T MXj 00910 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,product,Teuchos::null,MXj); 00911 00912 // Xj <- Xj - prevX prevX^T MXj 00913 // = Xj - prevX product 00914 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00915 *out << "Orthogonalizing X[" << j << "]...\n"; 00916 #endif 00917 MVT::MvTimesMatAddMv( -ONE, *prevX, product, ONE, *Xj ); 00918 00919 // Update MXj 00920 if (this->_hasOp) { 00921 // MXj <- Op*Xj_new 00922 // = Op*(Xj_old - prevX prevX^T MXj) 00923 // = MXj - prevMX product 00924 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00925 *out << "Updating MX[" << j << "]...\n"; 00926 #endif 00927 MVT::MvTimesMatAddMv( -ONE, *prevMX, product, ONE, *MXj ); 00928 } 00929 00930 // Compute new Op-norm 00931 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj); 00932 MagnitudeType product_norm = product.normOne(); 00933 00934 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00935 *out << "newNorm = " << newNorm[0] << "\n"; 00936 *out << "prodoct_norm = " << product_norm << "\n"; 00937 #endif 00938 00939 // Check if a correction is needed. 00940 if ( product_norm/newNorm[0] >= tol_ || newNorm[0] < eps_*origNorm[0]) { 00941 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00942 if (product_norm/newNorm[0] >= tol_) { 00943 *out << "product_norm/newNorm == " << product_norm/newNorm[0] << "... another step of Gram-Schmidt.\n"; 00944 } 00945 else { 00946 *out << "eps*origNorm == " << eps_*origNorm[0] << "... another step of Gram-Schmidt.\n"; 00947 } 00948 #endif 00949 // Apply the second step of Gram-Schmidt 00950 // This is the same as above 00951 Teuchos::SerialDenseMatrix<int,ScalarType> P2(numX,1); 00952 MatOrthoManager<ScalarType,MV,OP>::innerProdMat(*prevX,*Xj,P2,Teuchos::null,MXj); 00953 product += P2; 00954 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00955 *out << "Orthogonalizing X[" << j << "]...\n"; 00956 #endif 00957 MVT::MvTimesMatAddMv( -ONE, *prevX, P2, ONE, *Xj ); 00958 if ((this->_hasOp)) { 00959 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00960 *out << "Updating MX[" << j << "]...\n"; 00961 #endif 00962 MVT::MvTimesMatAddMv( -ONE, *prevMX, P2, ONE, *MXj ); 00963 } 00964 // Compute new Op-norms 00965 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm2,MXj); 00966 product_norm = P2.normOne(); 00967 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00968 *out << "newNorm2 = " << newNorm2[0] << "\n"; 00969 *out << "product_norm = " << product_norm << "\n"; 00970 #endif 00971 if ( product_norm/newNorm2[0] >= tol_ || newNorm2[0] < eps_*origNorm[0] ) { 00972 // we don't do another GS, we just set it to zero. 00973 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00974 if (product_norm/newNorm2[0] >= tol_) { 00975 *out << "product_norm/newNorm2 == " << product_norm/newNorm2[0] << "... setting vector to zero.\n"; 00976 } 00977 else if (newNorm[0] < newNorm2[0]) { 00978 *out << "newNorm2 > newNorm... setting vector to zero.\n"; 00979 } 00980 else { 00981 *out << "eps*origNorm == " << eps_*origNorm[0] << "... setting vector to zero.\n"; 00982 } 00983 #endif 00984 MVT::MvInit(*Xj,ZERO); 00985 if ((this->_hasOp)) { 00986 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 00987 *out << "Setting MX[" << j << "] to zero as well.\n"; 00988 #endif 00989 MVT::MvInit(*MXj,ZERO); 00990 } 00991 } 00992 } 00993 } // if (numX > 0) do GS 00994 00995 // save the coefficients, if we are working on the original vector and not a randomly generated one 00996 if (numTrials == 0) { 00997 for (int i=0; i<numX; i++) { 00998 B(i,j) = product(i,0); 00999 } 01000 } 01001 01002 // Check if Xj has any directional information left after the orthogonalization. 01003 MatOrthoManager<ScalarType,MV,OP>::normMat(*Xj,newNorm,MXj); 01004 if ( newNorm[0] != ZERO && newNorm[0] > SCT::sfmin() ) { 01005 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 01006 *out << "Normalizing X[" << j << "], norm(X[" << j << "]) = " << newNorm[0] << "\n"; 01007 #endif 01008 // Normalize Xj. 01009 // Xj <- Xj / norm 01010 MVT::MvScale( *Xj, ONE/newNorm[0]); 01011 if (this->_hasOp) { 01012 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 01013 *out << "Normalizing M*X[" << j << "]...\n"; 01014 #endif 01015 // Update MXj. 01016 MVT::MvScale( *MXj, ONE/newNorm[0]); 01017 } 01018 01019 // save it, if it corresponds to the original vector and not a randomly generated one 01020 if (numTrials == 0) { 01021 B(j,j) = newNorm[0]; 01022 } 01023 01024 // We are not rank deficient in this vector. Move on to the next vector in X. 01025 rankDef = false; 01026 break; 01027 } 01028 else { 01029 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 01030 *out << "Not normalizing M*X[" << j << "]...\n"; 01031 #endif 01032 // There was nothing left in Xj after orthogonalizing against previous columns in X. 01033 // X is rank deficient. 01034 // reflect this in the coefficients 01035 B(j,j) = ZERO; 01036 01037 if (completeBasis) { 01038 // Fill it with random information and keep going. 01039 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 01040 *out << "Inserting random vector in X[" << j << "]...\n"; 01041 #endif 01042 MVT::MvRandom( *Xj ); 01043 if (this->_hasOp) { 01044 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 01045 *out << "Updating M*X[" << j << "]...\n"; 01046 #endif 01047 OPT::Apply( *(this->_Op), *Xj, *MXj ); 01048 this->_OpCounter += MVT::GetNumberVecs(*Xj); 01049 } 01050 } 01051 else { 01052 rankDef = true; 01053 break; 01054 } 01055 } 01056 } // for (numTrials = 0; numTrials < 10; ++numTrials) 01057 01058 // if rankDef == true, then quit and notify user of rank obtained 01059 if (rankDef == true) { 01060 TEUCHOS_TEST_FOR_EXCEPTION( rankDef && completeBasis, OrthoError, 01061 "Anasazi::BasicOrthoManager::findBasis(): Unable to complete basis" ); 01062 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 01063 *out << "Returning early, rank " << j << " from Anasazi::BasicOrthoManager::findBasis(...)\n"; 01064 #endif 01065 return j; 01066 } 01067 01068 } // for (j = 0; j < xc; ++j) 01069 01070 #ifdef ANASAZI_BASIC_ORTHO_DEBUG 01071 *out << "Returning " << xc << " from Anasazi::BasicOrthoManager::findBasis(...)\n"; 01072 #endif 01073 return xc; 01074 } 01075 01076 } // namespace Anasazi 01077 01078 #endif // ANASAZI_BASIC_ORTHOMANAGER_HPP 01079
1.7.6.1
