Anasazi Version of the Day
Loading...
Searching...
No Matches
AnasaziTraceMinBaseSolMgr.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Anasazi: Block Eigensolvers Package
4//
5// Copyright 2004 NTESS and the Anasazi contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef ANASAZI_TraceMinBase_SOLMGR_HPP
11#define ANASAZI_TraceMinBase_SOLMGR_HPP
12
20#include "AnasaziBasicSort.hpp"
21#include "AnasaziConfigDefs.hpp"
29#include "AnasaziStatusTestSpecTrans.hpp"
34#include "AnasaziTypes.hpp"
35
36#include "Teuchos_TimeMonitor.hpp"
37#include "Teuchos_FancyOStream.hpp"
38
39using Teuchos::RCP;
40using Teuchos::rcp;
41
42namespace Anasazi {
43namespace Experimental {
44
45
78template<class ScalarType, class MV, class OP>
79class TraceMinBaseSolMgr : public SolverManager<ScalarType,MV,OP> {
80
81 private:
84 typedef Teuchos::ScalarTraits<ScalarType> SCT;
85 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
86 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
87
88 public:
89
91
92
139 Teuchos::ParameterList &pl );
140
144
146
147
150 return *problem_;
151 }
152
154 int getNumIters() const {
155 return numIters_;
156 }
157
165 Teuchos::Array<RCP<Teuchos::Time> > getTimers() const {
166 return Teuchos::tuple(_timerSolve, _timerRestarting, _timerLocking);
167 }
168
170
172
173
198
200 void setGlobalStatusTest(const RCP< StatusTest<ScalarType,MV,OP> > &global);
201
203 const RCP< StatusTest<ScalarType,MV,OP> > & getGlobalStatusTest() const;
204
206 void setLockingStatusTest(const RCP< StatusTest<ScalarType,MV,OP> > &locking);
207
209 const RCP< StatusTest<ScalarType,MV,OP> > & getLockingStatusTest() const;
210
212 void setDebugStatusTest(const RCP< StatusTest<ScalarType,MV,OP> > &debug);
213
215 const RCP< StatusTest<ScalarType,MV,OP> > & getDebugStatusTest() const;
216
218
219 protected:
220 RCP<Eigenproblem<ScalarType,MV,OP> > problem_;
221
222 int numIters_;
223
224 // Block variables
225 int blockSize_, numBlocks_, numRestartBlocks_;
226
227 // Output variables
228 RCP<OutputManager<ScalarType> > printer_;
229
230 // Convergence variables
231 MagnitudeType convTol_;
232 bool relConvTol_;
233 enum ResType convNorm_;
234
235 // Locking variables
236 MagnitudeType lockTol_;
237 int maxLocked_, lockQuorum_;
238 bool useLocking_, relLockTol_;
239 enum ResType lockNorm_;
240
241 // Shifting variables
242 enum WhenToShiftType whenToShift_;
243 MagnitudeType traceThresh_, shiftTol_;
244 enum HowToShiftType howToShift_;
245 bool useMultipleShifts_, relShiftTol_, considerClusters_;
246 std::string shiftNorm_;
247
248 // Other variables
249 int maxKrylovIter_;
250 std::string ortho_, which_;
251 enum SaddleSolType saddleSolType_;
252 bool projectAllVecs_, projectLockedVecs_, computeAllRes_, useRHSR_, useHarmonic_, noSort_;
253 MagnitudeType alpha_;
254
255 // Timers
256 RCP<Teuchos::Time> _timerSolve, _timerRestarting, _timerLocking;
257
258 // Status tests
259 RCP<StatusTest<ScalarType,MV,OP> > globalTest_;
260 RCP<StatusTest<ScalarType,MV,OP> > lockingTest_;
261 RCP<StatusTest<ScalarType,MV,OP> > debugTest_;
262
263 // TraceMin specific functions
264 void copyPartOfState(const TraceMinBaseState<ScalarType,MV>& oldState, TraceMinBaseState<ScalarType,MV>& newState, const std::vector<int> indToCopy) const;
265
266 void setParameters(Teuchos::ParameterList &pl) const;
267
268 void printParameters(std::ostream &os) const;
269
270 virtual RCP< TraceMinBase<ScalarType,MV,OP> > createSolver(
271 const RCP<SortManager<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> > &sorter,
272 const RCP<StatusTest<ScalarType,MV,OP> > &outputtest,
273 const RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
274 Teuchos::ParameterList &plist
275 ) =0;
276
277 virtual bool needToRestart(const RCP< TraceMinBase<ScalarType,MV,OP> > solver) =0;
278
279 virtual bool performRestart(int &numRestarts, RCP< TraceMinBase<ScalarType,MV,OP> > solver) =0;
280};
281
282
285// Constructor
286template<class ScalarType, class MV, class OP>
288 const RCP<Eigenproblem<ScalarType,MV,OP> > &problem,
289 Teuchos::ParameterList &pl ) :
290 problem_(problem)
291#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
292 , _timerSolve(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBaseSolMgr::solve()")),
293 _timerRestarting(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBaseSolMgr restarting")),
294 _timerLocking(Teuchos::TimeMonitor::getNewTimer("Anasazi: TraceMinBaseSolMgr locking"))
295#endif
296{
297 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
298 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(), std::invalid_argument, "Problem not set.");
299 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isHermitian(), std::invalid_argument, "Problem not symmetric.");
300 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getInitVec() == Teuchos::null,std::invalid_argument, "Problem does not contain initial vectors to clone from.");
301
302 std::string strtmp;
303
305 // Block parameters
306
307 // TODO: the default is different for TraceMin and TraceMin-Davidson
308 // block size: default is nev()
309// blockSize_ = pl.get("Block Size",problem_->getNEV());
310// TEUCHOS_TEST_FOR_EXCEPTION(blockSize_ <= 0, std::invalid_argument,
311// "Anasazi::TraceMinBaseSolMgr: \"Block Size\" must be strictly positive.");
312
313 // TODO: add Num Blocks as a parameter to both child classes, since they have different default values
314// numBlocks_ = pl.get("Num Blocks",5);
315// TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ < 2, std::invalid_argument,
316// "Anasazi::TraceMinBaseSolMgr: \"Num Blocks\" must be >= 2.");
317
319 // Output parameters
320
321 // Create a formatted output stream to print to.
322 // See if user requests output processor.
323 int osProc = pl.get("Output Processor", 0);
324
325 // If not passed in by user, it will be chosen based upon operator type.
326 Teuchos::RCP<Teuchos::FancyOStream> osp;
327
328 if (pl.isParameter("Output Stream")) {
329 osp = Teuchos::getParameter<Teuchos::RCP<Teuchos::FancyOStream> >(pl,"Output Stream");
330 }
331 else {
332 osp = OutputStreamTraits<OP>::getOutputStream (*problem_->getOperator(), osProc);
333 }
334
335 int verbosity = Anasazi::Errors;
336 if (pl.isParameter("Verbosity")) {
337 if (Teuchos::isParameterType<int>(pl,"Verbosity")) {
338 verbosity = pl.get("Verbosity", verbosity);
339 } else {
340 verbosity = (int)Teuchos::getParameter<Anasazi::MsgType>(pl,"Verbosity");
341 }
342 }
343 printer_ = rcp( new OutputManager<ScalarType>(verbosity,osp) );
344
345 // TODO: Add restart parameters to TraceMin-Davidson
346
348 // Convergence parameters
349 convTol_ = pl.get("Convergence Tolerance",MT::prec());
350 TEUCHOS_TEST_FOR_EXCEPTION(convTol_ < 0, std::invalid_argument,
351 "Anasazi::TraceMinBaseSolMgr: \"Convergence Tolerance\" must be nonnegative.");
352
353 relConvTol_ = pl.get("Relative Convergence Tolerance",true);
354 strtmp = pl.get("Convergence Norm",std::string("2"));
355 if (strtmp == "2") {
356 convNorm_ = RES_2NORM;
357 }
358 else if (strtmp == "M") {
359 convNorm_ = RES_ORTH;
360 }
361 else {
362 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
363 "Anasazi::TraceMinBaseSolMgr: Invalid Convergence Norm.");
364 }
365
367 // Locking parameters
368 useLocking_ = pl.get("Use Locking",true);
369 relLockTol_ = pl.get("Relative Locking Tolerance",true);
370 lockTol_ = pl.get("Locking Tolerance",convTol_/10);
371
372 TEUCHOS_TEST_FOR_EXCEPTION(relConvTol_ != relLockTol_, std::invalid_argument,
373 "Anasazi::TraceMinBaseSolMgr: \"Relative Convergence Tolerance\" and \"Relative Locking Tolerance\" have different values. If you set one, you should always set the other.");
374
375 TEUCHOS_TEST_FOR_EXCEPTION(lockTol_ < 0, std::invalid_argument,
376 "Anasazi::TraceMinBaseSolMgr: \"Locking Tolerance\" must be nonnegative.");
377
378 strtmp = pl.get("Locking Norm",std::string("2"));
379 if (strtmp == "2") {
380 lockNorm_ = RES_2NORM;
381 }
382 else if (strtmp == "M") {
383 lockNorm_ = RES_ORTH;
384 }
385 else {
386 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
387 "Anasazi::TraceMinBaseSolMgr: Invalid Locking Norm.");
388 }
389
390 // max locked: default is nev(), must satisfy maxLocked_ + blockSize_ >= nev
391 if (useLocking_) {
392 maxLocked_ = pl.get("Max Locked",problem_->getNEV());
393 TEUCHOS_TEST_FOR_EXCEPTION(maxLocked_ <= 0, std::invalid_argument,
394 "Anasazi::TraceMinBaseSolMgr: \"Max Locked\" must be strictly positive.");
395 }
396 else {
397 maxLocked_ = 0;
398 }
399
400 if (useLocking_) {
401 lockQuorum_ = pl.get("Locking Quorum",1);
402 TEUCHOS_TEST_FOR_EXCEPTION(lockQuorum_ <= 0, std::invalid_argument,
403 "Anasazi::TraceMinBaseSolMgr: \"Locking Quorum\" must be strictly positive.");
404 }
405
407 // Ritz shift parameters
408
409 // When to shift - what triggers a shift?
410 strtmp = pl.get("When To Shift", "Always");
411
412 if(strtmp == "Never")
413 whenToShift_ = NEVER_SHIFT;
414 else if(strtmp == "After Trace Levels")
415 whenToShift_ = SHIFT_WHEN_TRACE_LEVELS;
416 else if(strtmp == "Residual Becomes Small")
417 whenToShift_ = SHIFT_WHEN_RESID_SMALL;
418 else if(strtmp == "Always")
419 whenToShift_ = ALWAYS_SHIFT;
420 else
421 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
422 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"When To Shift\"; valid options are \"Never\", \"After Trace Levels\", \"Residual Becomes Small\", \"Always\".");
423
424
425 // How small does the % change in trace have to get before shifting?
426 traceThresh_ = pl.get("Trace Threshold", 0.02);
427
428 TEUCHOS_TEST_FOR_EXCEPTION(traceThresh_ < 0, std::invalid_argument,
429 "Anasazi::TraceMinBaseSolMgr: \"Trace Threshold\" must be nonnegative.");
430
431 // Shift threshold - if the residual of an eigenpair is less than this, then shift
432 shiftTol_ = pl.get("Shift Tolerance", 0.1);
433
434 TEUCHOS_TEST_FOR_EXCEPTION(shiftTol_ < 0, std::invalid_argument,
435 "Anasazi::TraceMinBaseSolMgr: \"Shift Tolerance\" must be nonnegative.");
436
437 // Use relative convergence tolerance - scale by eigenvalue?
438 relShiftTol_ = pl.get("Relative Shift Tolerance", true);
439
440 // Which norm to use in determining whether to shift
441 shiftNorm_ = pl.get("Shift Norm", "2");
442
443 TEUCHOS_TEST_FOR_EXCEPTION(shiftNorm_ != "2" && shiftNorm_ != "M", std::invalid_argument,
444 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Shift Norm\"; valid options are \"2\" and \"M\".");
445
446 noSort_ = pl.get("No Sorting", false);
447
448 // How to choose shift
449 strtmp = pl.get("How To Choose Shift", "Adjusted Ritz Values");
450
451 if(strtmp == "Largest Converged")
452 howToShift_ = LARGEST_CONVERGED_SHIFT;
453 else if(strtmp == "Adjusted Ritz Values")
454 howToShift_ = ADJUSTED_RITZ_SHIFT;
455 else if(strtmp == "Ritz Values")
456 howToShift_ = RITZ_VALUES_SHIFT;
457 else if(strtmp == "Experimental Shift")
458 howToShift_ = EXPERIMENTAL_SHIFT;
459 else
460 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
461 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"How To Choose Shift\"; valid options are \"Largest Converged\", \"Adjusted Ritz Values\", \"Ritz Values\".");
462
463 // Consider clusters - if all eigenvalues are in one cluster, it's not expecially safe to shift
464 considerClusters_ = pl.get("Consider Clusters", true);
465
466 // Use multiple shifts
467 useMultipleShifts_ = pl.get("Use Multiple Shifts", true);
468
470 // Other parameters
471
472 // which orthogonalization to use
473 ortho_ = pl.get("Orthogonalization", "SVQB");
474 TEUCHOS_TEST_FOR_EXCEPTION(ortho_ != "DGKS" && ortho_ != "SVQB" && ortho_ != "ICGS", std::invalid_argument,
475 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Orthogonalization\"; valid options are \"DGKS\", \"SVQB\", \"ICGS\".");
476
477 strtmp = pl.get("Saddle Solver Type", "Projected Krylov");
478
479 if(strtmp == "Projected Krylov")
480 saddleSolType_ = PROJECTED_KRYLOV_SOLVER;
481 else if(strtmp == "Schur Complement")
482 saddleSolType_ = SCHUR_COMPLEMENT_SOLVER;
483 else if(strtmp == "Block Diagonal Preconditioned Minres")
484 saddleSolType_ = BD_PREC_MINRES;
485 else if(strtmp == "HSS Preconditioned Gmres")
486 saddleSolType_ = HSS_PREC_GMRES;
487 else
488 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
489 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Saddle Solver Type\"; valid options are \"Projected Krylov\", \"Schur Complement\", and \"Block Diagonal Preconditioned Minres\".");
490
491 projectAllVecs_ = pl.get("Project All Vectors", true);
492 projectLockedVecs_ = pl.get("Project Locked Vectors", true);
493 computeAllRes_ = pl.get("Compute All Residuals", true);
494 useRHSR_ = pl.get("Use Residual as RHS", false);
495 alpha_ = pl.get("HSS: alpha", 1.0);
496
497 TEUCHOS_TEST_FOR_EXCEPTION(projectLockedVecs_ && ! projectAllVecs_, std::invalid_argument,
498 "Anasazi::TraceMinBaseSolMgr: If you want to project out the locked vectors, you should really project out ALL the vectors of X.");
499
500 // Maximum number of inner iterations
501 maxKrylovIter_ = pl.get("Maximum Krylov Iterations", 200);
502 TEUCHOS_TEST_FOR_EXCEPTION(maxKrylovIter_ < 1, std::invalid_argument,
503 "Anasazi::TraceMinBaseSolMgr: \"Maximum Krylov Iterations\" must be greater than 0.");
504
505 // Which eigenvalues we want to get
506 which_ = pl.get("Which", "SM");
507 TEUCHOS_TEST_FOR_EXCEPTION(which_ != "SM" && which_ != "LM", std::invalid_argument,
508 "Anasazi::TraceMinBaseSolMgr: Invalid value for \"Which\"; valid options are \"SM\" and \"LM\".");
509
510 // Test whether we are shifting without an operator K
511 // This is a really bad idea
512 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getOperator() == Teuchos::null && whenToShift_ != NEVER_SHIFT, std::invalid_argument,
513 "Anasazi::TraceMinBaseSolMgr: It is an exceptionally bad idea to use Ritz shifts when finding the largest eigenpairs of a standard eigenvalue problem. If we don't use Ritz shifts, it may take extra iterations to converge, but we NEVER have to solve a single linear system. Using Ritz shifts forces us to solve systems of the form (I + sigma A)x=f, and it probably doesn't benefit us enough to outweigh the extra cost. We may add support for this feature in the future, but for now, please set \"When To Shift\" to \"Never\".");
514
515#ifdef BELOS_PSEUDO_BLOCK_GMRES_SOLMGR_HPP
516 // Test whether we are using a projected preconditioner with multiple Ritz shifts
517 // We can't currently do this for reasons that are complicated and are explained in the user manual
518 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER && useMultipleShifts_, std::invalid_argument,
519 "Anasazi::TraceMinBaseSolMgr: When you use the projected Krylov solver with preconditioning, the preconditioner must be projected as well. In theory, if the preconditioner is SPD, the projected preconditioner will be SPSD, but in practice, it can have small negative eigenvalues, presumably due to machine arithmetic. This means we can't use TraceMin's built-in MINRES, and we are forced to use Belos for now. When you use multiple Ritz shifts, you are essentially using a different operator to solve each linear system. Belos can't handle this right now, but we're working on a solution. For now, please set \"Use Multiple Shifts\" to false.");
520#else
521 // Test whether we are using a projected preconditioner without Belos.
522 // P Prec P should be positive definite if Prec is positive-definite,
523 // but it tends not to be in practice, presumably due to machine arithmetic
524 // As a result, we have to use pseudo-block gmres for now.
525 // Make sure it's available.
526 TEUCHOS_TEST_FOR_EXCEPTION(problem_->getPrec() != Teuchos::null && saddleSolType_ == PROJECTED_KRYLOV_SOLVER, std::invalid_argument,
527 "Anasazi::TraceMinBaseSolMgr: When you use the projected Krylov solver with preconditioning, the preconditioner must be projected as well. In theory, if the preconditioner is SPD, the projected preconditioner will be SPSD, but in practice, it can have small negative eigenvalues, presumably due to machine arithmetic. This means we can't use TraceMin's built-in MINRES, and we are forced to use Belos for now. You didn't install Belos. You have three options to correct this problem:\n1. Reinstall Trilinos with Belos enabled\n2. Don't use a preconditioner\n3. Choose a different method for solving the saddle-point problem (Recommended)");
528
529
530#endif
531
532
533}
534
535
538// solve()
539template<class ScalarType, class MV, class OP>
542{
543 typedef SolverUtils<ScalarType,MV,OP> msutils;
544
545 const int nev = problem_->getNEV();
546
547#ifdef TEUCHOS_DEBUG
548 RCP<Teuchos::FancyOStream>
549 out = Teuchos::getFancyOStream(Teuchos::rcpFromRef(printer_->stream(Debug)));
550 out->setShowAllFrontMatter(false).setShowProcRank(true);
551 *out << "Entering Anasazi::TraceMinBaseSolMgr::solve()\n";
552#endif
553
555 // Sort manager
556 RCP<BasicSort<MagnitudeType> > sorter = rcp( new BasicSort<MagnitudeType>("SM") );
557
559 // Handle the spectral transformation if necessary
560 // TODO: Make sure we undo this before returning...
561 if(which_ == "LM")
562 {
563 RCP<const OP> swapHelper = problem_->getOperator();
564 problem_->setOperator(problem_->getM());
565 problem_->setM(swapHelper);
566 problem_->setProblem();
567 }
568
570 // Status tests
571 //
572 // convergence
573 RCP<StatusTest<ScalarType,MV,OP> > convtest;
574 if (globalTest_ == Teuchos::null) {
575 if(which_ == "SM")
576 convtest = rcp( new StatusTestResNorm<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_) );
577 else
578 convtest = rcp( new StatusTestSpecTrans<ScalarType,MV,OP>(convTol_,nev,convNorm_,relConvTol_,true,problem_->getOperator()) );
579 }
580 else {
581 convtest = globalTest_;
582 }
583 RCP<StatusTestWithOrdering<ScalarType,MV,OP> > ordertest
584 = rcp( new StatusTestWithOrdering<ScalarType,MV,OP>(convtest,sorter,nev) );
585 // locking
586 RCP<StatusTest<ScalarType,MV,OP> > locktest;
587 if (useLocking_) {
588 if (lockingTest_ == Teuchos::null) {
589 if(which_ == "SM")
590 locktest = rcp( new StatusTestResNorm<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_) );
591 else
592 locktest = rcp( new StatusTestSpecTrans<ScalarType,MV,OP>(lockTol_,lockQuorum_,lockNorm_,relLockTol_,true,problem_->getOperator()) );
593 }
594 else {
595 locktest = lockingTest_;
596 }
597 }
598 // for a non-short-circuited OR test, the order doesn't matter
599 Teuchos::Array<RCP<StatusTest<ScalarType,MV,OP> > > alltests;
600 alltests.push_back(ordertest);
601 if (locktest != Teuchos::null) alltests.push_back(locktest);
602 if (debugTest_ != Teuchos::null) alltests.push_back(debugTest_);
603
604 RCP<StatusTestCombo<ScalarType,MV,OP> > combotest
606 // printing StatusTest
607 RCP<StatusTestOutput<ScalarType,MV,OP> > outputtest;
608 if ( printer_->isVerbosity(Debug) ) {
609 outputtest = rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed+Failed+Undefined ) );
610 }
611 else {
612 outputtest = rcp( new StatusTestOutput<ScalarType,MV,OP>( printer_,combotest,1,Passed ) );
613 }
614
616 // Orthomanager
617 RCP<MatOrthoManager<ScalarType,MV,OP> > ortho;
618 if (ortho_=="SVQB") {
619 ortho = rcp( new SVQBOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
620 } else if (ortho_=="DGKS") {
621 ortho = rcp( new BasicOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
622 } else if (ortho_=="ICGS") {
623 ortho = rcp( new ICGSOrthoManager<ScalarType,MV,OP>(problem_->getM()) );
624 } else {
625 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): Invalid orthogonalization type.");
626 }
627
629 // Parameter list
630 Teuchos::ParameterList plist;
631 setParameters(plist);
632
634 // TraceMinBase solver
635 RCP<TraceMinBase<ScalarType,MV,OP> > tm_solver
636 = createSolver(sorter,outputtest,ortho,plist);
637 // set any auxiliary vectors defined in the problem
638 RCP< const MV > probauxvecs = problem_->getAuxVecs();
639 if (probauxvecs != Teuchos::null) {
640 tm_solver->setAuxVecs( Teuchos::tuple< RCP<const MV> >(probauxvecs) );
641 }
642
644 // Storage
645 //
646 // lockvecs will contain eigenvectors that have been determined "locked" by the status test
647 int curNumLocked = 0;
648 RCP<MV> lockvecs;
649 // lockvecs is used to hold the locked eigenvectors, as well as for temporary storage when locking.
650 // when locking, we will lock some number of vectors numnew, where numnew <= maxlocked - curlocked
651 // we will produce numnew random vectors, which will go into the space with the new basis.
652 // we will also need numnew storage for the image of these random vectors under A and M;
653 // columns [curlocked+1,curlocked+numnew] will be used for this storage
654 if (maxLocked_ > 0) {
655 lockvecs = MVT::Clone(*problem_->getInitVec(),maxLocked_);
656 }
657 std::vector<MagnitudeType> lockvals;
658 //
659 // Restarting occurs under two scenarios: when the basis is full and after locking.
660 //
661 // For the former, a new basis of size blockSize*numRestartBlocks is generated using the current basis
662 // and the most significant primitive Ritz vectors (projected eigenvectors).
663 // [S,L] = eig(KK)
664 // S = [Sr St] // some for "r"estarting, some are "t"runcated
665 // newV = V*Sr
666 // KK_new = newV'*K*newV = Sr'*V'*K*V*Sr = Sr'*KK*Sr
667 // Therefore, the only multivector operation needed is for the generation of newV.
668 //
669 // * If the multiplication is explicit, it requires a workspace of blockSize*numRestartBlocks vectors.
670 // This space must be specifically allocated for that task, as we don't have any space of that size.
671 // It (workMV) will be allocated at the beginning of solve()
672 // * Optionally, the multiplication can be performed implicitly, via a Householder QR factorization of
673 // Sr. This can be done in situ, using the basis multivector contained in the solver. This requires
674 // that we cast away the const on the multivector returned from getState(). Workspace for this approach
675 // is a single vector. the solver's internal storage must be preserved (X,MX,KX,R), requiring us to
676 // allocate this vector.
677 //
678 // For the latter (restarting after locking), the new basis is the same size as existing basis. If numnew
679 // vectors are locked, they are deflated from the current basis and replaced with randomly generated
680 // vectors.
681 // [S,L] = eig(KK)
682 // S = [Sl Su] // partitioned: "l"ocked and "u"nlocked
683 // newL = V*Sl = X(locked)
684 // defV = V*Su
685 // augV = rand(numnew) // orthogonal to oldL,newL,defV,auxvecs
686 // newV = [defV augV]
687 // Kknew = newV'*K*newV = [Su'*KK*Su defV'*K*augV]
688 // [augV'*K*defV augV'*K*augV]
689 // locked = [oldL newL]
690 // Clearly, this operation is more complicated than the previous.
691 // Here is a list of the significant computations that need to be performed:
692 // - newL will be put into space in lockvecs, but will be copied from getState().X at the end
693 // - defV,augV will be stored in workspace the size of the current basis.
694 // - M*augV and K*augV are needed; they will be stored in lockvecs. As a result, newL will
695 // not be put into lockvecs until the end.
696 //
697 // Therefore, we must allocate workMV when ((maxRestarts_ > 0) || (useLocking_ == true)) && inSituRestart == false
698 // It will be allocated to size (numBlocks-1)*blockSize
699 //
700
701 // some consts and utils
702 const ScalarType ONE = SCT::one();
703 const ScalarType ZERO = SCT::zero();
704
705 // go ahead and initialize the solution to nothing in case we throw an exception
707 sol.numVecs = 0;
708 problem_->setSolution(sol);
709
710 int numRestarts = 0;
711
712 // enter solve() iterations
713 {
714#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
715 Teuchos::TimeMonitor slvtimer(*_timerSolve);
716#endif
717
718 // tell tm_solver to iterate
719 while (1) {
720 try {
721 tm_solver->iterate();
722
724 //
725 //
727 if (debugTest_ != Teuchos::null && debugTest_->getStatus() == Passed) {
728 throw AnasaziError("Anasazi::TraceMinBaseSolMgr::solve(): User-specified debug status test returned Passed.");
729 }
731 //
732 // check convergence next
733 //
735 else if (ordertest->getStatus() == Passed ) {
736 // we have convergence
737 // ordertest->whichVecs() tells us which vectors from lockvecs and solver state are the ones we want
738 // ordertest->howMany() will tell us how many
739 break;
740 }
742 //
743 // check locking if we didn't converge or restart
744 //
746 else if (locktest != Teuchos::null && locktest->getStatus() == Passed) {
747
748#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
749 Teuchos::TimeMonitor lcktimer(*_timerLocking);
750#endif
751
752 //
753 // get current state
754 TraceMinBaseState<ScalarType,MV> state = tm_solver->getState();
755 const int curdim = state.curDim;
756
757 //
758 // get number,indices of vectors to be locked
759 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() <= 0,std::logic_error,
760 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() non-positive.");
761 TEUCHOS_TEST_FOR_EXCEPTION(locktest->howMany() != (int)locktest->whichVecs().size(),std::logic_error,
762 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: howMany() not consistent with whichVecs().");
763 // we should have room to lock vectors; otherwise, locking should have been deactivated
764 TEUCHOS_TEST_FOR_EXCEPTION(curNumLocked == maxLocked_,std::logic_error,
765 "Anasazi::TraceMinBaseSolMgr::solve(): status test mistake: locking not deactivated.");
766 //
767 // don't lock more than maxLocked_; we didn't allocate enough space.
768 std::vector<int> tmp_vector_int;
769 if (curNumLocked + locktest->howMany() > maxLocked_) {
770 // just use the first of them
771 for(int i=0; i<maxLocked_-curNumLocked; i++)
772 tmp_vector_int.push_back(locktest->whichVecs()[i]);
773// tmp_vector_int.insert(tmp_vector_int.begin(),locktest->whichVecs().begin(),locktest->whichVecs().begin()+maxLocked_-curNumLocked);
774 }
775 else {
776 tmp_vector_int = locktest->whichVecs();
777 }
778
779 const std::vector<int> lockind(tmp_vector_int);
780 const int numNewLocked = lockind.size();
781 //
782 // generate indices of vectors left unlocked
783 // curind = [0,...,curdim-1] = UNION( lockind, unlockind )
784 const int numUnlocked = curdim-numNewLocked;
785 tmp_vector_int.resize(curdim);
786 for (int i=0; i<curdim; i++) tmp_vector_int[i] = i;
787 const std::vector<int> curind(tmp_vector_int); // curind = [0 ... curdim-1]
788 tmp_vector_int.resize(numUnlocked);
789 std::set_difference(curind.begin(),curind.end(),lockind.begin(),lockind.end(),tmp_vector_int.begin());
790 const std::vector<int> unlockind(tmp_vector_int); // unlockind = [0 ... curdim-1] - lockind
791 tmp_vector_int.clear();
792
793 //
794 // debug printing
795 if (printer_->isVerbosity(Debug)) {
796 printer_->print(Debug,"Locking vectors: ");
797 for (unsigned int i=0; i<lockind.size(); i++) {printer_->stream(Debug) << " " << lockind[i];}
798 printer_->print(Debug,"\n");
799 printer_->print(Debug,"Not locking vectors: ");
800 for (unsigned int i=0; i<unlockind.size(); i++) {printer_->stream(Debug) << " " << unlockind[i];}
801 printer_->print(Debug,"\n");
802 }
803
804 // Copy eigenvalues we want to lock into lockvals
805 std::vector<Value<ScalarType> > allvals = tm_solver->getRitzValues();
806 for(unsigned int i=0; i<allvals.size(); i++)
807 printer_->stream(Debug) << "Ritz value[" << i << "] = " << allvals[i].realpart << std::endl;
808 for (int i=0; i<numNewLocked; i++) {
809 lockvals.push_back(allvals[lockind[i]].realpart);
810 }
811
812 // Copy vectors we want to lock into lockvecs
813 RCP<const MV> newLocked = MVT::CloneView(*tm_solver->getRitzVectors(),lockind);
814 std::vector<int> indlock(numNewLocked);
815 for (int i=0; i<numNewLocked; i++) indlock[i] = curNumLocked+i;
816 if(useHarmonic_)
817 {
818 RCP<MV> tempMV = MVT::CloneCopy(*newLocked);
819 ortho->normalizeMat(*tempMV);
820 MVT::SetBlock(*tempMV,indlock,*lockvecs);
821 }
822 else
823 {
824 MVT::SetBlock(*newLocked,indlock,*lockvecs);
825 }
826
827 // Tell the StatusTestWithOrdering that things have been locked
828 // This is VERY important
829 // If this set of lines is removed, the code does not terminate correctly
830 if(noSort_)
831 {
832 for(unsigned int aliciaInd=0; aliciaInd<lockvals.size(); aliciaInd++)
833 {
834 lockvals[aliciaInd] = 0.0;
835 }
836 }
837 ordertest->setAuxVals(lockvals);
838
839 // Set the auxiliary vectors so that we remain orthogonal to the ones we locked
840 // Remember to include any aux vecs provided by the user
841 curNumLocked += numNewLocked;
842
843 if(ordertest->getStatus() == Passed) break;
844
845 std::vector<int> curlockind(curNumLocked);
846 for (int i=0; i<curNumLocked; i++) curlockind[i] = i;
847 RCP<const MV> curlocked = MVT::CloneView(*lockvecs,curlockind);
848
849 Teuchos::Array< RCP<const MV> > aux;
850 if (probauxvecs != Teuchos::null) aux.push_back(probauxvecs);
851 aux.push_back(curlocked);
852 tm_solver->setAuxVecs(aux);
853
854 // Disable locking if we have locked the maximum number of things
855 printer_->stream(Debug) << "curNumLocked: " << curNumLocked << std::endl;
856 printer_->stream(Debug) << "maxLocked: " << maxLocked_ << std::endl;
857 if (curNumLocked == maxLocked_) {
858 // disabled locking now
859 combotest->removeTest(locktest);
860 locktest = Teuchos::null;
861 printer_->stream(Debug) << "Removed locking test\n";
862 }
863
864 int newdim = numRestartBlocks_*blockSize_;
866 if(newdim <= numUnlocked)
867 {
868 if(useHarmonic_)
869 {
870 std::vector<int> desiredSubscripts(newdim);
871 for(int i=0; i<newdim; i++)
872 {
873 desiredSubscripts[i] = unlockind[i];
874 printer_->stream(Debug) << "H desiredSubscripts[" << i << "] = " << desiredSubscripts[i] << std::endl;
875 }
876 newstate.V = MVT::CloneView(*tm_solver->getRitzVectors(),desiredSubscripts);
877 newstate.curDim = newdim;
878 }
879 else
880 {
881 std::vector<int> desiredSubscripts(newdim);
882 for(int i=0; i<newdim; i++)
883 {
884 desiredSubscripts[i] = unlockind[i];
885 printer_->stream(Debug) << "desiredSubscripts[" << i << "] = " << desiredSubscripts[i] << std::endl;
886 }
887
888 copyPartOfState(state, newstate, desiredSubscripts);
889 }
890 }
891 else
892 {
893 // TODO: Come back to this and make it more efficient
894
895 // Replace the converged eigenvectors with random ones
896 int nrandom = newdim-numUnlocked;
897
898 RCP<const MV> helperMV;
899 RCP<MV> totalV = MVT::Clone(*tm_solver->getRitzVectors(),newdim);
900
901 // Holds old things that we're keeping
902 tmp_vector_int.resize(numUnlocked);
903 for(int i=0; i<numUnlocked; i++) tmp_vector_int[i] = i;
904 RCP<MV> oldV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
905
906 // Copy over the old things
907 if(useHarmonic_)
908 helperMV = MVT::CloneView(*tm_solver->getRitzVectors(),unlockind);
909 else
910 helperMV = MVT::CloneView(*state.V,unlockind);
911 MVT::Assign(*helperMV,*oldV);
912
913 // Holds random vectors we're generating
914 tmp_vector_int.resize(nrandom);
915 for(int i=0; i<nrandom; i++) tmp_vector_int[i] = i+numUnlocked;
916 RCP<MV> newV = MVT::CloneViewNonConst(*totalV,tmp_vector_int);
917
918 // Create random things
919 MVT::MvRandom(*newV);
920
921 newstate.V = totalV;
922 newstate.curDim = newdim;
923
924 // Copy Ritz shifts
925 RCP< std::vector<ScalarType> > helperRS = rcp( new std::vector<ScalarType>(blockSize_) );
926 for(unsigned int i=0; i<(unsigned int)blockSize_; i++)
927 {
928 if(i < unlockind.size() && unlockind[i] < blockSize_)
929 (*helperRS)[i] = (*state.ritzShifts)[unlockind[i]];
930 else
931 (*helperRS)[i] = ZERO;
932 }
933 newstate.ritzShifts = helperRS;
934 }
935
936 // Determine the largest safe shift
937 newstate.largestSafeShift = std::abs(lockvals[0]);
938 for(size_t i=0; i<lockvals.size(); i++)
939 newstate.largestSafeShift = std::max(std::abs(lockvals[i]), newstate.largestSafeShift);
940
941 // Prepare new state, removing converged vectors
942 // TODO: Init will perform some unnecessary calculations; look into it
943 // TODO: The residual norms should be part of the state
944 newstate.NEV = state.NEV - numNewLocked;
945 tm_solver->initialize(newstate);
946 } // end of locking
948 //
949 // check for restarting before locking: if we need to lock, it will happen after the restart
950 //
952 else if ( needToRestart(tm_solver) ) {
953 // if performRestart returns false, we exceeded the maximum number of restarts
954 if(performRestart(numRestarts, tm_solver) == false)
955 break;
956 } // end of restarting
958 //
959 // we returned from iterate(), but none of our status tests Passed.
960 // something is wrong, and it is probably our fault.
961 //
963 else {
964 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): Invalid return from tm_solver::iterate().");
965 }
966 }
967 catch (const AnasaziError &err) {
968 printer_->stream(Errors)
969 << "Anasazi::TraceMinBaseSolMgr::solve() caught unexpected exception from Anasazi::TraceMinBase::iterate() at iteration " << tm_solver->getNumIters() << std::endl
970 << err.what() << std::endl
971 << "Anasazi::TraceMinBaseSolMgr::solve() returning Unconverged with no solutions." << std::endl;
972 return Unconverged;
973 }
974 }
975
976 sol.numVecs = ordertest->howMany();
977 if (sol.numVecs > 0) {
978 sol.Evecs = MVT::Clone(*problem_->getInitVec(),sol.numVecs);
979 sol.Espace = sol.Evecs;
980 sol.Evals.resize(sol.numVecs);
981 std::vector<MagnitudeType> vals(sol.numVecs);
982
983 // copy them into the solution
984 std::vector<int> which = ordertest->whichVecs();
985 // indices between [0,blockSize) refer to vectors/values in the solver
986 // indices between [-curNumLocked,-1] refer to locked vectors/values [0,curNumLocked)
987 // everything has already been ordered by the solver; we just have to partition the two references
988 std::vector<int> inlocked(0), insolver(0);
989 for (unsigned int i=0; i<which.size(); i++) {
990 if (which[i] >= 0) {
991 TEUCHOS_TEST_FOR_EXCEPTION(which[i] >= blockSize_,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): positive indexing mistake from ordertest.");
992 insolver.push_back(which[i]);
993 }
994 else {
995 // sanity check
996 TEUCHOS_TEST_FOR_EXCEPTION(which[i] < -curNumLocked,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): negative indexing mistake from ordertest.");
997 inlocked.push_back(which[i] + curNumLocked);
998 }
999 }
1000
1001 TEUCHOS_TEST_FOR_EXCEPTION(insolver.size() + inlocked.size() != (unsigned int)sol.numVecs,std::logic_error,"Anasazi::TraceMinBaseSolMgr::solve(): indexing mistake.");
1002
1003 // set the vecs,vals in the solution
1004 if (insolver.size() > 0) {
1005 // set vecs
1006 int lclnum = insolver.size();
1007 std::vector<int> tosol(lclnum);
1008 for (int i=0; i<lclnum; i++) tosol[i] = i;
1009 RCP<const MV> v = MVT::CloneView(*tm_solver->getRitzVectors(),insolver);
1010 MVT::SetBlock(*v,tosol,*sol.Evecs);
1011 // set vals
1012 std::vector<Value<ScalarType> > fromsolver = tm_solver->getRitzValues();
1013 for (unsigned int i=0; i<insolver.size(); i++) {
1014 vals[i] = fromsolver[insolver[i]].realpart;
1015 }
1016 }
1017
1018 // get the vecs,vals from locked storage
1019 if (inlocked.size() > 0) {
1020 int solnum = insolver.size();
1021 // set vecs
1022 int lclnum = inlocked.size();
1023 std::vector<int> tosol(lclnum);
1024 for (int i=0; i<lclnum; i++) tosol[i] = solnum + i;
1025 RCP<const MV> v = MVT::CloneView(*lockvecs,inlocked);
1026 MVT::SetBlock(*v,tosol,*sol.Evecs);
1027 // set vals
1028 for (unsigned int i=0; i<inlocked.size(); i++) {
1029 vals[i+solnum] = lockvals[inlocked[i]];
1030 }
1031 }
1032
1033 // undo the spectral transformation if necessary
1034 // if we really passed the solver Bx = \lambda A x, invert the eigenvalues
1035 if(which_ == "LM")
1036 {
1037 for(size_t i=0; i<vals.size(); i++)
1038 vals[i] = ONE/vals[i];
1039 }
1040
1041 // sort the eigenvalues and permute the eigenvectors appropriately
1042 {
1043 std::vector<int> order(sol.numVecs);
1044 sorter->sort(vals,Teuchos::rcpFromRef(order),sol.numVecs);
1045 // store the values in the Eigensolution
1046 for (int i=0; i<sol.numVecs; i++) {
1047 sol.Evals[i].realpart = vals[i];
1048 sol.Evals[i].imagpart = MT::zero();
1049 }
1050 // now permute the eigenvectors according to order
1051 msutils::permuteVectors(sol.numVecs,order,*sol.Evecs);
1052 }
1053
1054 // setup sol.index, remembering that all eigenvalues are real so that index = {0,...,0}
1055 sol.index.resize(sol.numVecs,0);
1056 }
1057 }
1058
1059 // print final summary
1060 tm_solver->currentStatus(printer_->stream(FinalSummary));
1061
1062 printParameters(printer_->stream(FinalSummary));
1063
1064 // print timing information
1065#ifdef ANASAZI_TEUCHOS_TIME_MONITOR
1066 if ( printer_->isVerbosity( TimingDetails ) ) {
1067 Teuchos::TimeMonitor::summarize( printer_->stream( TimingDetails ) );
1068 }
1069#endif
1070
1071 problem_->setSolution(sol);
1072 printer_->stream(Debug) << "Returning " << sol.numVecs << " eigenpairs to eigenproblem." << std::endl;
1073
1074 // get the number of iterations taken for this call to solve().
1075 numIters_ = tm_solver->getNumIters();
1076
1077 if (sol.numVecs < nev) {
1078 return Unconverged; // return from TraceMinBaseSolMgr::solve()
1079 }
1080 return Converged; // return from TraceMinBaseSolMgr::solve()
1081}
1082
1083
1084template <class ScalarType, class MV, class OP>
1085void
1087 const RCP< StatusTest<ScalarType,MV,OP> > &global)
1088{
1089 globalTest_ = global;
1090}
1091
1092template <class ScalarType, class MV, class OP>
1093const RCP< StatusTest<ScalarType,MV,OP> > &
1095{
1096 return globalTest_;
1097}
1098
1099template <class ScalarType, class MV, class OP>
1100void
1102 const RCP< StatusTest<ScalarType,MV,OP> > &debug)
1103{
1104 debugTest_ = debug;
1105}
1106
1107template <class ScalarType, class MV, class OP>
1108const RCP< StatusTest<ScalarType,MV,OP> > &
1110{
1111 return debugTest_;
1112}
1113
1114template <class ScalarType, class MV, class OP>
1115void
1117 const RCP< StatusTest<ScalarType,MV,OP> > &locking)
1118{
1119 lockingTest_ = locking;
1120}
1121
1122template <class ScalarType, class MV, class OP>
1123const RCP< StatusTest<ScalarType,MV,OP> > &
1125{
1126 return lockingTest_;
1127}
1128
1129template <class ScalarType, class MV, class OP>
1130void TraceMinBaseSolMgr<ScalarType,MV,OP>::copyPartOfState(const TraceMinBaseState<ScalarType,MV>& oldState, TraceMinBaseState<ScalarType,MV>& newState, const std::vector<int> indToCopy) const
1131{
1132 const ScalarType ONE = Teuchos::ScalarTraits<MagnitudeType>::one();
1133 const ScalarType ZERO = Teuchos::ScalarTraits<MagnitudeType>::zero();
1134
1135 newState.curDim = indToCopy.size();
1136 std::vector<int> fullIndices(oldState.curDim);
1137 for(int i=0; i<oldState.curDim; i++) fullIndices[i] = i;
1138
1139 // Initialize with X.
1140 // Note that we didn't compute enough vectors of X, but we can very easily using the Ritz vectors.
1141 // That's why they're part of the state.
1142 // Note that there will ALWAYS be enough vectors
1143
1144 // Helpful vectors for computing views and whatnot
1145 std::vector<int> oldIndices;
1146 std::vector<int> newIndices;
1147 for(int i=0; i<newState.curDim; i++)
1148 {
1149 if(indToCopy[i] < blockSize_)
1150 oldIndices.push_back(indToCopy[i]);
1151 else
1152 newIndices.push_back(indToCopy[i]);
1153 }
1154
1155 int olddim = oldIndices.size();
1156 int newdim = newIndices.size();
1157
1158 // If there are no new vectors being copied
1159 if(computeAllRes_)
1160 {
1161 newState.V = MVT::CloneView(*oldState.X, indToCopy);
1162 newState.R = MVT::CloneView(*oldState.R, indToCopy);
1163 newState.X = newState.V;
1164
1165 if(problem_->getOperator() != Teuchos::null)
1166 {
1167 newState.KV = MVT::CloneView(*oldState.KX, indToCopy);
1168 newState.KX = newState.KV;
1169 }
1170 else
1171 {
1172 newState.KV = Teuchos::null;
1173 newState.KX = Teuchos::null;
1174 }
1175
1176 if(problem_->getM() != Teuchos::null)
1177 {
1178 newState.MopV = MVT::CloneView(*oldState.MX, indToCopy);
1179 newState.MX = newState.MopV;
1180 }
1181 else
1182 {
1183 newState.MopV = Teuchos::null;
1184 newState.MX = Teuchos::null;
1185 }
1186 }
1187 else if(newdim == 0)
1188 {
1189 std::vector<int> blockind(blockSize_);
1190 for(int i=0; i<blockSize_; i++)
1191 blockind[i] = i;
1192
1193 // Initialize with X
1194 newState.V = MVT::CloneView(*oldState.X, blockind);
1195 newState.KV = MVT::CloneView(*oldState.KX, blockind);
1196 newState.R = MVT::CloneView(*oldState.R, blockind);
1197 newState.X = MVT::CloneView(*newState.V, blockind);
1198 newState.KX = MVT::CloneView(*newState.KV, blockind);
1199
1200 if(problem_->getM() != Teuchos::null)
1201 {
1202 newState.MopV = MVT::CloneView(*oldState.MX, blockind);
1203 newState.MX = MVT::CloneView(*newState.MopV, blockind);
1204 }
1205 else
1206 {
1207 newState.MopV = Teuchos::null;
1208 newState.MX = Teuchos::null;
1209 }
1210 }
1211 else
1212 {
1213 // More helpful vectors
1214 std::vector<int> oldPart(olddim);
1215 for(int i=0; i<olddim; i++) oldPart[i] = i;
1216 std::vector<int> newPart(newdim);
1217 for(int i=0; i<newdim; i++) newPart[i] = olddim+i;
1218
1219 // Helpful multivectors for views and whatnot
1220 RCP<MV> helper = MVT::Clone(*oldState.V,newState.curDim);
1221 RCP<MV> oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1222 RCP<MV> newHelper = MVT::CloneViewNonConst(*helper,newPart);
1223 RCP<const MV> viewHelper;
1224
1225 // Get the parts of the Ritz vectors we are interested in.
1226 Teuchos::SerialDenseMatrix<int,ScalarType> newRV(oldState.curDim,newdim);
1227 for(int r=0; r<oldState.curDim; r++)
1228 {
1229 for(int c=0; c<newdim; c++)
1230 newRV(r,c) = (*oldState.RV)(r,newIndices[c]);
1231 }
1232
1233 // We're going to compute X as V*RitzVecs
1234 viewHelper = MVT::CloneView(*oldState.V,fullIndices);
1235 MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1236 viewHelper = MVT::CloneView(*oldState.X,oldIndices);
1237 MVT::Assign(*viewHelper,*oldHelper);
1238 newState.V = MVT::CloneCopy(*helper);
1239
1240 // Also compute KX as KV*RitzVecs
1241 viewHelper = MVT::CloneView(*oldState.KV,fullIndices);
1242 MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1243 viewHelper = MVT::CloneView(*oldState.KX,oldIndices);
1244 MVT::Assign(*viewHelper,*oldHelper);
1245 newState.KV = MVT::CloneCopy(*helper);
1246
1247 // Do the same with MX if necessary
1248 if(problem_->getM() != Teuchos::null)
1249 {
1250 viewHelper = MVT::CloneView(*oldState.MopV,fullIndices);
1251 MVT::MvTimesMatAddMv(ONE,*viewHelper,newRV,ZERO,*newHelper);
1252 viewHelper = MVT::CloneView(*oldState.MX,oldIndices);
1253 MVT::Assign(*viewHelper,*oldHelper);
1254 newState.MopV = MVT::CloneCopy(*helper);
1255 }
1256 else
1257 newState.MopV = newState.V;
1258
1259 // Get X, MX, KX
1260 std::vector<int> blockVec(blockSize_);
1261 for(int i=0; i<blockSize_; i++) blockVec[i] = i;
1262 newState.X = MVT::CloneView(*newState.V,blockVec);
1263 newState.KX = MVT::CloneView(*newState.KV,blockVec);
1264 newState.MX = MVT::CloneView(*newState.MopV,blockVec);
1265
1266 // Update the residuals
1267 if(blockSize_-oldIndices.size() > 0)
1268 {
1269 // There are vectors we have not computed the residual for yet
1270 newPart.resize(blockSize_-oldIndices.size());
1271 helper = MVT::Clone(*oldState.V,blockSize_);
1272 oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1273 newHelper = MVT::CloneViewNonConst(*helper,newPart);
1274
1275 RCP<MV> scaledMV = MVT::CloneCopy(*newState.MX,newPart);
1276 RCP<const MV> localKV = MVT::CloneView(*newState.KX,newPart);
1277 std::vector<ScalarType> scalarVec(blockSize_-oldIndices.size());
1278 for(unsigned int i=0; i<(unsigned int)blockSize_-oldIndices.size(); i++) scalarVec[i] = (*oldState.T)[newPart[i]];
1279 MVT::MvScale(*scaledMV,scalarVec);
1280
1281 helper = MVT::Clone(*oldState.V,blockSize_);
1282 oldHelper = MVT::CloneViewNonConst(*helper,oldPart);
1283 newHelper = MVT::CloneViewNonConst(*helper,newPart);
1284 MVT::MvAddMv(ONE,*localKV,-ONE,*scaledMV,*newHelper);
1285 viewHelper = MVT::CloneView(*oldState.R,oldIndices);
1286 MVT::Assign(*viewHelper,*oldHelper);
1287 newState.R = MVT::CloneCopy(*helper);
1288 }
1289 else
1290 newState.R = oldState.R;
1291 }
1292
1293 // Since we are setting V:=X, V is orthonormal
1294 newState.isOrtho = true;
1295
1296 // Get the first eigenvalues
1297 RCP< std::vector<ScalarType> > helperT = rcp( new std::vector<ScalarType>(newState.curDim) );
1298 for(int i=0; i<newState.curDim; i++) (*helperT)[i] = (*oldState.T)[indToCopy[i]];
1299 newState.T = helperT;
1300
1301 // X'KX is diag(T)
1302 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > newKK = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newState.curDim,newState.curDim) );
1303 for(int i=0; i<newState.curDim; i++)
1304 (*newKK)(i,i) = (*newState.T)[i];
1305 newState.KK = newKK;
1306
1307 // The associated Ritz vectors are I
1308 RCP< Teuchos::SerialDenseMatrix<int,ScalarType> > newRV = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(newState.curDim,newState.curDim) );
1309 for(int i=0; i<newState.curDim; i++)
1310 (*newRV)(i,i) = ONE;
1311 newState.RV = newRV;
1312
1313 // Get the Ritz shifts
1314 RCP< std::vector<ScalarType> > helperRS = rcp( new std::vector<ScalarType>(blockSize_) );
1315 for(int i=0; i<blockSize_; i++)
1316 {
1317 if(indToCopy[i] < blockSize_)
1318 (*helperRS)[i] = (*oldState.ritzShifts)[indToCopy[i]];
1319 else
1320 (*helperRS)[i] = ZERO;
1321 }
1322 newState.ritzShifts = helperRS;
1323}
1324
1325
1326template <class ScalarType, class MV, class OP>
1327void TraceMinBaseSolMgr<ScalarType,MV,OP>::setParameters(Teuchos::ParameterList &pl) const
1328{
1329 pl.set("Block Size", blockSize_);
1330 pl.set("Num Blocks", numBlocks_);
1331 pl.set("Num Restart Blocks", numRestartBlocks_);
1332 pl.set("When To Shift", whenToShift_);
1333 pl.set("Trace Threshold", traceThresh_);
1334 pl.set("Shift Tolerance", shiftTol_);
1335 pl.set("Relative Shift Tolerance", relShiftTol_);
1336 pl.set("Shift Norm", shiftNorm_);
1337 pl.set("How To Choose Shift", howToShift_);
1338 pl.set("Consider Clusters", considerClusters_);
1339 pl.set("Use Multiple Shifts", useMultipleShifts_);
1340 pl.set("Saddle Solver Type", saddleSolType_);
1341 pl.set("Project All Vectors", projectAllVecs_);
1342 pl.set("Project Locked Vectors", projectLockedVecs_);
1343 pl.set("Compute All Residuals", computeAllRes_);
1344 pl.set("Use Residual as RHS", useRHSR_);
1345 pl.set("Use Harmonic Ritz Values", useHarmonic_);
1346 pl.set("Maximum Krylov Iterations", maxKrylovIter_);
1347 pl.set("HSS: alpha", alpha_);
1348}
1349
1350
1351template <class ScalarType, class MV, class OP>
1352void TraceMinBaseSolMgr<ScalarType,MV,OP>::printParameters(std::ostream &os) const
1353{
1354 os << "\n\n\n";
1355 os << "========================================\n";
1356 os << "========= TraceMin parameters ==========\n";
1357 os << "========================================\n";
1358 os << "=========== Block parameters ===========\n";
1359 os << "Block Size: " << blockSize_ << std::endl;
1360 os << "Num Blocks: " << numBlocks_ << std::endl;
1361 os << "Num Restart Blocks: " << numRestartBlocks_ << std::endl;
1362 os << "======== Convergence parameters ========\n";
1363 os << "Convergence Tolerance: " << convTol_ << std::endl;
1364 os << "Relative Convergence Tolerance: " << relConvTol_ << std::endl;
1365 os << "========== Locking parameters ==========\n";
1366 os << "Use Locking: " << useLocking_ << std::endl;
1367 os << "Locking Tolerance: " << lockTol_ << std::endl;
1368 os << "Relative Locking Tolerance: " << relLockTol_ << std::endl;
1369 os << "Max Locked: " << maxLocked_ << std::endl;
1370 os << "Locking Quorum: " << lockQuorum_ << std::endl;
1371 os << "========== Shifting parameters =========\n";
1372 os << "When To Shift: ";
1373 if(whenToShift_ == NEVER_SHIFT) os << "Never\n";
1374 else if(whenToShift_ == SHIFT_WHEN_TRACE_LEVELS) os << "After Trace Levels\n";
1375 else if(whenToShift_ == SHIFT_WHEN_RESID_SMALL) os << "Residual Becomes Small\n";
1376 else if(whenToShift_ == ALWAYS_SHIFT) os << "Always\n";
1377 os << "Consider Clusters: " << considerClusters_ << std::endl;
1378 os << "Trace Threshohld: " << traceThresh_ << std::endl;
1379 os << "Shift Tolerance: " << shiftTol_ << std::endl;
1380 os << "Relative Shift Tolerance: " << relShiftTol_ << std::endl;
1381 os << "How To Choose Shift: ";
1382 if(howToShift_ == LARGEST_CONVERGED_SHIFT) os << "Largest Converged\n";
1383 else if(howToShift_ == ADJUSTED_RITZ_SHIFT) os << "Adjusted Ritz Values\n";
1384 else if(howToShift_ == RITZ_VALUES_SHIFT) os << "Ritz Values\n";
1385 os << "Use Multiple Shifts: " << useMultipleShifts_ << std::endl;
1386 os << "=========== Other parameters ===========\n";
1387 os << "Orthogonalization: " << ortho_ << std::endl;
1388 os << "Saddle Solver Type: ";
1389 if(saddleSolType_ == PROJECTED_KRYLOV_SOLVER) os << "Projected Krylov\n";
1390 else if(saddleSolType_ == SCHUR_COMPLEMENT_SOLVER) os << "Schur Complement\n";
1391 os << "Project All Vectors: " << projectAllVecs_ << std::endl;
1392 os << "Project Locked Vectors: " << projectLockedVecs_ << std::endl;
1393 os << "Compute All Residuals: " << computeAllRes_ << std::endl;
1394 os << "========================================\n\n\n";
1395}
1396
1397
1398}} // end Anasazi namespace
1399
1400#endif /* ANASAZI_TraceMinBase_SOLMGR_HPP */
Basic implementation of the Anasazi::OrthoManager class.
Basic implementation of the Anasazi::SortManager class.
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...
Basic implementation of the Anasazi::OrthoManager class.
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...
Abstract base class for trace minimization eigensolvers.
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...
This class defines the interface required by an eigensolver and status test class to compute solution...
The Anasazi::TraceMinBaseSolMgr provides an abstract base class for the TraceMin series of solver man...
const RCP< StatusTest< ScalarType, MV, OP > > & getLockingStatusTest() const
Get the status test defining locking.
Teuchos::Array< RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
const RCP< StatusTest< ScalarType, MV, OP > > & getDebugStatusTest() const
Get the status test for debugging.
ReturnType solve()
This method performs possibly repeated calls to the underlying eigensolver's iterate() routine until ...
void setDebugStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &debug)
Set the status test for debugging.
void setGlobalStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &global)
Set the status test defining global convergence.
TraceMinBaseSolMgr(const RCP< Eigenproblem< ScalarType, MV, OP > > &problem, Teuchos::ParameterList &pl)
Basic constructor for TraceMinBaseSolMgr.
int getNumIters() const
Get the iteration count for the most recent call to solve().
void setLockingStatusTest(const RCP< StatusTest< ScalarType, MV, OP > > &locking)
Set the status test defining locking.
const RCP< StatusTest< ScalarType, MV, OP > > & getGlobalStatusTest() const
Get the status test defining global convergence.
const Eigenproblem< ScalarType, MV, OP > & getProblem() const
Return the eigenvalue problem.
This is an abstract base class for the trace minimization eigensolvers.
An implementation of the Anasazi::GenOrthoManager that performs orthogonalization using iterated clas...
Anasazi's templated virtual class for providing routines for orthogonalization and orthonormalization...
Traits class which defines basic operations on multivectors.
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.
Anasazi's templated pure virtual class for managing the sorting of approximate eigenvalues computed b...
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.
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.
Structure to contain pointers to TraceMinBase state variables.
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > KK
The current projected K matrix.
RCP< const MV > MX
The image of the current eigenvectors under M, or Teuchos::null if M was not specified.
RCP< const MV > KX
The image of the current eigenvectors under K.
ScalarType largestSafeShift
Largest safe shift.
int curDim
The current dimension of the solver.
RCP< const MV > MopV
The image of V under M, or Teuchos::null if M was not specified.
RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > RV
The current Ritz vectors.
RCP< const MV > V
The current basis.
bool isOrtho
Whether V has been projected and orthonormalized already.
RCP< const std::vector< typename Teuchos::ScalarTraits< ScalarType >::magnitudeType > > T
The current Ritz values. This vector is a copy of the internal data.
RCP< const MV > X
The current eigenvectors.
RCP< const MV > R
The current residual vectors.
int NEV
Number of unconverged eigenvalues.
RCP< const MV > KV
The image of V under K.
RCP< const std::vector< ScalarType > > ritzShifts
Current Ritz shifts.
Output managers remove the need for the eigensolver to know any information about the required output...