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