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