Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziBlockDavidsonSolMgr.hpp
Go to the documentation of this file.
1// @HEADER
2// ***********************************************************************
3//
4// Anasazi: Block Eigensolvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ***********************************************************************
40// @HEADER
41
42#ifndef ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
43#define ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP
44
49#include "AnasaziConfigDefs.hpp"
50#include "AnasaziTypes.hpp"
51
55
57#include "AnasaziBasicSort.hpp"
66#include "Teuchos_BLAS.hpp"
67#include "Teuchos_LAPACK.hpp"
68#include "Teuchos_TimeMonitor.hpp"
69#include "Teuchos_FancyOStream.hpp"
70
71
85namespace Anasazi {
86
87
120template<class ScalarType, class MV, class OP>
121class BlockDavidsonSolMgr : public SolverManager<ScalarType,MV,OP> {
122
123 private:
126 typedef Teuchos::ScalarTraits<ScalarType> SCT;
127 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
128 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
129
130 public:
131
133
134
160 BlockDavidsonSolMgr( const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
161 Teuchos::ParameterList &pl );
162
166
168
169
172 return *problem_;
173 }
174
176 int getNumIters() const {
177 return numIters_;
178 }
179
187 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
188 return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
189 }
190
192
194
195
218
220 void setGlobalStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global);
221
223 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
224
226 void setLockingStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking);
227
229 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
230
232 void setDebugStatusTest(const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug);
233
235 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
236
238
239 private:
240 Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
241
242 std::string whch_, ortho_;
243
244 MagnitudeType convtol_, locktol_;
245 int maxRestarts_;
246 bool useLocking_;
247 bool relconvtol_, rellocktol_;
248 int blockSize_, numBlocks_, numIters_;
249 int maxLocked_;
250 int lockQuorum_;
251 bool inSituRestart_;
252 int numRestartBlocks_;
253 enum ResType convNorm_, lockNorm_;
254
255 Teuchos::RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
256
257 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
258 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
259 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
260
261 Teuchos::RCP<OutputManager<ScalarType> > printer_;
262};
263
264
265// Constructor
266template<class ScalarType, class MV, class OP>
268 const Teuchos::RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
269 Teuchos::ParameterList &pl ) :
270 problem_(problem),
271 whch_("SR"),
272 ortho_("SVQB"),
273 convtol_(MT::prec()),
274 maxRestarts_(20),
275 useLocking_(false),
276 relconvtol_(true),
277 rellocktol_(true),
278 blockSize_(0),
279 numBlocks_(0),
280 numIters_(0),
281 maxLocked_(0),
282 lockQuorum_(1),
283 inSituRestart_(false),
284 numRestartBlocks_(1)
285#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
286 , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr::solve()")),
287 _timerRestarting(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr restarting")),
288 _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: BlockDavidsonSolMgr locking"))
289#endif
290{
291 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
292 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
293 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
294 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
295
296 std::string strtmp;
297
298 // which values to solve for
299 whch_ = pl.get("Which",whch_);
300 TEUCHOS_TEST_FOR_EXCEPTION(whch_ != "SM" && whch_ != "LM" && whch_ != "SR" && whch_ != "LR",std::invalid_argument, "Invalid sorting string.");
301
302 // which orthogonalization to use
303 ortho_ = pl.get("Orthogonalization",ortho_);
304 if (ortho_ != "DGKS" && ortho_ != "SVQB") {
305 ortho_ = "SVQB";
306 }
307
308 // convergence tolerance
309 convtol_ = pl.get("Convergence Tolerance",convtol_);
310 relconvtol_ = pl.get("Relative Convergence Tolerance",relconvtol_);
311 strtmp = pl.get("Convergence Norm",std::string("2"));
312 if (strtmp == "2") {
313 convNorm_ = RES_2NORM;
314 }
315 else if (strtmp == "M") {
316 convNorm_ = RES_ORTH;
317 }
318 else {
319 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
320 "Anasazi::BlockDavidsonSolMgr: Invalid Convergence Norm.");
321 }
322
323 // locking tolerance
324 useLocking_ = pl.get("Use Locking",useLocking_);
325 rellocktol_ = pl.get("Relative Locking Tolerance",rellocktol_);
326 // default: should be less than convtol_
327 locktol_ = convtol_/10;
328 locktol_ = pl.get("Locking Tolerance",locktol_);
329 strtmp = pl.get("Locking Norm",std::string("2"));
330 if (strtmp == "2") {
331 lockNorm_ = RES_2NORM;
332 }
333 else if (strtmp == "M") {
334 lockNorm_ = RES_ORTH;
335 }
336 else {
337 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
338 "Anasazi::BlockDavidsonSolMgr: Invalid Locking Norm.");
339 }
340
341 // maximum number of restarts
342 maxRestarts_ = pl.get("Maximum Restarts",maxRestarts_);
343
344 // block size: default is nev()
345 blockSize_ = pl.get("Block Size",problem_->getNEV());
346 TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
347 "Anasazi::BlockDavidsonSolMgr: \"Block Size\" must be strictly positive.");
348 numBlocks_ = pl.get("Num Blocks",2);
349 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 1, std::invalid_argument,
350 "Anasazi::BlockDavidsonSolMgr: \"Num Blocks\" must be >= 1.");
351
352 // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
353 if (useLocking_) {
354 maxLocked_ = pl.get("Max Locked",problem_->getNEV());
355 }
356 else {
357 maxLocked_ = 0;
358 }
359 if (maxLocked_ == 0) {
360 useLocking_ = false;
361 }
362 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ < 0, std::invalid_argument,
363 "Anasazi::BlockDavidsonSolMgr: \"Max Locked\" must be positive.");
364 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ + blockSize_ < problem_->getNEV(),
365 std::invalid_argument,
366 "Anasazi::BlockDavidsonSolMgr: Not enough storage space for requested number of eigenpairs.");
367 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_)*blockSize_ + maxLocked_ > MVT::GetGlobalLength(*problem_->getInitVec()),
368 std::invalid_argument,
369 "Anasazi::BlockDavidsonSolMgr: Potentially impossible orthogonality requests. Reduce basis size or locking size.");
370
371 if (useLocking_) {
372 lockQuorum_ = pl.get("Locking Quorum",lockQuorum_);
373 TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0,
374 std::invalid_argument,
375 "Anasazi::BlockDavidsonSolMgr: \"Locking Quorum\" must be strictly positive.");
376 }
377
378 // restart size
379 numRestartBlocks_ = pl.get("Num Restart Blocks",numRestartBlocks_);
380 TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ <= 0, std::invalid_argument,
381 "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly positive.");
382 TEUCHOS_TEST_FOR_EXCEPTION(numRestartBlocks_ >= numBlocks_, std::invalid_argument,
383 "Anasazi::BlockDavidsonSolMgr: \"Num Restart Blocks\" must be strictly less than \"Num Blocks\".");
384
385 // restarting technique: V*Q or applyHouse(V,H,tau)
386 if (pl.isParameter("In Situ Restarting")) {
387 if (Teuchos::isParameterType<bool>(pl,"In Situ Restarting")) {
388 inSituRestart_ = pl.get("In Situ Restarting",inSituRestart_);
389 } else {
390 inSituRestart_ = ( Teuchos::getParameter<int>(pl,"In Situ Restarting") != 0 );
391 }
392 }
393
394 // Output manager
395 // Create a formatted output stream to print to.
396 // See if user requests output processor.
397 int osProc = pl.get("Output Processor", 0);
398
399 // If not passed in by user, it will be chosen based upon operator type.
400 Teuchos::RCP<Teuchos::FancyOStream> osp;
401
402 // output stream
403 if (pl.isParameter("Output Stream")) {
404 osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
405 }
406 else {
407 osp = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc);
408 }
409
410 // verbosity
411 int verbosity = Anasazi::Errors;
412 if (pl.isParameter("Verbosity")) {
413 if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
414 verbosity = pl.get("Verbosity", verbosity);
415 } else {
416 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
417 }
418 }
419 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity,osp) );
420
421}
422
423
424// solve()
425template<class ScalarType, class MV, class OP>
428
429 typedef SolverUtils<ScalarType,MV,OP> msutils;
430
431 const int nev = problem_->getNEV();
432
433#ifdef TEUCHOS_DEBUG
434 Teuchos::RCP<Teuchos::FancyOStream>
435 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
436 out->setShowAllFrontMatter(false).setShowProcRank(true);
437 *out << "Entering Anasazi::BlockDavidsonSolMgr::solve()\n";
438#endif
439
441 // Sort manager
442 Teuchos::RCP<BasicSort<MagnitudeType> > sorter = Teuchos::rcp( new BasicSort<MagnitudeType>(whch_) );
443
445 // Status tests
446 //
447 // convergence
448 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convtest;
449 if (globalTest_ == Teuchos::null) {
450 convtest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(convtol_,nev,convNorm_,relconvtol_) );
451 }
452 else {
453 convtest = globalTest_;
454 }
455 Teuchos::RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
456 = Teuchos::rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
457 // locking
458 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > locktest;
459 if (useLocking_) {
460 if (lockingTest_ == Teuchos::null) {
461 locktest = Teuchos::rcp( new StatusTestResNorm<ScalarType,MV,OP>(locktol_,lockQuorum_,lockNorm_,rellocktol_) );
462 }
463 else {
464 locktest = lockingTest_;
465 }
466 }
467 // for a non-short-circuited OR test, the order doesn't matter
468 Teuchos::Array<Teuchos::RCP<StatusTest<ScalarType,MV,OP> > > alltests;
469 alltests.push_back(ordertest);
470 if (locktest != Teuchos::null) alltests.push_back(locktest);
471 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
472
473 Teuchos::RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
475 // printing StatusTest
476 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
477 if ( printer_->isVerbosity(Debug) ) {
478 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
479 }
480 else {
481 outputtest = Teuchos::rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
482 }
483
485 // Orthomanager
486 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
487 if (ortho_=="SVQB") {
488 ortho = Teuchos::rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
489 } else if (ortho_=="DGKS") {
490 ortho = Teuchos::rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
491 } else {
492 TEUCHOS_TEST_FOR_EXCEPTION(ortho_!="SVQB"&&ortho_!="DGKS",std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid orthogonalization type.");
493 }
494
496 // Parameter list
497 Teuchos::ParameterList plist;
498 plist.set("Block Size",blockSize_);
499 plist.set("Num Blocks",numBlocks_);
500
502 // BlockDavidson solver
503 Teuchos::RCP<BlockDavidson<ScalarType,MV,OP> > bd_solver
504 = Teuchos::rcp( new BlockDavidson<ScalarType,MV,OP>(problem_,sorter,printer_,outputtest,ortho,plist) );
505 // set any auxiliary vectors defined in the problem
506 Teuchos::RCP< const MV > probauxvecs = problem_->getAuxVecs();
507 if (probauxvecs != Teuchos::null) {
508 bd_solver->setAuxVecs( Teuchos::tuple< Teuchos::RCP<const MV> >(probauxvecs) );
509 }
510
512 // Storage
513 //
514 // lockvecs will contain eigenvectors that have been determined "locked" by the status test
515 int curNumLocked = 0;
516 Teuchos::RCP<MV> lockvecs;
517 // lockvecs is used to hold the locked eigenvectors, as well as for temporary storage when locking.
518 // when locking, we will lock some number of vectors numnew, where numnew <= maxlocked - curlocked
519 // we will produce numnew random vectors, which will go into the space with the new basis.
520 // we will also need numnew storage for the image of these random vectors under A and M;
521 // columns [curlocked+1,curlocked+numnew] will be used for this storage
522 if (maxLocked_ > 0) {
523 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
524 }
525 std::vector<MagnitudeType> lockvals;
526 //
527 // Restarting occurs under two scenarios: when the basis is full and after locking.
528 //
529 // For the former, a new basis of size blockSize*numRestartBlocks is generated using the current basis
530 // and the most significant primitive Ritz vectors (projected eigenvectors).
531 // [S,L] = eig(KK)
532 // S = [Sr St] // some for "r"estarting, some are "t"runcated
533 // newV = V*Sr
534 // KK_new = newV'*K*newV = Sr'*V'*K*V*Sr = Sr'*KK*Sr
535 // Therefore, the only multivector operation needed is for the generation of newV.
536 //
537 // * If the multiplication is explicit, it requires a workspace of blockSize*numRestartBlocks vectors.
538 // This space must be specifically allocated for that task, as we don't have any space of that size.
539 // It (workMV) will be allocated at the beginning of solve()
540 // * Optionally, the multiplication can be performed implicitly, via a Householder QR factorization of
541 // Sr. This can be done in situ, using the basis multivector contained in the solver. This requires
542 // that we cast away the const on the multivector returned from getState(). Workspace for this approach
543 // is a single vector. the solver's internal storage must be preserved (X,MX,KX,R), requiring us to
544 // allocate this vector.
545 //
546 // For the latter (restarting after locking), the new basis is the same size as existing basis. If numnew
547 // vectors are locked, they are deflated from the current basis and replaced with randomly generated
548 // vectors.
549 // [S,L] = eig(KK)
550 // S = [Sl Su] // partitioned: "l"ocked and "u"nlocked
551 // newL = V*Sl = X(locked)
552 // defV = V*Su
553 // augV = rand(numnew) // orthogonal to oldL,newL,defV,auxvecs
554 // newV = [defV augV]
555 // Kknew = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
556 // [augV'*K*defV augV'*K*augV]
557 // locked = [oldL newL]
558 // Clearly, this operation is more complicated than the previous.
559 // Here is a list of the significant computations that need to be performed:
560 // - newL will be put into space in lockvecs, but will be copied from getState().X at the end
561 // - defV,augV will be stored in workspace the size of the current basis.
562 // - If inSituRestart==true, we compute defV in situ in bd_solver::V_ and
563 // put augV at the end of bd_solver::V_
564 // - If inSituRestart==false, we must have curDim vectors available for
565 // defV and augV; we will allocate a multivector (workMV) at the beginning of solve()
566 // for this purpose.
567 // - M*augV and K*augV are needed; they will be stored in lockvecs. As a result, newL will
568 // not be put into lockvecs until the end.
569 //
570 // Therefore, we must allocate workMV when ((maxRestarts_ > 0) || (useLocking_ == true)) && inSituRestart == false
571 // It will be allocated to size (numBlocks-1)*blockSize
572 //
573 Teuchos::RCP<MV> workMV;
574 if (inSituRestart_ == false) {
575 // we need storage space to restart, either if we may lock or if may restart after a full basis
576 if (useLocking_==true || maxRestarts_ > 0) {
577 workMV = MVT::Clone(*problem_->getInitVec(),(numBlocks_-1)*blockSize_);
578 }
579 else {
580 // we will never need to restart.
581 workMV = Teuchos::null;
582 }
583 }
584 else { // inSituRestart_ == true
585 // we will restart in situ, if we need to restart
586 // three situation remain:
587 // - never restart => no space needed
588 // - only restart for locking (i.e., never restart full) => no space needed
589 // - restart for full basis => need one vector
590 if (maxRestarts_ > 0) {
591 workMV = MVT::Clone(*problem_->getInitVec(),1);
592 }
593 else {
594 workMV = Teuchos::null;
595 }
596 }
597
598 // some consts and utils
599 const ScalarType ONE = SCT::one();
600 const ScalarType ZERO = SCT::zero();
601 Teuchos::LAPACK<int,ScalarType> lapack;
602 Teuchos::BLAS<int,ScalarType> blas;
603
604 // go ahead and initialize the solution to nothing in case we throw an exception
606 sol.numVecs = 0;
607 problem_->setSolution(sol);
608
609 int numRestarts = 0;
610
611 // enter solve() iterations
612 {
613#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
614 Teuchos::TimeMonitor slvtimer(*_timerSolve);
615#endif
616
617 // tell bd_solver to iterate
618 while (1) {
619 try {
620 bd_solver->iterate();
621
623 //
624 // check user-specified debug test; if it passed, return an exception
625 //
627 if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
628 throw AnasaziError("Anasazi::BlockDavidsonSolMgr::solve(): User-specified debug status test returned Passed.");
629 }
631 //
632 // check convergence next
633 //
635 else if (ordertest->getStatus() == Passed ) {
636 // we have convergence
637 // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
638 // ordertest->howMany() will tell us how many
639 break;
640 }
642 //
643 // check for restarting before locking: if we need to lock, it will happen after the restart
644 //
646 else if ( bd_solver->getCurSubspaceDim() == bd_solver->getMaxSubspaceDim() ) {
647
648#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
649 Teuchos::TimeMonitor restimer(*_timerRestarting);
650#endif
651
652 if ( numRestarts >= maxRestarts_ ) {
653 break; // break from while(1){bd_solver->iterate()}
654 }
655 numRestarts++;
656
657 printer_->stream(IterationDetails) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
658
659 BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
660 int curdim = state.curDim;
661 int newdim = numRestartBlocks_*blockSize_;
662
663# ifdef TEUCHOS_DEBUG
664 {
665 std::vector<Value<ScalarType> > ritzvalues = bd_solver->getRitzValues();
666 *out << "Ritz values from solver:\n";
667 for (unsigned int i=0; i<ritzvalues.size(); i++) {
668 *out << ritzvalues[i].realpart << " ";
669 }
670 *out << "\n";
671 }
672# endif
673
674 //
675 // compute eigenvectors of the projected problem
676 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
677 std::vector<MagnitudeType> theta(curdim);
678 int rank = curdim;
679# ifdef TEUCHOS_DEBUG
680 {
681 std::stringstream os;
682 os << "KK before HEEV...\n"
683 << printMat(*state.KK) << "\n";
684 *out << os.str();
685 }
686# endif
687 int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
688 TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
689 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver."); // this should never happen
690 TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
691 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
692
693# ifdef TEUCHOS_DEBUG
694 {
695 std::stringstream os;
696 *out << "Ritz values from heev(KK):\n";
697 for (unsigned int i=0; i<theta.size(); i++) *out << theta[i] << " ";
698 os << "\nRitz vectors from heev(KK):\n"
699 << printMat(S) << "\n";
700 *out << os.str();
701 }
702# endif
703
704 //
705 // sort the eigenvalues (so that we can order the eigenvectors)
706 {
707 std::vector<int> order(curdim);
708 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
709 //
710 // apply the same ordering to the primitive ritz vectors
711 msutils::permuteVectors(order,S);
712 }
713# ifdef TEUCHOS_DEBUG
714 {
715 std::stringstream os;
716 *out << "Ritz values from heev(KK) after sorting:\n";
717 std::copy(theta.begin(), theta.end(), std::ostream_iterator<ScalarType>(*out, " "));
718 os << "\nRitz vectors from heev(KK) after sorting:\n"
719 << printMat(S) << "\n";
720 *out << os.str();
721 }
722# endif
723 //
724 // select the significant primitive ritz vectors
725 Teuchos::SerialDenseMatrix<int,ScalarType> Sr(Teuchos::View,S,curdim,newdim);
726# ifdef TEUCHOS_DEBUG
727 {
728 std::stringstream os;
729 os << "Significant primitive Ritz vectors:\n"
730 << printMat(Sr) << "\n";
731 *out << os.str();
732 }
733# endif
734 //
735 // generate newKK = Sr'*KKold*Sr
736 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(newdim,newdim);
737 {
738 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,newdim),
739 KKold(Teuchos::View,*state.KK,curdim,curdim);
740 int teuchosRet;
741 // KKtmp = KKold*Sr
742 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Sr,ZERO);
743 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
744 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
745 // newKK = Sr'*KKtmp = Sr'*KKold*Sr
746 teuchosRet = newKK.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Sr,KKtmp,ZERO);
747 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
748 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
749 // make it Hermitian in memory
750 for (int j=0; j<newdim-1; ++j) {
751 for (int i=j+1; i<newdim; ++i) {
752 newKK(i,j) = SCT::conjugate(newKK(j,i));
753 }
754 }
755 }
756# ifdef TEUCHOS_DEBUG
757 {
758 std::stringstream os;
759 os << "Sr'*KK*Sr:\n"
760 << printMat(newKK) << "\n";
761 *out << os.str();
762 }
763# endif
764
765 // prepare new state
767 rstate.curDim = newdim;
768 rstate.KK = Teuchos::rcpFromRef(newKK);
769 //
770 // we know that newX = newV*Sr(:,1:bS) = oldV*S(:1:bS) = oldX
771 // the restarting preserves the Ritz vectors and residual
772 // for the Ritz values, we want all of the values associated with newV.
773 // these have already been placed at the beginning of theta
774 rstate.X = state.X;
775 rstate.KX = state.KX;
776 rstate.MX = state.MX;
777 rstate.R = state.R;
778 rstate.T = Teuchos::rcp( new std::vector<MagnitudeType>(theta.begin(),theta.begin()+newdim) );
779
780 if (inSituRestart_ == true) {
781# ifdef TEUCHOS_DEBUG
782 *out << "Beginning in-situ restart...\n";
783# endif
784 //
785 // get non-const pointer to solver's basis so we can work in situ
786 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
787 //
788 // perform Householder QR of Sr = Q [D;0], where D is unit diag.
789 // WARNING: this will overwrite Sr; however, we do not need Sr anymore after this
790 std::vector<ScalarType> tau(newdim), work(newdim);
791 int geqrf_info;
792 lapack.GEQRF(curdim,newdim,Sr.values(),Sr.stride(),&tau[0],&work[0],work.size(),&geqrf_info);
793 TEUCHOS_TEST_FOR_EXCEPTION(geqrf_info != 0,std::logic_error,
794 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
795 if (printer_->isVerbosity(Debug)) {
796 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,Sr,newdim,newdim);
797 for (int j=0; j<newdim; j++) {
798 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
799 for (int i=j+1; i<newdim; i++) {
800 R(i,j) = ZERO;
801 }
802 }
803 printer_->stream(Debug) << "||Triangular factor of Sr - I||: " << R.normFrobenius() << std::endl;
804 }
805 //
806 // perform implicit oldV*Sr
807 // this actually performs oldV*[Sr Su*M] = [newV truncV], for some unitary M
808 // we are actually interested in only the first newdim vectors of the result
809 {
810 std::vector<int> curind(curdim);
811 for (int i=0; i<curdim; i++) curind[i] = i;
812 Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
813 msutils::applyHouse(newdim,*oldV,Sr,tau,workMV);
814 }
815 //
816 // put the new basis into the state for initialize()
817 // the new basis is contained in the the first newdim columns of solverbasis
818 // initialize() will recognize that pointer bd_solver.V_ == pointer rstate.V, and will neglect the copy.
819 rstate.V = solverbasis;
820 }
821 else { // inSituRestart == false)
822# ifdef TEUCHOS_DEBUG
823 *out << "Beginning ex-situ restart...\n";
824# endif
825 // newV = oldV*Sr, explicitly. workspace is in workMV
826 std::vector<int> curind(curdim), newind(newdim);
827 for (int i=0; i<curdim; i++) curind[i] = i;
828 for (int i=0; i<newdim; i++) newind[i] = i;
829 Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
830 Teuchos::RCP<MV> newV = MVT::CloneViewNonConst(*workMV ,newind);
831
832 MVT::MvTimesMatAddMv(ONE,*oldV,Sr,ZERO,*newV);
833 //
834 // put the new basis into the state for initialize()
835 rstate.V = newV;
836 }
837
838 //
839 // send the new state to the solver
840 bd_solver->initialize(rstate);
841 } // end of restarting
843 //
844 // check locking if we didn't converge or restart
845 //
847 else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
848
849#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
850 Teuchos::TimeMonitor lcktimer(*_timerLocking);
851#endif
852
853 //
854 // get current state
855 BlockDavidsonState<ScalarType,MV> state = bd_solver->getState();
856 const int curdim = state.curDim;
857
858 //
859 // get number,indices of vectors to be locked
860 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
861 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() non-positive.");
862 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
863 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
864 // we should have room to lock vectors; otherwise, locking should have been deactivated
865 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
866 "Anasazi::BlockDavidsonSolMgr::solve(): status test mistake: locking not deactivated.");
867 //
868 // don't lock more than maxLocked_; we didn't allocate enough space.
869 std::vector<int> tmp_vector_int;
870 if (curNumLocked + locktest->howMany() > maxLocked_) {
871 // just use the first of them
872 tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
873 }
874 else {
875 tmp_vector_int = locktest->whichVecs();
876 }
877 const std::vector<int> lockind(tmp_vector_int);
878 const int numNewLocked = lockind.size();
879 //
880 // generate indices of vectors left unlocked
881 // curind = [0,...,curdim-1] = UNION( lockind, unlockind )
882 const int numUnlocked = curdim-numNewLocked;
883 tmp_vector_int.resize(curdim);
884 for (int i=0; i<curdim; i++) tmp_vector_int[i] = i;
885 const std::vector<int> curind(tmp_vector_int); // curind = [0 ... curdim-1]
886 tmp_vector_int.resize(numUnlocked);
887 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
888 const std::vector<int> unlockind(tmp_vector_int); // unlockind = [0 ... curdim-1] - lockind
889 tmp_vector_int.clear();
890
891 //
892 // debug printing
893 if (printer_->isVerbosity(Debug)) {
894 printer_->print(Debug,"Locking vectors: ");
895 for (unsigned int i=0; i<lockind.size(); i++) {printer_->stream(Debug) << " " << lockind[i];}
896 printer_->print(Debug,"\n");
897 printer_->print(Debug,"Not locking vectors: ");
898 for (unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(Debug) << " " << unlockind[i];}
899 printer_->print(Debug,"\n");
900 }
901
902 //
903 // we need primitive ritz vectors/values:
904 // [S,L] = eig(oldKK)
905 //
906 // this will be partitioned as follows:
907 // locked: Sl = S(lockind) // we won't actually need Sl
908 // unlocked: Su = S(unlockind)
909 //
910 Teuchos::SerialDenseMatrix<int,ScalarType> S(curdim,curdim);
911 std::vector<MagnitudeType> theta(curdim);
912 {
913 int rank = curdim;
914 int info = msutils::directSolver(curdim,*state.KK,Teuchos::null,S,theta,rank,10);
915 TEUCHOS_TEST_FOR_EXCEPTION(info != 0 ,std::logic_error,
916 "Anasazi::BlockDavidsonSolMgr::solve(): error calling SolverUtils::directSolver."); // this should never happen
917 TEUCHOS_TEST_FOR_EXCEPTION(rank != curdim,std::logic_error,
918 "Anasazi::BlockDavidsonSolMgr::solve(): direct solve did not compute all eigenvectors."); // this should never happen
919 //
920 // sort the eigenvalues (so that we can order the eigenvectors)
921 std::vector<int> order(curdim);
922 sorter->sort(theta,Teuchos::rcpFromRef(order),curdim);
923 //
924 // apply the same ordering to the primitive ritz vectors
925 msutils::permuteVectors(order,S);
926 }
927 //
928 // select the unlocked ritz vectors
929 // the indexing in unlockind is relative to the ordered primitive ritz vectors
930 // (this is why we ordered theta,S above)
931 Teuchos::SerialDenseMatrix<int,ScalarType> Su(curdim,numUnlocked);
932 for (int i=0; i<numUnlocked; i++) {
933 blas.COPY(curdim, S[unlockind[i]], 1, Su[i], 1);
934 }
935
936
937 //
938 // newV has the following form:
939 // newV = [defV augV]
940 // - defV will be of size curdim - numNewLocked, and contain the generated basis: defV = oldV*Su
941 // - augV will be of size numNewLocked, and contain random directions to make up for the lost space
942 //
943 // we will need a pointer to defV below to generate the off-diagonal block of newKK
944 // go ahead and setup pointer to augV
945 //
946 Teuchos::RCP<MV> defV, augV;
947 if (inSituRestart_ == true) {
948 //
949 // get non-const pointer to solver's basis so we can work in situ
950 Teuchos::RCP<MV> solverbasis = Teuchos::rcp_const_cast<MV>(state.V);
951 //
952 // perform Householder QR of Su = Q [D;0], where D is unit diag.
953 // work on a copy of Su, since we need Su below to build newKK
954 Teuchos::SerialDenseMatrix<int,ScalarType> copySu(Su);
955 std::vector<ScalarType> tau(numUnlocked), work(numUnlocked);
956 int info;
957 lapack.GEQRF(curdim,numUnlocked,copySu.values(),copySu.stride(),&tau[0],&work[0],work.size(),&info);
958 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,std::logic_error,
959 "Anasazi::BlockDavidsonSolMgr::solve(): error calling GEQRF during restarting.");
960 if (printer_->isVerbosity(Debug)) {
961 Teuchos::SerialDenseMatrix<int,ScalarType> R(Teuchos::Copy,copySu,numUnlocked,numUnlocked);
962 for (int j=0; j<numUnlocked; j++) {
963 R(j,j) = SCT::magnitude(R(j,j)) - 1.0;
964 for (int i=j+1; i<numUnlocked; i++) {
965 R(i,j) = ZERO;
966 }
967 }
968 printer_->stream(Debug) << "||Triangular factor of Su - I||: " << R.normFrobenius() << std::endl;
969 }
970 //
971 // perform implicit oldV*Su
972 // this actually performs oldV*[Su Sl*M] = [defV lockV], for some unitary M
973 // we are actually interested in only the first numUnlocked vectors of the result
974 {
975 Teuchos::RCP<MV> oldV = MVT::CloneViewNonConst(*solverbasis,curind);
976 msutils::applyHouse(numUnlocked,*oldV,copySu,tau,workMV);
977 }
978 std::vector<int> defind(numUnlocked), augind(numNewLocked);
979 for (int i=0; i<numUnlocked ; i++) defind[i] = i;
980 for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
981 defV = MVT::CloneViewNonConst(*solverbasis,defind);
982 augV = MVT::CloneViewNonConst(*solverbasis,augind);
983 }
984 else { // inSituRestart == false)
985 // defV = oldV*Su, explicitly. workspace is in workMV
986 std::vector<int> defind(numUnlocked), augind(numNewLocked);
987 for (int i=0; i<numUnlocked ; i++) defind[i] = i;
988 for (int i=0; i<numNewLocked; i++) augind[i] = numUnlocked+i;
989 Teuchos::RCP<const MV> oldV = MVT::CloneView(*state.V,curind);
990 defV = MVT::CloneViewNonConst(*workMV,defind);
991 augV = MVT::CloneViewNonConst(*workMV,augind);
992
993 MVT::MvTimesMatAddMv(ONE,*oldV,Su,ZERO,*defV);
994 }
995
996 //
997 // lockvecs will be partitioned as follows:
998 // lockvecs = [curlocked augTmp ...]
999 // - augTmp will be used for the storage of M*augV and K*augV
1000 // later, the locked vectors (stored in state.X and referenced via const MV view newLocked)
1001 // will be moved into lockvecs on top of augTmp when it is no longer needed as workspace.
1002 // - curlocked will be used in orthogonalization of augV
1003 //
1004 // newL is the new locked vectors; newL = oldV*Sl = RitzVectors(lockind)
1005 // we will not produce them, but instead retrieve them from RitzVectors
1006 //
1007 Teuchos::RCP<const MV> curlocked, newLocked;
1008 Teuchos::RCP<MV> augTmp;
1009 {
1010 // setup curlocked
1011 if (curNumLocked > 0) {
1012 std::vector<int> curlockind(curNumLocked);
1013 for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
1014 curlocked = MVT::CloneView(*lockvecs,curlockind);
1015 }
1016 else {
1017 curlocked = Teuchos::null;
1018 }
1019 // setup augTmp
1020 std::vector<int> augtmpind(numNewLocked);
1021 for (int i=0; i<numNewLocked; i++) augtmpind[i] = curNumLocked+i;
1022 augTmp = MVT::CloneViewNonConst(*lockvecs,augtmpind);
1023 // setup newLocked
1024 newLocked = MVT::CloneView(*bd_solver->getRitzVectors(),lockind);
1025 }
1026
1027 //
1028 // generate augV and perform orthogonalization
1029 //
1030 MVT::MvRandom(*augV);
1031 //
1032 // orthogonalize it against auxvecs, defV, and all locked vectors (new and current)
1033 // use augTmp as storage for M*augV, if hasM
1034 {
1035 Teuchos::Array<Teuchos::RCP<const MV> > against;
1036 Teuchos::Array<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > dummyC;
1037 if (probauxvecs != Teuchos::null) against.push_back(probauxvecs);
1038 if (curlocked != Teuchos::null) against.push_back(curlocked);
1039 against.push_back(newLocked);
1040 against.push_back(defV);
1041 if (problem_->getM() != Teuchos::null) {
1042 OPT::Apply(*problem_->getM(),*augV,*augTmp);
1043 }
1044 ortho->projectAndNormalizeMat(*augV,against,dummyC,Teuchos::null,augTmp);
1045 }
1046
1047 //
1048 // form newKK
1049 //
1050 // newKK = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
1051 // [augV'*K*defV augV'*K*augV]
1052 //
1053 // first, generate the principal submatrix, the projection of K onto the unlocked portion of oldV
1054 //
1055 Teuchos::SerialDenseMatrix<int,ScalarType> newKK(curdim,curdim);
1056 {
1057 Teuchos::SerialDenseMatrix<int,ScalarType> KKtmp(curdim,numUnlocked),
1058 KKold(Teuchos::View,*state.KK,curdim,curdim),
1059 KK11(Teuchos::View,newKK,numUnlocked,numUnlocked);
1060 int teuchosRet;
1061 // KKtmp = KKold*Su
1062 teuchosRet = KKtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,ONE,KKold,Su,ZERO);
1063 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1064 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1065 // KK11 = Su'*KKtmp = Su'*KKold*Su
1066 teuchosRet = KK11.multiply(Teuchos::CONJ_TRANS,Teuchos::NO_TRANS,ONE,Su,KKtmp,ZERO);
1067 TEUCHOS_TEST_FOR_EXCEPTION(teuchosRet != 0,std::logic_error,
1068 "Anasazi::BlockDavidsonSolMgr::solve(): Logic error calling SerialDenseMatrix::multiply.");
1069 }
1070 //
1071 // project the stiffness matrix on augV
1072 {
1073 OPT::Apply(*problem_->getOperator(),*augV,*augTmp);
1074 Teuchos::SerialDenseMatrix<int,ScalarType> KK12(Teuchos::View,newKK,numUnlocked,numNewLocked,0,numUnlocked),
1075 KK22(Teuchos::View,newKK,numNewLocked,numNewLocked,numUnlocked,numUnlocked);
1076 MVT::MvTransMv(ONE,*defV,*augTmp,KK12);
1077 MVT::MvTransMv(ONE,*augV,*augTmp,KK22);
1078 }
1079 //
1080 // done with defV,augV
1081 defV = Teuchos::null;
1082 augV = Teuchos::null;
1083 //
1084 // make it hermitian in memory (fill in KK21)
1085 for (int j=0; j<curdim; ++j) {
1086 for (int i=j+1; i<curdim; ++i) {
1087 newKK(i,j) = SCT::conjugate(newKK(j,i));
1088 }
1089 }
1090 //
1091 // we are done using augTmp as storage
1092 // put newLocked into lockvecs, new values into lockvals
1093 augTmp = Teuchos::null;
1094 {
1095 std::vector<Value<ScalarType> > allvals = bd_solver->getRitzValues();
1096 for (int i=0; i<numNewLocked; i++) {
1097 lockvals.push_back(allvals[lockind[i]].realpart);
1098 }
1099
1100 std::vector<int> indlock(numNewLocked);
1101 for (int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
1102 MVT::SetBlock(*newLocked,indlock,*lockvecs);
1103 newLocked = Teuchos::null;
1104
1105 curNumLocked += numNewLocked;
1106 std::vector<int> curlockind(curNumLocked);
1107 for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
1108 curlocked = MVT::CloneView(*lockvecs,curlockind);
1109 }
1110 // add locked vecs as aux vecs, along with aux vecs from problem
1111 // add lockvals to ordertest
1112 // disable locktest if curNumLocked == maxLocked
1113 {
1114 ordertest->setAuxVals(lockvals);
1115
1116 Teuchos::Array< Teuchos::RCP<const MV> > aux;
1117 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
1118 aux.push_back(curlocked);
1119 bd_solver->setAuxVecs(aux);
1120
1121 if (curNumLocked == maxLocked_) {
1122 // disabled locking now
1123 combotest->removeTest(locktest);
1124 }
1125 }
1126
1127 //
1128 // prepare new state
1130 rstate.curDim = curdim;
1131 if (inSituRestart_) {
1132 // data is already in the solver's memory
1133 rstate.V = state.V;
1134 }
1135 else {
1136 // data is in workspace and will be copied to solver memory
1137 rstate.V = workMV;
1138 }
1139 rstate.KK = Teuchos::rcpFromRef(newKK);
1140 //
1141 // pass new state to the solver
1142 bd_solver->initialize(rstate);
1143 } // end of locking
1145 //
1146 // we returned from iterate(), but none of our status tests Passed.
1147 // something is wrong, and it is probably our fault.
1148 //
1150 else {
1151 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): Invalid return from bd_solver::iterate().");
1152 }
1153 }
1154 catch (const AnasaziError &err) {
1155 printer_->stream(Errors)
1156 << "Anasazi::BlockDavidsonSolMgr::solve() caught unexpected exception from Anasazi::BlockDavidson::iterate() at iteration " << bd_solver->getNumIters() << std::endl
1157 << err.what() << std::endl
1158 << "Anasazi::BlockDavidsonSolMgr::solve() returning Unconverged with no solutions." << std::endl;
1159 return Unconverged;
1160 }
1161 }
1162
1163 // clear temp space
1164 workMV = Teuchos::null;
1165
1166 sol.numVecs = ordertest->howMany();
1167 if (sol.numVecs > 0) {
1168 sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
1169 sol.Espace = sol.Evecs;
1170 sol.Evals.resize(sol.numVecs);
1171 std::vector<MagnitudeType> vals(sol.numVecs);
1172
1173 // copy them into the solution
1174 std::vector<int> which = ordertest->whichVecs();
1175 // indices between [0,blockSize) refer to vectors/values in the solver
1176 // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
1177 // everything has already been ordered by the solver; we just have to partition the two references
1178 std::vector<int> inlocked(0), insolver(0);
1179 for (unsigned int i=0; i<which.size(); i++) {
1180 if (which[i] >= 0) {
1181 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): positive indexing mistake from ordertest.");
1182 insolver.push_back(which[i]);
1183 }
1184 else {
1185 // sanity check
1186 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): negative indexing mistake from ordertest.");
1187 inlocked.push_back(which[i] + curNumLocked);
1188 }
1189 }
1190
1191 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::BlockDavidsonSolMgr::solve(): indexing mistake.");
1192
1193 // set the vecs,vals in the solution
1194 if (insolver.size() > 0) {
1195 // set vecs
1196 int lclnum = insolver.size();
1197 std::vector<int> tosol(lclnum);
1198 for (int i=0; i<lclnum; i++) tosol[i] = i;
1199 Teuchos::RCP<const MV> v = MVT::CloneView(*bd_solver->getRitzVectors(),insolver);
1200 MVT::SetBlock(*v,tosol,*sol.Evecs);
1201 // set vals
1202 std::vector<Value<ScalarType> > fromsolver = bd_solver->getRitzValues();
1203 for (unsigned int i=0; i<insolver.size(); i++) {
1204 vals[i] = fromsolver[insolver[i]].realpart;
1205 }
1206 }
1207
1208 // get the vecs,vals from locked storage
1209 if (inlocked.size() > 0) {
1210 int solnum = insolver.size();
1211 // set vecs
1212 int lclnum = inlocked.size();
1213 std::vector<int> tosol(lclnum);
1214 for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1215 Teuchos::RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1216 MVT::SetBlock(*v,tosol,*sol.Evecs);
1217 // set vals
1218 for (unsigned int i=0; i<inlocked.size(); i++) {
1219 vals[i+solnum] = lockvals[inlocked[i]];
1220 }
1221 }
1222
1223 // sort the eigenvalues and permute the eigenvectors appropriately
1224 {
1225 std::vector<int> order(sol.numVecs);
1226 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
1227 // store the values in the Eigensolution
1228 for (int i=0; i<sol.numVecs; i++) {
1229 sol.Evals[i].realpart = vals[i];
1230 sol.Evals[i].imagpart = MT::zero();
1231 }
1232 // now permute the eigenvectors according to order
1233 msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1234 }
1235
1236 // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1237 sol.index.resize(sol.numVecs,0);
1238 }
1239 }
1240
1241 // print final summary
1242 bd_solver->currentStatus(printer_->stream(FinalSummary));
1243
1244 // print timing information
1245#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1246 if ( printer_->isVerbosity( TimingDetails ) ) {
1247 Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
1248 }
1249#endif
1250
1251 problem_->setSolution(sol);
1252 printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1253
1254 // get the number of iterations taken for this call to solve().
1255 numIters_ = bd_solver->getNumIters();
1256
1257 if (sol.numVecs < nev) {
1258 return Unconverged; // return from BlockDavidsonSolMgr::solve()
1259 }
1260 return Converged; // return from BlockDavidsonSolMgr::solve()
1261}
1262
1263
1264template <class ScalarType, class MV, class OP>
1265void
1267 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &global)
1268{
1269 globalTest_ = global;
1270}
1271
1272template <class ScalarType, class MV, class OP>
1273const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1275{
1276 return globalTest_;
1277}
1278
1279template <class ScalarType, class MV, class OP>
1280void
1282 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &debug)
1283{
1284 debugTest_ = debug;
1285}
1286
1287template <class ScalarType, class MV, class OP>
1288const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1290{
1291 return debugTest_;
1292}
1293
1294template <class ScalarType, class MV, class OP>
1295void
1297 const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &locking)
1298{
1299 lockingTest_ = locking;
1300}
1301
1302template <class ScalarType, class MV, class OP>
1303const Teuchos::RCP< StatusTest<ScalarType,MV,OP> > &
1305{
1306 return lockingTest_;
1307}
1308
1309} // end Anasazi namespace
1310
1311#endif /* ANASAZI_BLOCKDAVIDSON_SOLMGR_HPP */
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::SortManager class.
Implementation of the block Davidson method.
Anasazi header file which uses auto-configuration information to include necessary C++ headers.
Abstract base class which defines the interface required by an eigensolver and status test class to c...
Abstract class definition for Anasazi Output Managers.
Abstract class definition for Anasazi output stream.
Orthogonalization manager based on the SVQB technique described in "A Block Orthogonalization Procedu...
Pure virtual base class which describes the basic interface for a solver manager.
Class which provides internal utilities for the Anasazi solvers.
Status test for forming logical combinations of other status tests.
Special StatusTest for printing status tests.
A status test for testing the norm of the eigenvectors residuals.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Types and exceptions used within Anasazi solvers and interfaces.
An exception class parent to all Anasazi exceptions.
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using (potentially)...
An implementation of the Anasazi::SortManager that performs a collection of common sorting techniques...
The BlockDavidsonSolMgr provides a powerful solver manager over the BlockDavidson eigensolver.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
BlockDavidsonSolMgr(const Teuchos::RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for BlockDavidsonSolMgr.
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
void setLockingStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
int getNumIters() const
Get the iteration count for the most recent call to solve().
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
void setGlobalStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
This class implements a Block Davidson iteration, a preconditioned iteration for solving linear Hermi...
This class defines the interface required by an eigensolver and status test class to compute solution...
Traits class which defines basic operations on multivectors.
static ptrdiff_t GetGlobalLength(const MV &mv)
Return the number of rows in the given multivector mv.
Virtual base class which defines basic traits for the operator type.
Output managers remove the need for the eigensolver to know any information about the required output...
An implementation of the Anasazi::MatOrthoManager that performs orthogonalization using the SVQB iter...
The Anasazi::SolverManager is a templated virtual base class that defines the basic interface that an...
Anasazi's templated, static class providing utilities for the solvers.
Status test for forming logical combinations of other status tests.
A special StatusTest for printing other status tests.
A status test for testing the norm of the eigenvectors residuals.
A status test for testing the norm of the eigenvectors residuals along with a set of auxiliary eigenv...
Common interface of stopping criteria for Anasazi's solvers.
Namespace Anasazi contains the classes, structs, enums and utilities used by the Anasazi package.
ReturnType
Enumerated type used to pass back information from a solver manager.
ResType
Enumerated type used to specify which residual norm used by residual norm status tests.
@ IterationDetails
Structure to contain pointers to BlockDavidson state variables.
int curDim
The current dimension of the solver.
Teuchos::RCP< const MV > V
The basis for the Krylov space.
Teuchos::RCP< const MV > X
The current eigenvectors.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
Teuchos::RCP< const MV > R
The current residual vectors.
Teuchos::RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
Teuchos::RCP< const MV > KX
The image of the current eigenvectors under K.
Teuchos::RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
Struct for storing an eigenproblem solution.
Teuchos::RCP< MV > Evecs
The computed eigenvectors.
std::vector< int > index
An index into Evecs to allow compressed storage of eigenvectors for real, non-Hermitian problems.
int numVecs
The number of computed eigenpairs.
Teuchos::RCP< MV > Espace
An orthonormal basis for the computed eigenspace.
std::vector< Value< ScalarType > > Evals
The computed eigenvalues.
Output managers remove the need for the eigensolver to know any information about the required output...