00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #ifndef ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
00030 #define ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
00031
00036 #include "AnasaziConfigDefs.hpp"
00037 #include "AnasaziTypes.hpp"
00038
00039 #include "AnasaziEigenproblem.hpp"
00040 #include "AnasaziSolverManager.hpp"
00041 #include "AnasaziSolverUtils.hpp"
00042
00043 #include "AnasaziBlockDavidson.hpp"
00044 #include "AnasaziBasicSort.hpp"
00045 #include "AnasaziSVQBOrthoManager.hpp"
00046 #include "AnasaziBasicOrthoManager.hpp"
00047 #include "AnasaziStatusTestResNorm.hpp"
00048 #include "AnasaziStatusTestWithOrdering.hpp"
00049 #include "AnasaziStatusTestCombo.hpp"
00050 #include "AnasaziStatusTestOutput.hpp"
00051 #include "AnasaziBasicOutputManager.hpp"
00052 #include "Teuchos_BLAS.hpp"
00053 #include "Teuchos_LAPACK.hpp"
00054 #include "Teuchos_TimeMonitor.hpp"
00055 #ifdef TEUCHOS_DEBUG
00056 # include <Teuchos_FancyOStream.hpp>
00057 #endif
00058 #ifdef HAVE_MPI
00059 #include <mpi.h>
00060 #endif
00061
00062
00063
00077 namespace Anasazi {
00078
00079
00112 template<class ScalarType, class MV, class OP>
00113 class BlockDavidsonSolMgr : public SolverManager<ScalarType,MV,OP> {
00114
00115 private:
00116 typedef MultiVecTraits<ScalarType,MV> MVT;
00117 typedef OperatorTraits<ScalarType,MV,OP> OPT;
00118 typedef Teuchos::ScalarTraits<ScalarType> SCT;
00119 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
00120 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
00121
00122 public:
00123
00125
00126
00149 BlockDavidsonSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00150 Teuchos::ParameterList &pl );
00151
00153 virtual ~BlockDavidsonSolMgr() {};
00155
00157
00158
00160 const Eigenproblem<ScalarType,MV,OP>& getProblem() const {
00161 return *problem_;
00162 }
00163
00165 int getNumIters() const {
00166 return numIters_;
00167 }
00168
00176 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
00177 return tuple(_timerSolve, _timerRestarting, _timerLocking);
00178 }
00179
00181
00183
00184
00206 ReturnType solve();
00207
00209 void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
00210
00212 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
00213
00215 void setLockingStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking);
00216
00218 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
00219
00221 void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
00222
00224 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
00225
00227
00228 private:
00229 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
00230
00231 std::string whch_, ortho_;
00232
00233 MagnitudeType convtol_, locktol_;
00234 int maxRestarts_;
00235 bool useLocking_;
00236 bool relconvtol_, rellocktol_;
00237 int blockSize_, numBlocks_, numIters_;
00238 int maxLocked_;
00239 int lockQuorum_;
00240 bool inSituRestart_;
00241 int numRestartBlocks_;
00242 enum StatusTestResNorm<ScalarType,MV,OP>::ResType convNorm_, lockNorm_;
00243
00244 Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
00245
00246 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
00247 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
00248 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
00249
00250 Teuchos::RCP<BasicOutputManager<ScalarType> > printer_;
00251 };
00252
00253
00254
00255 template<class ScalarType, class MV, class OP>
00256 BlockDavidsonSolMgr<ScalarType,MV,OP>::BlockDavidsonSolMgr(
00257 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
00258 Teuchos::ParameterList &pl ) :
00259 problem_(problem),
00260 whch_("SR"),
00261 ortho_("SVQB"),
00262 convtol_(MT::prec()),
00263 maxRestarts_(20),
00264 useLocking_(false),
00265 relconvtol_(true),
00266 rellocktol_(true),
00267 blockSize_(0),
00268 numBlocks_(0),
00269 numIters_(0),
00270 maxLocked_(0),
00271 lockQuorum_(1),
00272 inSituRestart_(false),
00273 numRestartBlocks_(1),
00274 _timerSolve(Teuchos::TimeMonitor::getNewTimer("BlockDavidsonSolMgr::solve()")),
00275 _timerRestarting(Teuchos::TimeMonitor::getNewTimer("BlockDavidsonSolMgr restarting")),
00276 _timerLocking(Teuchos::TimeMonitor::getNewTimer("BlockDavidsonSolMgr locking"))
00277 {
00278 TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
00279 TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
00280 TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
00281 TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
00282
00283 std::string strtmp;
00284
00285
00286 whch_ = pl.get("Which",whch_);
00287 TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",std::invalid_argument, "Invalid sorting string.");
00288
00289
00290 ortho_ = pl.get("Orthogonalization",ortho_);
00291 if (ortho_ != "DGKS" && ortho_ != "SVQB") {
00292 ortho_ = "SVQB";
00293 }
00294
00295
00296 convtol_ = pl.get("Convergence Tolerance",convtol_);
00297 relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
00298 strtmp = pl.get("Convergence Norm",std::string("2"));
00299 if (strtmp == "2") {
00300 convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00301 }
00302 else if (strtmp == "M") {
00303 convNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00304 }
00305 else {
00306 TEST_FOR_EXCEPTION(true, std::invalid_argument,
00307 "Anasazi::BlockDavidsonSolMgr: Invalid Convergence Norm.");
00308 }
00309
00310
00311 useLocking_ = pl.get("Use Locking",useLocking_);
00312 rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
00313
00314 locktol_ = convtol_/10;
00315 locktol_ = pl.get("Locking Tolerance",locktol_);
00316 strtmp = pl.get("Locking Norm",std::string("2"));
00317 if (strtmp == "2") {
00318 lockNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_2NORM;
00319 }
00320 else if (strtmp == "M") {
00321 lockNorm_ = StatusTestResNorm<ScalarType,MV,OP>::RES_ORTH;
00322 }
00323 else {
00324 TEST_FOR_EXCEPTION(true, std::invalid_argument,
00325 "Anasazi::BlockDavidsonSolMgr: Invalid Locking Norm.");
00326 }
00327
00328
00329 maxRestarts_ = pl.get("Maximum Restarts",maxRestarts_);
00330
00331
00332 blockSize_ = pl.get("Block Size",problem_->getNEV());
00333 TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
00334 "Anasazi::BlockDavidsonSolMgr: \"Block Size\" must be strictly positive.");
00335 numBlocks_ = pl.get("Num Blocks",2);
00336 TEST_FOR_EXCEPTION(numBlocks_ <= 1, std::invalid_argument,
00337 "Anasazi::BlockDavidsonSolMgr: \"Num Blocks\" must be >= 1.");
00338
00339
00340 if (useLocking_) {
00341 maxLocked_ = pl.get("Max Locked",problem_->getNEV());
00342 }
00343 else {
00344 maxLocked_ = 0;
00345 }
00346 if (maxLocked_ == 0) {
00347 useLocking_ = false;
00348 }
00349 TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
00350 "Anasazi::BlockDavidsonSolMgr: \"Max Locked\" must be positive.");
00351 TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
00352 std::invalid_argument,
00353 "Anasazi::BlockDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
00354 TEST_FOR_EXCEPTION(numBlocks_*blockSize_ + maxLocked_ > MVT::GetVecLength(*problem_->getInitVec()),
00355 std::invalid_argument,
00356 "Anasazi::BlockDavidsonSolMgr: Potentially impossible orthogonality requests. Reduce basis size or locking size.");
00357
00358 if (useLocking_) {
00359 lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
00360 TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
00361 std::invalid_argument,
00362 "Anasazi::BlockDavidsonSolMgr: \"Locking Quorum\" must be strictly positive.");
00363 }
00364
00365
00366 numRestartBlocks_ = pl.get("Num Restart Blocks",numRestartBlocks_);
00367 TEST_FOR_EXCEPTION(numRestartBlocks_ <= 0, std::invalid_argument,
00368 "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly positive.");
00369 TEST_FOR_EXCEPTION(numRestartBlocks_ >= numBlocks_, std::invalid_argument,
00370 "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly less than \"Num Blocks\".");
00371
00372
00373 if (pl.isParameter("In Situ Restarting")) {
00374 if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
00375 inSituRestart_ = pl.get("In Situ Restarting",inSituRestart_);
00376 } else {
00377 inSituRestart_ = (bool)Teuchos::getParameter<int>(pl,"In Situ Restarting");
00378 }
00379 }
00380
00381
00382 std::string fntemplate = "";
00383 bool allProcs = false;
00384 if (pl.isParameter("Output on all processors")) {
00385 if (Teuchos::isParameterType<bool>(pl,"Output on all processors")) {
00386 allProcs = pl.get("Output on all processors",allProcs);
00387 } else {
00388 allProcs = (bool)Teuchos::getParameter<int>(pl,"Output on all processors");
00389 }
00390 }
00391 fntemplate = pl.get("Output filename template",fntemplate);
00392 int MyPID;
00393 # ifdef HAVE_MPI
00394
00395 int mpiStarted = 0;
00396 MPI_Initialized(&mpiStarted);
00397 if (mpiStarted) MPI_Comm_rank(MPI_COMM_WORLD, &MyPID);
00398 else MyPID=0;
00399 # else
00400 MyPID = 0;
00401 # endif
00402 if (fntemplate != "") {
00403 std::ostringstream MyPIDstr;
00404 MyPIDstr << MyPID;
00405
00406 int pos, start=0;
00407 while ( (pos = fntemplate.find("%d",start)) != -1 ) {
00408 fntemplate.replace(pos,2,MyPIDstr.str());
00409 start = pos+2;
00410 }
00411 }
00412 Teuchos::RCP<ostream> osp;
00413 if (fntemplate != "") {
00414 osp = Teuchos::rcp( new std::ofstream(fntemplate.c_str(),std::ios::out | std::ios::app) );
00415 if (!*osp) {
00416 osp = Teuchos::rcpFromRef(std::cout);
00417 std::cout << "Anasazi::BlockDavidsonSolMgr::constructor(): Could not open file for write: " << fntemplate << std::endl;
00418 }
00419 }
00420 else {
00421 osp = Teuchos::rcpFromRef(std::cout);
00422 }
00423
00424 int verbosity = Anasazi::Errors;
00425 if (pl.isParameter("Verbosity")) {
00426 if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
00427 verbosity = pl.get("Verbosity", verbosity);
00428 } else {
00429 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
00430 }
00431 }
00432 if (allProcs) {
00433
00434 printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,MyPID) );
00435 }
00436 else {
00437
00438 printer_ = Teuchos::rcp( new BasicOutputManager<ScalarType>(verbosity,osp,0) );
00439 }
00440 }
00441
00442
00443
00444 template<class ScalarType, class MV, class OP>
00445 ReturnType
00446 BlockDavidsonSolMgr<ScalarType,MV,OP>::solve() {
00447
00448 typedef SolverUtils<ScalarType,MV,OP> msutils;
00449
00450 const int nev = problem_->getNEV();
00451
00452 #ifdef TEUCHOS_DEBUG
00453 Teuchos::RCP<Teuchos::FancyOStream>
00454 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
00455 out->setShowAllFrontMatter(false).setShowProcRank(true);
00456 *out << "Entering Anasazi::BlockDavidsonSolMgr::solve()\n";
00457 #endif
00458
00460
00461 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
00462
00464
00465
00466
00467 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
00468 if (globalTest_ == Teuchos::null) {
00469 convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
00470 }
00471 else {
00472 convtest = globalTest_;
00473 }
00474 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
00475 = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
00476
00477 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
00478 if (useLocking_) {
00479 if (lockingTest_ == Teuchos::null) {
00480 locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
00481 }
00482 else {
00483 locktest = lockingTest_;
00484 }
00485 }
00486
00487 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
00488 alltests.push_back(ordertest);
00489 if (locktest != Teuchos::null) alltests.push_back(locktest);
00490 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
00491
00492 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
00493 = Teuchos::rcp( new StatusTestCombo<ScalarType,MV,OP>( StatusTestCombo<ScalarType,MV,OP>::OR, alltests) );
00494
00495 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
00496 if ( printer_->isVerbosity(Debug) ) {
00497 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
00498 }
00499 else {
00500 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
00501 }
00502
00504
00505 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
00506 if (ortho_=="SVQB") {
00507 ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00508 } else if (ortho_=="DGKS") {
00509 ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
00510 } else {
00511 TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid orthogonalization type.");
00512 }
00513
00515
00516 Teuchos::ParameterList plist;
00517 plist.set("Block Size",blockSize_);
00518 plist.set("Num Blocks",numBlocks_);
00519
00521
00522 Teuchos::RCP<BlockDavidson<ScalarType,MV,OP> > bd_solver
00523 = Teuchos::rcp( new BlockDavidson<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,plist) );
00524
00525 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
00526 if (probauxvecs != Teuchos::null) {
00527 bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
00528 }
00529
00531
00532
00533
00534 int curNumLocked = 0;
00535 Teuchos::RCP<MV> lockvecs;
00536
00537
00538
00539
00540
00541 if (maxLocked_ > 0) {
00542 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
00543 }
00544 std::vector<MagnitudeType> lockvals;
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592 Teuchos::RCP<MV> workMV;
00593 if (inSituRestart_ == false) {
00594
00595 if (useLocking_==true || maxRestarts_ > 0) {
00596 workMV = MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
00597 }
00598 else {
00599
00600 workMV = Teuchos::null;
00601 }
00602 }
00603 else {
00604
00605
00606
00607
00608
00609 if (maxRestarts_ > 0) {
00610 workMV = MVT::Clone(*problem_->getInitVec(),1);
00611 }
00612 else {
00613 workMV = Teuchos::null;
00614 }
00615 }
00616
00617
00618 const ScalarType ONE = SCT::one();
00619 const ScalarType ZERO = SCT::zero();
00620 Teuchos::LAPACK<int,ScalarType> lapack;
00621 Teuchos::BLAS<int,ScalarType> blas;
00622
00623
00624 Eigensolution<ScalarType,MV> sol;
00625 sol.numVecs = 0;
00626 problem_->setSolution(sol);
00627
00628 int numRestarts = 0;
00629
00630
00631 {
00632 Teuchos::TimeMonitor slvtimer(*_timerSolve);
00633
00634
00635 while (1) {
00636 try {
00637 bd_solver->iterate();
00638
00640
00641
00642
00644 if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
00645 throw AnasaziError("Anasazi::BlockDavidsonSolMgr::solve(): User-specified debug status test returned Passed.");
00646 }
00648
00649
00650
00652 else if (ordertest->getStatus() == Passed ) {
00653
00654
00655
00656 break;
00657 }
00659
00660
00661
00663 else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
00664
00665 Teuchos::TimeMonitor restimer(*_timerRestarting);
00666
00667 if ( numRestarts >= maxRestarts_ ) {
00668 break;
00669 }
00670 numRestarts++;
00671
00672 printer_->stream(IterationDetails) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
00673
00674 BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
00675 int curdim = state.curDim;
00676 int newdim = numRestartBlocks_*blockSize_;
00677
00678 # ifdef TEUCHOS_DEBUG
00679 {
00680 std::vector<Value<ScalarType> > ritzvalues = bd_solver->getRitzValues();
00681 *out << "Ritz values from solver:\n";
00682 for (unsigned int i=0; i<ritzvalues.size(); i++) {
00683 *out << ritzvalues[i].realpart << " ";
00684 }
00685 *out << "\n";
00686 }
00687 # endif
00688
00689
00690
00691 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
00692 std::vector<MagnitudeType> theta(curdim);
00693 int rank = curdim;
00694 # ifdef TEUCHOS_DEBUG
00695 {
00696 std::stringstream os;
00697 os << "KK before HEEV...\n"
00698 << *state.KK << "\n";
00699 *out << os.str();
00700 }
00701 # endif
00702 int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
00703 TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
00704 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");
00705 TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
00706 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors.");
00707
00708 # ifdef TEUCHOS_DEBUG
00709 {
00710 std::stringstream os;
00711 *out << "Ritz values from heev(KK):\n";
00712 for (unsigned int i=0; i<theta.size(); i++) *out << theta[i] << " ";
00713 os << "\nRitz vectors from heev(KK):\n"
00714 << S << "\n";
00715 *out << os.str();
00716 }
00717 # endif
00718
00719
00720
00721 {
00722 std::vector<int> order(curdim);
00723 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
00724
00725
00726 msutils::permuteVectors(order,S);
00727 }
00728 # ifdef TEUCHOS_DEBUG
00729 {
00730 std::stringstream os;
00731 *out << "Ritz values from heev(KK) after sorting:\n";
00732 std::copy(theta.begin(), theta.end(), std::ostream_iterator<ScalarType>(*out, " "));
00733 os << "\nRitz vectors from heev(KK) after sorting:\n"
00734 << S << "\n";
00735 *out << os.str();
00736 }
00737 # endif
00738
00739
00740 Teuchos::SerialDenseMatrix<int,ScalarType> Sr(Teuchos::View,S,curdim,newdim);
00741 # ifdef TEUCHOS_DEBUG
00742 {
00743 std::stringstream os;
00744 os << "Significant primitive Ritz vectors:\n"
00745 << Sr << "\n";
00746 *out << os.str();
00747 }
00748 # endif
00749
00750
00751 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(newdim,newdim);
00752 {
00753 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,newdim),
00754 KKold(Teuchos::View,*state.KK,curdim,curdim);
00755 int teuchosRet;
00756
00757 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Sr,ZERO);
00758 TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
00759 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
00760
00761 teuchosRet = newKK.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sr,KKtmp,ZERO);
00762 TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
00763 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
00764
00765 for (int j=0; j<newdim-1; ++j) {
00766 for (int i=j+1; i<newdim; ++i) {
00767 newKK(i,j) = SCT::conjugate(newKK(j,i));
00768 }
00769 }
00770 }
00771 # ifdef TEUCHOS_DEBUG
00772 {
00773 std::stringstream os;
00774 os << "Sr'*KK*Sr:\n"
00775 << newKK << "\n";
00776 *out << os.str();
00777 }
00778 # endif
00779
00780
00781 BlockDavidsonState<ScalarType,MV> rstate;
00782 rstate.curDim = newdim;
00783 rstate.KK = Teuchos::rcpFromRef(newKK);
00784
00785
00786
00787
00788
00789 rstate.X = state.X;
00790 rstate.KX = state.KX;
00791 rstate.MX = state.MX;
00792 rstate.R = state.R;
00793 rstate.T = Teuchos::rcp( new std::vector<MagnitudeType>(theta.begin(),theta.begin()+newdim) );
00794
00795 if (inSituRestart_ == true) {
00796 # ifdef TEUCHOS_DEBUG
00797 *out << "Beginning in-situ restart...\n";
00798 # endif
00799
00800
00801 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
00802
00803
00804
00805 std::vector<ScalarType> tau(newdim), work(newdim);
00806 int geqrf_info;
00807 lapack.GEQRF(curdim,newdim,Sr.values(),Sr.stride(),&tau[0],&work[0],work.size(),&geqrf_info);
00808 TEST_FOR_EXCEPTION(geqrf_info != 0,std::logic_error,
00809 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
00810 if (printer_->isVerbosity(Debug)) {
00811 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,Sr,newdim,newdim);
00812 for (int j=0; j<newdim; j++) {
00813 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
00814 for (int i=j+1; i<newdim; i++) {
00815 R(i,j) = ZERO;
00816 }
00817 }
00818 printer_->stream(Debug) << "||Triangular factor of Sr - I||: " << R.normFrobenius() << std::endl;
00819 }
00820
00821
00822
00823
00824 {
00825 std::vector<int> curind(curdim);
00826 for (int i=0; i<curdim; i++) curind[i] = i;
00827 Teuchos::RCP<MV> oldV = MVT::CloneView(*solverbasis,curind);
00828 msutils::applyHouse(newdim,*oldV,Sr,tau,workMV);
00829 }
00830
00831
00832
00833
00834 rstate.V = solverbasis;
00835 }
00836 else {
00837 # ifdef TEUCHOS_DEBUG
00838 *out << "Beginning ex-situ restart...\n";
00839 # endif
00840
00841 std::vector<int> curind(curdim), newind(newdim);
00842 for (int i=0; i<curdim; i++) curind[i] = i;
00843 for (int i=0; i<newdim; i++) newind[i] = i;
00844 Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
00845 Teuchos::RCP<MV> newV = MVT::CloneView(*workMV ,newind);
00846
00847 MVT::MvTimesMatAddMv(ONE,*oldV,Sr,ZERO,*newV);
00848
00849
00850 rstate.V = newV;
00851 }
00852
00853
00854
00855 bd_solver->initialize(rstate);
00856 }
00858
00859
00860
00862 else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
00863
00864 Teuchos::TimeMonitor lcktimer(*_timerLocking);
00865
00866
00867
00868 BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
00869 const int curdim = state.curDim;
00870
00871
00872
00873 TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
00874 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() non-positive.");
00875 TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
00876 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
00877
00878 TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
00879 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: locking not deactivated.");
00880
00881
00882 std::vector<int> tmp_vector_int;
00883 if (curNumLocked + locktest->howMany() > maxLocked_) {
00884
00885 tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
00886 }
00887 else {
00888 tmp_vector_int = locktest->whichVecs();
00889 }
00890 const std::vector<int> lockind(tmp_vector_int);
00891 const int numNewLocked = lockind.size();
00892
00893
00894
00895 const int numUnlocked = curdim-numNewLocked;
00896 tmp_vector_int.resize(curdim);
00897 for (int i=0; i<curdim; i++) tmp_vector_int[i] = i;
00898 const std::vector<int> curind(tmp_vector_int);
00899 tmp_vector_int.resize(numUnlocked);
00900 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
00901 const std::vector<int> unlockind(tmp_vector_int);
00902 tmp_vector_int.clear();
00903
00904
00905
00906 if (printer_->isVerbosity(Debug)) {
00907 printer_->print(Debug,"Locking vectors: ");
00908 for (unsigned int i=0; i<lockind.size(); i++) {printer_->stream(Debug) << " " << lockind[i];}
00909 printer_->print(Debug,"\n");
00910 printer_->print(Debug,"Not locking vectors: ");
00911 for (unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(Debug) << " " << unlockind[i];}
00912 printer_->print(Debug,"\n");
00913 }
00914
00915
00916
00917
00918
00919
00920
00921
00922
00923 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
00924 std::vector<MagnitudeType> theta(curdim);
00925 {
00926 int rank = curdim;
00927 int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
00928 TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
00929 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver.");
00930 TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
00931 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors.");
00932
00933
00934 std::vector<int> order(curdim);
00935 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
00936
00937
00938 msutils::permuteVectors(order,S);
00939 }
00940
00941
00942
00943
00944 Teuchos::SerialDenseMatrix<int,ScalarType> Su(curdim,numUnlocked);
00945 for (int i=0; i<numUnlocked; i++) {
00946 blas.COPY(curdim, S[unlockind[i]], 1, Su[i], 1);
00947 }
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
00959 Teuchos::RCP<MV> defV, augV;
00960 if (inSituRestart_ == true) {
00961
00962
00963 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
00964
00965
00966
00967 Teuchos::SerialDenseMatrix<int,ScalarType> copySu(Su);
00968 std::vector<ScalarType> tau(numUnlocked), work(numUnlocked);
00969 int info;
00970 lapack.GEQRF(curdim,numUnlocked,copySu.values(),copySu.stride(),&tau[0],&work[0],work.size(),&info);
00971 TEST_FOR_EXCEPTION(info != 0,std::logic_error,
00972 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
00973 if (printer_->isVerbosity(Debug)) {
00974 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copySu,numUnlocked,numUnlocked);
00975 for (int j=0; j<numUnlocked; j++) {
00976 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
00977 for (int i=j+1; i<numUnlocked; i++) {
00978 R(i,j) = ZERO;
00979 }
00980 }
00981 printer_->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
00982 }
00983
00984
00985
00986
00987 {
00988 Teuchos::RCP<MV> oldV = MVT::CloneView(*solverbasis,curind);
00989 msutils::applyHouse(numUnlocked,*oldV,copySu,tau,workMV);
00990 }
00991 std::vector<int> defind(numUnlocked), augind(numNewLocked);
00992 for (int i=0; i<numUnlocked ; i++) defind[i] = i;
00993 for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
00994 defV = MVT::CloneView(*solverbasis,defind);
00995 augV = MVT::CloneView(*solverbasis,augind);
00996 }
00997 else {
00998
00999 std::vector<int> defind(numUnlocked), augind(numNewLocked);
01000 for (int i=0; i<numUnlocked ; i++) defind[i] = i;
01001 for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
01002 Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
01003 defV = MVT::CloneView(*workMV,defind);
01004 augV = MVT::CloneView(*workMV,augind);
01005
01006 MVT::MvTimesMatAddMv(ONE,*oldV,Su,ZERO,*defV);
01007 }
01008
01009
01010
01011
01012
01013
01014
01015
01016
01017
01018
01019
01020 Teuchos::RCP<const MV> curlocked, newLocked;
01021 Teuchos::RCP<MV> augTmp;
01022 {
01023
01024 if (curNumLocked > 0) {
01025 std::vector<int> curlockind(curNumLocked);
01026 for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
01027 curlocked = MVT::CloneView(*lockvecs,curlockind);
01028 }
01029 else {
01030 curlocked = Teuchos::null;
01031 }
01032
01033 std::vector<int> augtmpind(numNewLocked);
01034 for (int i=0; i<numNewLocked; i++) augtmpind[i] = curNumLocked+i;
01035 augTmp = MVT::CloneView(*lockvecs,augtmpind);
01036
01037 newLocked = MVT::CloneView(*bd_solver->getRitzVectors(),lockind);
01038 }
01039
01040
01041
01042
01043 MVT::MvRandom(*augV);
01044
01045
01046
01047 {
01048 Teuchos::Array<Teuchos::RCP<const MV> > against;
01049 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
01050 if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
01051 if (curlocked != Teuchos::null) against.push_back(curlocked);
01052 against.push_back(newLocked);
01053 against.push_back(defV);
01054 if (problem_->getM() != Teuchos::null) {
01055 OPT::Apply(*problem_->getM(),*augV,*augTmp);
01056 }
01057 ortho->projectAndNormalizeMat(*augV,against,dummyC,Teuchos::null,augTmp);
01058 }
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(curdim,curdim);
01069 {
01070 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,numUnlocked),
01071 KKold(Teuchos::View,*state.KK,curdim,curdim),
01072 KK11(Teuchos::View,newKK,numUnlocked,numUnlocked);
01073 int teuchosRet;
01074
01075 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Su,ZERO);
01076 TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
01077 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
01078
01079 teuchosRet = KK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Su,KKtmp,ZERO);
01080 TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
01081 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
01082 }
01083
01084
01085 {
01086 OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
01087 Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,newKK,numUnlocked,numNewLocked,0,numUnlocked),
01088 KK22(Teuchos::View,newKK,numNewLocked,numNewLocked,numUnlocked,numUnlocked);
01089 MVT::MvTransMv(ONE,*defV,*augTmp,KK12);
01090 MVT::MvTransMv(ONE,*augV,*augTmp,KK22);
01091 }
01092
01093
01094 defV = Teuchos::null;
01095 augV = Teuchos::null;
01096
01097
01098 for (int j=0; j<curdim; ++j) {
01099 for (int i=j+1; i<curdim; ++i) {
01100 newKK(i,j) = SCT::conjugate(newKK(j,i));
01101 }
01102 }
01103
01104
01105
01106 augTmp = Teuchos::null;
01107 {
01108 std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
01109 for (int i=0; i<numNewLocked; i++) {
01110 lockvals.push_back(allvals[lockind[i]].realpart);
01111 }
01112
01113 std::vector<int> indlock(numNewLocked);
01114 for (int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
01115 MVT::SetBlock(*newLocked,indlock,*lockvecs);
01116 newLocked = Teuchos::null;
01117
01118 curNumLocked += numNewLocked;
01119 std::vector<int> curlockind(curNumLocked);
01120 for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
01121 curlocked = MVT::CloneView(*lockvecs,curlockind);
01122 }
01123
01124
01125
01126 {
01127 ordertest->setAuxVals(lockvals);
01128
01129 Teuchos::Array< Teuchos::RCP<const MV> > aux;
01130 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
01131 aux.push_back(curlocked);
01132 bd_solver->setAuxVecs(aux);
01133
01134 if (curNumLocked == maxLocked_) {
01135
01136 combotest->removeTest(locktest);
01137 }
01138 }
01139
01140
01141
01142 BlockDavidsonState<ScalarType,MV> rstate;
01143 rstate.curDim = curdim;
01144 if (inSituRestart_) {
01145
01146 rstate.V = state.V;
01147 }
01148 else {
01149
01150 rstate.V = workMV;
01151 }
01152 rstate.KK = Teuchos::rcpFromRef(newKK);
01153
01154
01155 bd_solver->initialize(rstate);
01156 }
01158
01159
01160
01161
01163 else {
01164 TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
01165 }
01166 }
01167 catch (const AnasaziError &err) {
01168 printer_->stream(Errors)
01169 << "Anasazi::BlockDavidsonSolMgr::solve() caught unexpected exception from Anasazi::BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << std::endl
01170 << err.what() << std::endl
01171 << "Anasazi::BlockDavidsonSolMgr::solve() returning Unconverged with no solutions." << std::endl;
01172 return Unconverged;
01173 }
01174 }
01175
01176
01177 workMV = Teuchos::null;
01178
01179 sol.numVecs = ordertest->howMany();
01180 if (sol.numVecs > 0) {
01181 sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
01182 sol.Espace = sol.Evecs;
01183 sol.Evals.resize(sol.numVecs);
01184 std::vector<MagnitudeType> vals(sol.numVecs);
01185
01186
01187 std::vector<int> which = ordertest->whichVecs();
01188
01189
01190
01191 std::vector<int> inlocked(0), insolver(0);
01192 for (unsigned int i=0; i<which.size(); i++) {
01193 if (which[i] >= 0) {
01194 TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): positive indexing mistake from ordertest.");
01195 insolver.push_back(which[i]);
01196 }
01197 else {
01198
01199 TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): negative indexing mistake from ordertest.");
01200 inlocked.push_back(which[i] + curNumLocked);
01201 }
01202 }
01203
01204 TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
01205
01206
01207 if (insolver.size() > 0) {
01208
01209 int lclnum = insolver.size();
01210 std::vector<int> tosol(lclnum);
01211 for (int i=0; i<lclnum; i++) tosol[i] = i;
01212 Teuchos::RCP<const MV> v = MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
01213 MVT::SetBlock(*v,tosol,*sol.Evecs);
01214
01215 std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
01216 for (unsigned int i=0; i<insolver.size(); i++) {
01217 vals[i] = fromsolver[insolver[i]].realpart;
01218 }
01219 }
01220
01221
01222 if (inlocked.size() > 0) {
01223 int solnum = insolver.size();
01224
01225 int lclnum = inlocked.size();
01226 std::vector<int> tosol(lclnum);
01227 for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
01228 Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
01229 MVT::SetBlock(*v,tosol,*sol.Evecs);
01230
01231 for (unsigned int i=0; i<inlocked.size(); i++) {
01232 vals[i+solnum] = lockvals[inlocked[i]];
01233 }
01234 }
01235
01236
01237 {
01238 std::vector<int> order(sol.numVecs);
01239 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
01240
01241 for (int i=0; i<sol.numVecs; i++) {
01242 sol.Evals[i].realpart = vals[i];
01243 sol.Evals[i].imagpart = MT::zero();
01244 }
01245
01246 msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
01247 }
01248
01249
01250 sol.index.resize(sol.numVecs,0);
01251 }
01252 }
01253
01254
01255 bd_solver->currentStatus(printer_->stream(FinalSummary));
01256
01257
01258 Teuchos::TimeMonitor::summarize(printer_->stream(TimingDetails));
01259
01260 problem_->setSolution(sol);
01261 printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
01262
01263
01264 numIters_ = bd_solver->getNumIters();
01265
01266 if (sol.numVecs < nev) {
01267 return Unconverged;
01268 }
01269 return Converged;
01270 }
01271
01272
01273 template <class ScalarType, class MV, class OP>
01274 void
01275 BlockDavidsonSolMgr<ScalarType,MV,OP>::setGlobalStatusTest(
01276 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
01277 {
01278 globalTest_ = global;
01279 }
01280
01281 template <class ScalarType, class MV, class OP>
01282 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
01283 BlockDavidsonSolMgr<ScalarType,MV,OP>::getGlobalStatusTest() const
01284 {
01285 return globalTest_;
01286 }
01287
01288 template <class ScalarType, class MV, class OP>
01289 void
01290 BlockDavidsonSolMgr<ScalarType,MV,OP>::setDebugStatusTest(
01291 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
01292 {
01293 debugTest_ = debug;
01294 }
01295
01296 template <class ScalarType, class MV, class OP>
01297 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
01298 BlockDavidsonSolMgr<ScalarType,MV,OP>::getDebugStatusTest() const
01299 {
01300 return debugTest_;
01301 }
01302
01303 template <class ScalarType, class MV, class OP>
01304 void
01305 BlockDavidsonSolMgr<ScalarType,MV,OP>::setLockingStatusTest(
01306 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
01307 {
01308 lockingTest_ = locking;
01309 }
01310
01311 template <class ScalarType, class MV, class OP>
01312 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
01313 BlockDavidsonSolMgr<ScalarType,MV,OP>::getLockingStatusTest() const
01314 {
01315 return lockingTest_;
01316 }
01317
01318 }
01319
01320 #endif