Belos Version of the Day
Loading...
Searching...
No Matches
BelosRCGSolMgr.hpp
Go to the documentation of this file.
1// @HEADER
2// *****************************************************************************
3// Belos: Block Linear Solvers Package
4//
5// Copyright 2004-2016 NTESS and the Belos contributors.
6// SPDX-License-Identifier: BSD-3-Clause
7// *****************************************************************************
8// @HEADER
9
10#ifndef BELOS_RCG_SOLMGR_HPP
11#define BELOS_RCG_SOLMGR_HPP
12
16
17#include "BelosConfigDefs.hpp"
18#include "BelosTypes.hpp"
19
20#include "BelosLinearProblem.hpp"
21#include "BelosSolverManager.hpp"
22
23#include "BelosRCGIter.hpp"
24#include "BelosStatusTestMaxIters.hpp"
25#include "BelosStatusTestGenResNorm.hpp"
26#include "BelosStatusTestCombo.hpp"
27#include "BelosStatusTestOutputFactory.hpp"
28#include "BelosOutputManager.hpp"
29#include "Teuchos_BLAS.hpp"
30#include "Teuchos_LAPACK.hpp"
31#include "Teuchos_as.hpp"
32#ifdef BELOS_TEUCHOS_TIME_MONITOR
33#include "Teuchos_TimeMonitor.hpp"
34#endif
35
42
81
82namespace Belos {
83
85
86
93 class RCGSolMgrLinearProblemFailure : public BelosError {public:
94 RCGSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg)
95 {}};
96
103 class RCGSolMgrLAPACKFailure : public BelosError {public:
104 RCGSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
105 {}};
106
108
109
110 // Partial specialization for unsupported ScalarType types.
111 // This contains a stub implementation.
112 template<class ScalarType, class MV, class OP, class DM = DefaultDenseMatrix<int,ScalarType>,
113 const bool supportsScalarType =
114 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
115 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
116 class RCGSolMgr :
117 public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, DM,
118 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
119 ! Teuchos::ScalarTraits<ScalarType>::isComplex>
120 {
121 static const bool scalarTypeIsSupported =
122 Belos::Details::LapackSupportsScalar<ScalarType>::value &&
123 ! Teuchos::ScalarTraits<ScalarType>::isComplex;
124 typedef Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, DM,
125 scalarTypeIsSupported> base_type;
126
127 public:
128 RCGSolMgr () :
129 base_type ()
130 {}
131 RCGSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP,DM> > &problem,
132 const Teuchos::RCP<Teuchos::ParameterList> &pl) :
133 base_type ()
134 {}
135 virtual ~RCGSolMgr () {}
136
138 Teuchos::RCP<SolverManager<ScalarType, MV, OP, DM> > clone () const override {
139 return Teuchos::rcp(new RCGSolMgr<ScalarType,MV,OP,DM,supportsScalarType>);
140 }
141 };
142
143 // Partial specialization for real ScalarType.
144 // This contains the actual working implementation of RCG.
145 // See discussion in the class documentation above.
146 template<class ScalarType, class MV, class OP, class DM>
147 class RCGSolMgr<ScalarType, MV, OP, DM, true> :
148 public Details::SolverManagerRequiresRealLapack<ScalarType, MV, OP, DM, true> {
149 private:
150 typedef MultiVecTraits<ScalarType,MV,DM> MVT;
151 typedef OperatorTraits<ScalarType,MV,OP> OPT;
152 typedef DenseMatTraits<ScalarType,DM> DMT;
153 typedef Teuchos::ScalarTraits<ScalarType> SCT;
154 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
155 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
156
157 public:
158
160
161
167 RCGSolMgr();
168
189
190 RCGSolMgr( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP,DM> > &problem,
191 const Teuchos::RCP<Teuchos::ParameterList> &pl );
192
194 virtual ~RCGSolMgr() {};
195
197 Teuchos::RCP<SolverManager<ScalarType, MV, OP, DM> > clone () const override {
198 return Teuchos::rcp(new RCGSolMgr<ScalarType,MV,OP,DM>);
199 }
200
201
203
204
205 const LinearProblem<ScalarType,MV,OP,DM>& getProblem() const override {
206 return *problem_;
207 }
208
210 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
211
213 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override { return params_; }
214
220 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
221 return Teuchos::tuple(timerSolve_);
222 }
223
228 MagnitudeType achievedTol() const override {
229 return achievedTol_;
230 }
231
233 int getNumIters() const override {
234 return numIters_;
235 }
236
238 bool isLOADetected() const override { return false; }
239
241
243
244
246 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP,DM> > &problem ) override { problem_ = problem; }
247
249 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
250
252 void setDebugStatusTest( const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &debugStatusTest ) override {
253 debugStatusTest_ = debugStatusTest;
254 // Force the cached status-test tree (and its output wrapper) to be rebuilt
255 // on the next solve so the debug test is wired into sTest_. This manager
256 // caches sTest_/outputTest_ behind null checks, so clear them and re-run
257 // setParameters.
258 sTest_ = Teuchos::null;
259 outputTest_ = Teuchos::null;
260 params_Set_ = false;
261 }
262
264
266
267
273 void reset( const ResetType type ) override {
274 if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) problem_->setProblem();
275 else if (type & Belos::RecycleSubspace) existU_ = false;
276 }
277
278
280
281
299 ReturnType solve() override;
300
302
305
307 std::string description() const override;
308
310
311 private:
312
313 // Called by all constructors; Contains init instructions common to all constructors
314 void init();
315
316 // Computes harmonic eigenpairs of projected matrix created during one cycle.
317 // Y contains the harmonic Ritz vectors corresponding to the recycleBlocks eigenvalues of smallest magnitude.
318 void getHarmonicVecs(const DM& F,
319 const DM& G,
320 DM& Y);
321
322 // Sort list of n floating-point numbers and return permutation vector
323 void sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm);
324
325 // Initialize solver state storage
326 void initializeStateStorage();
327
328 // Linear problem.
329 Teuchos::RCP<LinearProblem<ScalarType,MV,OP,DM> > problem_;
330
331 // Output manager.
332 Teuchos::RCP<OutputManager<ScalarType> > printer_;
333 Teuchos::RCP<std::ostream> outputStream_;
334
335 // Status test.
336 Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > sTest_;
337 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP,DM> > maxIterTest_;
338 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP,DM> > convTest_;
339 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP,DM> > outputTest_;
340 Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > debugStatusTest_;
341
342 // Current parameter list.
343 Teuchos::RCP<Teuchos::ParameterList> params_;
344
345 // Default solver values.
346 static constexpr int maxIters_default_ = 1000;
347 static constexpr int blockSize_default_ = 1;
348 static constexpr int numBlocks_default_ = 25;
349 static constexpr int recycleBlocks_default_ = 3;
350 static constexpr bool showMaxResNormOnly_default_ = false;
351 static constexpr int verbosity_default_ = Belos::Errors;
352 static constexpr int outputStyle_default_ = Belos::General;
353 static constexpr int outputFreq_default_ = -1;
354 static constexpr const char * label_default_ = "Belos";
355
356 //
357 // Current solver values.
358 //
359
361 MagnitudeType convtol_;
362
367 MagnitudeType achievedTol_;
368
370 int maxIters_;
371
373 int numIters_;
374
375 int numBlocks_, recycleBlocks_;
376 bool showMaxResNormOnly_;
377 int verbosity_, outputStyle_, outputFreq_;
378
380 // Solver State Storage
382 // Search vectors
383 Teuchos::RCP<MV> P_;
384 //
385 // A times current search direction
386 Teuchos::RCP<MV> Ap_;
387 //
388 // Residual vector
389 Teuchos::RCP<MV> r_;
390 //
391 // Preconditioned residual
392 Teuchos::RCP<MV> z_;
393 //
394 // Flag indicating that the recycle space should be used
395 bool existU_;
396 //
397 // Flag indicating that the updated recycle space has been created
398 bool existU1_;
399 //
400 // Recycled subspace and its image
401 Teuchos::RCP<MV> U_, AU_;
402 //
403 // Recycled subspace for next system and its image
404 Teuchos::RCP<MV> U1_;
405 //
406 // Coefficients arising in RCG iteration
407 Teuchos::RCP<std::vector<ScalarType> > Beta_;
408 Teuchos::RCP<std::vector<ScalarType> > Alpha_;
409 Teuchos::RCP<std::vector<ScalarType> > D_;
410 //
411 // Solutions to local least-squares problems
412 Teuchos::RCP<DM> Delta_;
413 //
414 // The matrix U^T A U
415 Teuchos::RCP<DM> UTAU_;
416 //
417 // LU factorization of U^T A U
418 Teuchos::RCP<DM> LUUTAU_;
419 //
420 // Data from LU factorization of UTAU
421 Teuchos::RCP<std::vector<int> > ipiv_;
422 //
423 // The matrix (AU)^T AU
424 Teuchos::RCP<DM> AUTAU_;
425 //
426 // The scalar r'*z
427 Teuchos::RCP<std::vector<ScalarType> > rTz_old_;
428 //
429 // Matrices needed for calculation of harmonic Ritz eigenproblem
430 Teuchos::RCP<DM> F_,G_,Y_;
431 //
432 // Matrices needed for updating recycle space
433 Teuchos::RCP<DM> L2_,DeltaL2_,AU1TUDeltaL2_;
434 Teuchos::RCP<DM> AU1TAU1_, AU1TU1_, AU1TAP_;
435 Teuchos::RCP<DM> FY_,GY_;
436 Teuchos::RCP<DM> APTAP_;
437 Teuchos::RCP<MV> U1Y1_, PY2_;
438 Teuchos::RCP<DM> AUTAP_, AU1TU_;
439 ScalarType dold;
441
442 // Timers.
443 std::string label_;
444 Teuchos::RCP<Teuchos::Time> timerSolve_;
445
446 // Internal state variables.
447 bool params_Set_;
448 };
449
450
451// Empty Constructor
452template<class ScalarType, class MV, class OP, class DM>
453RCGSolMgr<ScalarType,MV,OP,DM,true>::RCGSolMgr():
454 achievedTol_(0.0),
455 numIters_(0)
456{
457 init();
458}
459
460// Basic Constructor
461template<class ScalarType, class MV, class OP, class DM>
462RCGSolMgr<ScalarType,MV,OP,DM,true>::RCGSolMgr(
463 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP,DM> > &problem,
464 const Teuchos::RCP<Teuchos::ParameterList> &pl ) :
465 problem_(problem),
466 achievedTol_(0.0),
467 numIters_(0)
468{
469 init();
470 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null, std::invalid_argument, "Problem not given to solver manager.");
471
472 // If the parameter list pointer is null, then set the current parameters to the default parameter list.
473 if ( !is_null(pl) ) {
474 setParameters( pl );
475 }
476}
477
478// Common instructions executed in all constructors
479template<class ScalarType, class MV, class OP, class DM>
480void RCGSolMgr<ScalarType,MV,OP,DM,true>::init()
481{
482 outputStream_ = Teuchos::rcpFromRef(std::cout);
483 convtol_ = DefaultSolverParameters::convTol;
484 maxIters_ = maxIters_default_;
485 numBlocks_ = numBlocks_default_;
486 recycleBlocks_ = recycleBlocks_default_;
487 verbosity_ = verbosity_default_;
488 outputStyle_= outputStyle_default_;
489 outputFreq_= outputFreq_default_;
490 showMaxResNormOnly_ = showMaxResNormOnly_default_;
491 label_ = label_default_;
492 params_Set_ = false;
493 P_ = Teuchos::null;
494 Ap_ = Teuchos::null;
495 r_ = Teuchos::null;
496 z_ = Teuchos::null;
497 existU_ = false;
498 existU1_ = false;
499 U_ = Teuchos::null;
500 AU_ = Teuchos::null;
501 U1_ = Teuchos::null;
502 Alpha_ = Teuchos::null;
503 Beta_ = Teuchos::null;
504 D_ = Teuchos::null;
505 Delta_ = Teuchos::null;
506 UTAU_ = Teuchos::null;
507 LUUTAU_ = Teuchos::null;
508 ipiv_ = Teuchos::null;
509 AUTAU_ = Teuchos::null;
510 rTz_old_ = Teuchos::null;
511 F_ = Teuchos::null;
512 G_ = Teuchos::null;
513 Y_ = Teuchos::null;
514 L2_ = Teuchos::null;
515 DeltaL2_ = Teuchos::null;
516 AU1TUDeltaL2_ = Teuchos::null;
517 AU1TAU1_ = Teuchos::null;
518 AU1TU1_ = Teuchos::null;
519 AU1TAP_ = Teuchos::null;
520 FY_ = Teuchos::null;
521 GY_ = Teuchos::null;
522 APTAP_ = Teuchos::null;
523 U1Y1_ = Teuchos::null;
524 PY2_ = Teuchos::null;
525 AUTAP_ = Teuchos::null;
526 AU1TU_ = Teuchos::null;
527 dold = 0.;
528}
529
530template<class ScalarType, class MV, class OP, class DM>
531void RCGSolMgr<ScalarType,MV,OP,DM,true>::setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params )
532{
533 // Create the internal parameter list if ones doesn't already exist.
534 if (params_ == Teuchos::null) {
535 params_ = Teuchos::rcp( new Teuchos::ParameterList(*getValidParameters()) );
536 }
537 else {
538 params->validateParameters(*getValidParameters());
539 }
540
541 // Check for maximum number of iterations
542 if (params->isParameter("Maximum Iterations")) {
543 maxIters_ = params->get("Maximum Iterations",maxIters_default_);
544
545 // Update parameter in our list and in status test.
546 params_->set("Maximum Iterations", maxIters_);
547 if (maxIterTest_!=Teuchos::null)
548 maxIterTest_->setMaxIters( maxIters_ );
549 }
550
551 // Check for the maximum number of blocks.
552 if (params->isParameter("Num Blocks")) {
553 numBlocks_ = params->get("Num Blocks",numBlocks_default_);
554 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
555 "Belos::RCGSolMgr: \"Num Blocks\" must be strictly positive.");
556
557 // Update parameter in our list.
558 params_->set("Num Blocks", numBlocks_);
559 }
560
561 // Check for the maximum number of blocks.
562 if (params->isParameter("Num Recycled Blocks")) {
563 recycleBlocks_ = params->get("Num Recycled Blocks",recycleBlocks_default_);
564 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ <= 0, std::invalid_argument,
565 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be strictly positive.");
566
567 TEUCHOS_TEST_FOR_EXCEPTION(recycleBlocks_ >= numBlocks_, std::invalid_argument,
568 "Belos::RCGSolMgr: \"Num Recycled Blocks\" must be less than \"Num Blocks\".");
569
570 // Update parameter in our list.
571 params_->set("Num Recycled Blocks", recycleBlocks_);
572 }
573
574 // Check to see if the timer label changed.
575 if (params->isParameter("Timer Label")) {
576 std::string tempLabel = params->get("Timer Label", label_default_);
577
578 // Update parameter in our list and solver timer
579 if (tempLabel != label_) {
580 label_ = tempLabel;
581 params_->set("Timer Label", label_);
582 std::string solveLabel = label_ + ": RCGSolMgr total solve time";
583#ifdef BELOS_TEUCHOS_TIME_MONITOR
584 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
585#endif
586 }
587 }
588
589 // Check for a change in verbosity level
590 if (params->isParameter("Verbosity")) {
591 if (Teuchos::isParameterType<int>(*params,"Verbosity")) {
592 verbosity_ = params->get("Verbosity", verbosity_default_);
593 } else {
594 verbosity_ = (int)Teuchos::getParameter<Belos::MsgType>(*params,"Verbosity");
595 }
596
597 // Update parameter in our list.
598 params_->set("Verbosity", verbosity_);
599 if (printer_ != Teuchos::null)
600 printer_->setVerbosity(verbosity_);
601 }
602
603 // Check for a change in output style
604 if (params->isParameter("Output Style")) {
605 if (Teuchos::isParameterType<int>(*params,"Output Style")) {
606 outputStyle_ = params->get("Output Style", outputStyle_default_);
607 } else {
608 outputStyle_ = (int)Teuchos::getParameter<Belos::OutputType>(*params,"Output Style");
609 }
610
611 // Reconstruct the convergence test if the explicit residual test is not being used.
612 params_->set("Output Style", outputStyle_);
613 outputTest_ = Teuchos::null;
614 }
615
616 // output stream
617 if (params->isParameter("Output Stream")) {
618 outputStream_ = Teuchos::getParameter<Teuchos::RCP<std::ostream> >(*params,"Output Stream");
619
620 // Update parameter in our list.
621 params_->set("Output Stream", outputStream_);
622 if (printer_ != Teuchos::null)
623 printer_->setOStream( outputStream_ );
624 }
625
626 // frequency level
627 if (verbosity_ & Belos::StatusTestDetails) {
628 if (params->isParameter("Output Frequency")) {
629 outputFreq_ = params->get("Output Frequency", outputFreq_default_);
630 }
631
632 // Update parameter in out list and output status test.
633 params_->set("Output Frequency", outputFreq_);
634 if (outputTest_ != Teuchos::null)
635 outputTest_->setOutputFrequency( outputFreq_ );
636 }
637
638 // Create output manager if we need to.
639 if (printer_ == Teuchos::null) {
640 printer_ = Teuchos::rcp( new OutputManager<ScalarType>(verbosity_, outputStream_) );
641 }
642
643 // Convergence
644 typedef Belos::StatusTestCombo<ScalarType,MV,OP,DM> StatusTestCombo_t;
645 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP,DM> StatusTestResNorm_t;
646
647 // Check for convergence tolerance
648 if (params->isParameter("Convergence Tolerance")) {
649 if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
650 convtol_ = params->get ("Convergence Tolerance",
651 static_cast<MagnitudeType> (DefaultSolverParameters::convTol));
652 }
653 else {
654 convtol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
655 }
656
657 // Update parameter in our list and residual tests.
658 params_->set("Convergence Tolerance", convtol_);
659 if (convTest_ != Teuchos::null)
660 convTest_->setTolerance( convtol_ );
661 }
662
663 if (params->isParameter("Show Maximum Residual Norm Only")) {
664 showMaxResNormOnly_ = Teuchos::getParameter<bool>(*params,"Show Maximum Residual Norm Only");
665
666 // Update parameter in our list and residual tests
667 params_->set("Show Maximum Residual Norm Only", showMaxResNormOnly_);
668 if (convTest_ != Teuchos::null)
669 convTest_->setShowMaxResNormOnly( showMaxResNormOnly_ );
670 }
671
672 // Create status tests if we need to.
673
674 // Basic test checks maximum iterations and native residual.
675 if (maxIterTest_ == Teuchos::null)
676 maxIterTest_ = Teuchos::rcp( new StatusTestMaxIters<ScalarType,MV,OP,DM>( maxIters_ ) );
677
678 // Implicit residual test, using the native residual to determine if convergence was achieved.
679 if (convTest_ == Teuchos::null)
680 convTest_ = Teuchos::rcp( new StatusTestResNorm_t( convtol_, 1 ) );
681
682 if (sTest_ == Teuchos::null) {
683 sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, maxIterTest_, convTest_ ) );
684 if (Teuchos::nonnull(debugStatusTest_)) {
685 // Add the debug convergence test, if it exists.
686 sTest_ = Teuchos::rcp( new StatusTestCombo_t( StatusTestCombo_t::OR, sTest_, debugStatusTest_ ) );
687 }
688 }
689
690 if (outputTest_ == Teuchos::null) {
691
692 // Create the status test output class.
693 // This class manages and formats the output from the status test.
694 StatusTestOutputFactory<ScalarType,MV,OP,DM> stoFactory( outputStyle_ );
695 outputTest_ = stoFactory.create( printer_, sTest_, outputFreq_, Passed+Failed+Undefined );
696
697 // Set the solver string for the output test
698 std::string solverDesc = " Recycling CG ";
699 outputTest_->setSolverDesc( solverDesc );
700 }
701
702 // Create the timer if we need to.
703 if (timerSolve_ == Teuchos::null) {
704 std::string solveLabel = label_ + ": RCGSolMgr total solve time";
705#ifdef BELOS_TEUCHOS_TIME_MONITOR
706 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
707#endif
708 }
709
710 // Inform the solver manager that the current parameters were set.
711 params_Set_ = true;
712}
713
714
715template<class ScalarType, class MV, class OP, class DM>
716Teuchos::RCP<const Teuchos::ParameterList>
717RCGSolMgr<ScalarType,MV,OP,DM,true>::getValidParameters() const
718{
719 static Teuchos::RCP<const Teuchos::ParameterList> validPL;
720
721 // Set all the valid parameters and their default values.
722 if(is_null(validPL)) {
723 Teuchos::RCP<Teuchos::ParameterList> pl = Teuchos::parameterList();
724 pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
725 "The relative residual tolerance that needs to be achieved by the\n"
726 "iterative solver in order for the linear system to be declared converged.");
727 pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
728 "The maximum number of block iterations allowed for each\n"
729 "set of RHS solved.");
730 pl->set("Block Size", static_cast<int>(blockSize_default_),
731 "Block Size Parameter -- currently must be 1 for RCG");
732 pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
733 "The length of a cycle (and this max number of search vectors kept)\n");
734 pl->set("Num Recycled Blocks", static_cast<int>(recycleBlocks_default_),
735 "The number of vectors in the recycle subspace.");
736 pl->set("Verbosity", static_cast<int>(verbosity_default_),
737 "What type(s) of solver information should be outputted\n"
738 "to the output stream.");
739 pl->set("Output Style", static_cast<int>(outputStyle_default_),
740 "What style is used for the solver information outputted\n"
741 "to the output stream.");
742 pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
743 "How often convergence information should be outputted\n"
744 "to the output stream.");
745 pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
746 "A reference-counted pointer to the output stream where all\n"
747 "solver output is sent.");
748 pl->set("Show Maximum Residual Norm Only", static_cast<bool>(showMaxResNormOnly_default_),
749 "When convergence information is printed, only show the maximum\n"
750 "relative residual norm when the block size is greater than one.");
751 pl->set("Timer Label", static_cast<const char *>(label_default_),
752 "The string to use as a prefix for the timer labels.");
753 validPL = pl;
754 }
755 return validPL;
756}
757
758// initializeStateStorage
759template<class ScalarType, class MV, class OP, class DM>
760void RCGSolMgr<ScalarType,MV,OP,DM,true>::initializeStateStorage() {
761
762 // Check if there is any multivector to clone from.
763 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
764 if (rhsMV == Teuchos::null) {
765 // Nothing to do
766 return;
767 }
768 else {
769
770 // Initialize the state storage
771 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
772 "Belos::RCGSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
773
774 // If the subspace has not been initialized before, generate it using the RHS from lp_.
775 if (P_ == Teuchos::null) {
776 P_ = MVT::Clone( *rhsMV, numBlocks_+2 );
777 }
778 else {
779 // Generate P_ by cloning itself ONLY if more space is needed.
780 if (MVT::GetNumberVecs(*P_) < numBlocks_+2) {
781 Teuchos::RCP<const MV> tmp = P_;
782 P_ = MVT::Clone( *tmp, numBlocks_+2 );
783 }
784 }
785
786 // Generate Ap_ only if it doesn't exist
787 if (Ap_ == Teuchos::null)
788 Ap_ = MVT::Clone( *rhsMV, 1 );
789
790 // Generate r_ only if it doesn't exist
791 if (r_ == Teuchos::null)
792 r_ = MVT::Clone( *rhsMV, 1 );
793
794 // Generate z_ only if it doesn't exist
795 if (z_ == Teuchos::null)
796 z_ = MVT::Clone( *rhsMV, 1 );
797
798 // If the recycle space has not been initialized before, generate it using the RHS from lp_.
799 if (U_ == Teuchos::null) {
800 U_ = MVT::Clone( *rhsMV, recycleBlocks_ );
801 }
802 else {
803 // Generate U_ by cloning itself ONLY if more space is needed.
804 if (MVT::GetNumberVecs(*U_) < recycleBlocks_) {
805 Teuchos::RCP<const MV> tmp = U_;
806 U_ = MVT::Clone( *tmp, recycleBlocks_ );
807 }
808 }
809
810 // If the recycle space has not be initialized before, generate it using the RHS from lp_.
811 if (AU_ == Teuchos::null) {
812 AU_ = MVT::Clone( *rhsMV, recycleBlocks_ );
813 }
814 else {
815 // Generate AU_ by cloning itself ONLY if more space is needed.
816 if (MVT::GetNumberVecs(*AU_) < recycleBlocks_) {
817 Teuchos::RCP<const MV> tmp = AU_;
818 AU_ = MVT::Clone( *tmp, recycleBlocks_ );
819 }
820 }
821
822 // If the recycle space has not been initialized before, generate it using the RHS from lp_.
823 if (U1_ == Teuchos::null) {
824 U1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
825 }
826 else {
827 // Generate U1_ by cloning itself ONLY if more space is needed.
828 if (MVT::GetNumberVecs(*U1_) < recycleBlocks_) {
829 Teuchos::RCP<const MV> tmp = U1_;
830 U1_ = MVT::Clone( *tmp, recycleBlocks_ );
831 }
832 }
833
834 // Generate Alpha_ only if it doesn't exist, otherwise resize it.
835 if (Alpha_ == Teuchos::null)
836 Alpha_ = Teuchos::rcp( new std::vector<ScalarType>( numBlocks_, 1 ) );
837 else {
838 if ( (int)Alpha_->size() != numBlocks_ )
839 Alpha_->resize( numBlocks_, 1 );
840 }
841
842 // Generate Beta_ only if it doesn't exist, otherwise resize it.
843 if (Beta_ == Teuchos::null)
844 Beta_ = Teuchos::rcp( new std::vector<ScalarType>( numBlocks_ + 1 ) );
845 else {
846 if ( ((int)Beta_->size() != (numBlocks_+1)) )
847 Beta_->resize( numBlocks_ + 1 );
848 }
849
850 // Generate D_ only if it doesn't exist, otherwise resize it.
851 if (D_ == Teuchos::null)
852 D_ = Teuchos::rcp( new std::vector<ScalarType>( numBlocks_ ) );
853 else {
854 if ( (int)D_->size() != numBlocks_ )
855 D_->resize( numBlocks_ );
856 }
857
858 // Generate Delta_ only if it doesn't exist, otherwise resize it.
859 if (Delta_ == Teuchos::null)
860 Delta_ = DMT::Create( recycleBlocks_, numBlocks_ + 1 );
861 else {
862 if ( (DMT::GetNumRows(*Delta_) != recycleBlocks_) || (DMT::GetNumCols(*Delta_)!= (numBlocks_ + 1)) )
863 DMT::Reshape( *Delta_, recycleBlocks_, numBlocks_ + 1 );
864 }
865
866 // Generate UTAU_ only if it doesn't exist, otherwise resize it.
867 if (UTAU_ == Teuchos::null)
868 UTAU_ = DMT::Create( recycleBlocks_, recycleBlocks_ );
869 else {
870 if ( (DMT::GetNumRows(*UTAU_) != recycleBlocks_) || (DMT::GetNumCols(*UTAU_) != recycleBlocks_) )
871 DMT::Reshape( *UTAU_, recycleBlocks_, recycleBlocks_ );
872 }
873
874 // Generate LUUTAU_ only if it doesn't exist, otherwise resize it.
875 if (LUUTAU_ == Teuchos::null)
876 LUUTAU_ = DMT::Create( recycleBlocks_, recycleBlocks_ );
877 else {
878 if ( (DMT::GetNumRows(*LUUTAU_) != recycleBlocks_) || (DMT::GetNumCols(*LUUTAU_) != recycleBlocks_) )
879 DMT::Reshape( *LUUTAU_, recycleBlocks_, recycleBlocks_ );
880 }
881
882 // Generate ipiv_ only if it doesn't exist, otherwise resize it.
883 if (ipiv_ == Teuchos::null)
884 ipiv_ = Teuchos::rcp( new std::vector<int>(recycleBlocks_) );
885 else {
886 if ( (int)ipiv_->size() != recycleBlocks_ ) // if ipiv not correct size, always resize it
887 ipiv_->resize(recycleBlocks_);
888 }
889
890 // Generate AUTAU_ only if it doesn't exist, otherwise resize it.
891 if (AUTAU_ == Teuchos::null)
892 AUTAU_ = DMT::Create( recycleBlocks_, recycleBlocks_ );
893 else {
894 if ( (DMT::GetNumRows(*AUTAU_) != recycleBlocks_) || (DMT::GetNumCols(*AUTAU_) != recycleBlocks_) )
895 DMT::Reshape( *AUTAU_, recycleBlocks_, recycleBlocks_ );
896 }
897
898 // Generate rTz_old_ only if it doesn't exist
899 if (rTz_old_ == Teuchos::null)
900 rTz_old_ = Teuchos::rcp( new std::vector<ScalarType>(1) );
901 else {
902 if ( (rTz_old_->size() != 1) )
903 rTz_old_->resize( 1 );
904 }
905
906 // Generate F_ only if it doesn't exist
907 if (F_ == Teuchos::null)
908 F_ = DMT::Create( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
909 else {
910 if ( (DMT::GetNumRows(*F_) != (numBlocks_+recycleBlocks_)) || (DMT::GetNumCols(*F_) != numBlocks_+recycleBlocks_) )
911 DMT::Reshape( *F_, numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
912 }
913
914 // Generate G_ only if it doesn't exist
915 if (G_ == Teuchos::null)
916 G_ = DMT::Create( numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
917 else {
918 if ( (DMT::GetNumRows(*G_) != (numBlocks_+recycleBlocks_)) || (DMT::GetNumCols(*G_) != numBlocks_+recycleBlocks_) )
919 DMT::Reshape( *G_, numBlocks_+recycleBlocks_, numBlocks_+recycleBlocks_ );
920 }
921
922 // Generate Y_ only if it doesn't exist
923 if (Y_ == Teuchos::null)
924 Y_ = DMT::Create( numBlocks_+recycleBlocks_, recycleBlocks_ );
925 else {
926 if ( (DMT::GetNumRows(*Y_) != (numBlocks_+recycleBlocks_)) || (DMT::GetNumCols(*Y_) != recycleBlocks_) )
927 DMT::Reshape( *Y_, numBlocks_+recycleBlocks_, recycleBlocks_ );
928 }
929
930 // Generate L2_ only if it doesn't exist
931 if (L2_ == Teuchos::null)
932 L2_ = DMT::Create( numBlocks_+1, numBlocks_ );
933 else {
934 if ( (DMT::GetNumRows(*L2_) != (numBlocks_+1)) || (DMT::GetNumCols(*L2_) != numBlocks_) )
935 DMT::Reshape( *L2_, numBlocks_+1, numBlocks_ );
936 }
937
938 // Generate DeltaL2_ only if it doesn't exist
939 if (DeltaL2_ == Teuchos::null)
940 DeltaL2_ = DMT::Create( recycleBlocks_, numBlocks_ );
941 else {
942 if ( (DMT::GetNumRows(*DeltaL2_) != recycleBlocks_) || (DMT::GetNumCols(*DeltaL2_) != numBlocks_ ) )
943 DMT::Reshape( *DeltaL2_, recycleBlocks_, numBlocks_ );
944 }
945
946 // Generate AU1TUDeltaL2_ only if it doesn't exist
947 if (AU1TUDeltaL2_ == Teuchos::null)
948 AU1TUDeltaL2_ = DMT::Create( recycleBlocks_, numBlocks_ );
949 else {
950 if ( (DMT::GetNumRows(*AU1TUDeltaL2_) != recycleBlocks_) || (DMT::GetNumCols(*AU1TUDeltaL2_) != numBlocks_ ) )
951 DMT::Reshape( *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
952 }
953
954 // Generate AU1TAU1_ only if it doesn't exist
955 if (AU1TAU1_ == Teuchos::null)
956 AU1TAU1_ = DMT::Create( recycleBlocks_, recycleBlocks_ );
957 else {
958 if ( (DMT::GetNumRows(*AU1TAU1_) != recycleBlocks_) || (DMT::GetNumCols(*AU1TAU1_) != recycleBlocks_ ) )
959 DMT::Reshape( *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
960 }
961
962 // Generate GY_ only if it doesn't exist
963 if (GY_ == Teuchos::null)
964 GY_ = DMT::Create( numBlocks_ + recycleBlocks_, recycleBlocks_ );
965 else {
966 if ( (DMT::GetNumRows(*GY_) != (numBlocks_ + recycleBlocks_)) || (DMT::GetNumCols(*GY_) != recycleBlocks_ ) )
967 DMT::Reshape( *GY_, numBlocks_+recycleBlocks_, recycleBlocks_ );
968 }
969
970 // Generate AU1TU1_ only if it doesn't exist
971 if (AU1TU1_ == Teuchos::null)
972 AU1TU1_ = DMT::Create( recycleBlocks_, recycleBlocks_ );
973 else {
974 if ( (DMT::GetNumRows(*AU1TU1_) != recycleBlocks_) || (DMT::GetNumCols(*AU1TU1_) != recycleBlocks_ ) )
975 DMT::Reshape( *AU1TU1_, recycleBlocks_, recycleBlocks_ );
976 }
977
978 // Generate FY_ only if it doesn't exist
979 if (FY_ == Teuchos::null)
980 FY_ = DMT::Create( numBlocks_ + recycleBlocks_, recycleBlocks_ );
981 else {
982 if ( (DMT::GetNumRows(*FY_) != (numBlocks_ + recycleBlocks_)) || (DMT::GetNumCols(*FY_) != recycleBlocks_ ) )
983 DMT::Reshape( *FY_, numBlocks_+recycleBlocks_, recycleBlocks_ );
984 }
985
986 // Generate AU1TAP_ only if it doesn't exist
987 if (AU1TAP_ == Teuchos::null)
988 AU1TAP_ = DMT::Create( recycleBlocks_, numBlocks_ );
989 else {
990 if ( (DMT::GetNumRows(*AU1TAP_) != recycleBlocks_) || (DMT::GetNumCols(*AU1TAP_) != numBlocks_ ) )
991 DMT::Reshape( *AU1TAP_, recycleBlocks_, numBlocks_ );
992 }
993
994 // Generate APTAP_ only if it doesn't exist
995 if (APTAP_ == Teuchos::null)
996 APTAP_ = DMT::Create( numBlocks_, numBlocks_ );
997 else {
998 if ( (DMT::GetNumRows(*APTAP_) != numBlocks_) || (DMT::GetNumCols(*APTAP_) != (numBlocks_) ) )
999 DMT::Reshape( *APTAP_, numBlocks_, numBlocks_ );
1000 }
1001
1002 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1003 if (U1Y1_ == Teuchos::null) {
1004 U1Y1_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1005 }
1006 else {
1007 // Generate U1Y1_ by cloning itself ONLY if more space is needed.
1008 if (MVT::GetNumberVecs(*U1Y1_) < recycleBlocks_) {
1009 Teuchos::RCP<const MV> tmp = U1Y1_;
1010 U1Y1_ = MVT::Clone( *tmp, recycleBlocks_ );
1011 }
1012 }
1013
1014 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1015 if (PY2_ == Teuchos::null) {
1016 PY2_ = MVT::Clone( *rhsMV, recycleBlocks_ );
1017 }
1018 else {
1019 // Generate PY2_ by cloning itself ONLY if more space is needed.
1020 if (MVT::GetNumberVecs(*PY2_) < recycleBlocks_) {
1021 Teuchos::RCP<const MV> tmp = PY2_;
1022 PY2_ = MVT::Clone( *tmp, recycleBlocks_ );
1023 }
1024 }
1025
1026 // Generate AUTAP_ only if it doesn't exist
1027 if (AUTAP_ == Teuchos::null)
1028 AUTAP_ = DMT::Create( recycleBlocks_, numBlocks_ );
1029 else {
1030 if ( (DMT::GetNumRows(*AUTAP_) != recycleBlocks_) || (DMT::GetNumCols(*AUTAP_) != numBlocks_ ) )
1031 DMT::Reshape( *AUTAP_, recycleBlocks_, numBlocks_ );
1032 }
1033
1034 // Generate AU1TU_ only if it doesn't exist
1035 if (AU1TU_ == Teuchos::null)
1036 AU1TU_ = DMT::Create( recycleBlocks_, recycleBlocks_ );
1037 else {
1038 if ( (DMT::GetNumRows(*AU1TU_) != recycleBlocks_) || (DMT::GetNumCols(*AU1TU_) != recycleBlocks_ ) )
1039 DMT::Reshape( *AU1TU_, recycleBlocks_, recycleBlocks_ );
1040 }
1041
1042
1043 }
1044}
1045
1046template<class ScalarType, class MV, class OP, class DM>
1047ReturnType RCGSolMgr<ScalarType,MV,OP,DM,true>::solve() {
1048 ReturnType retType = Undetermined;
1049
1050 Teuchos::BLAS<int,ScalarType> blas;
1051 Teuchos::LAPACK<int,ScalarType> lapack;
1052 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1053 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1054
1055 std::vector<int> index(1), rindex(recycleBlocks_), nindex(numBlocks_);
1056 for (int i=0; i<recycleBlocks_; ++i) { rindex[i] = i; }
1057 for (int i=0; i<numBlocks_; ++i) { nindex[i] = i; }
1058
1059 // Count of number of cycles performed on current rhs
1060 int cycle = 0;
1061
1062 // Set the current parameters if they were not set before.
1063 // NOTE: This may occur if the user generated the solver manager with the default constructor and
1064 // then didn't set any parameters using setParameters().
1065 if (!params_Set_) {
1066 setParameters(Teuchos::parameterList(*getValidParameters()));
1067 }
1068
1069 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,RCGSolMgrLinearProblemFailure,
1070 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1071 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),RCGSolMgrLinearProblemFailure,
1072 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1073 TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1074 RCGSolMgrLinearProblemFailure,
1075 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1076
1077 // Grab the preconditioning object
1078 Teuchos::RCP<OP> precObj;
1079 if (problem_->getLeftPrec() != Teuchos::null) {
1080 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1081 }
1082 else if (problem_->getRightPrec() != Teuchos::null) {
1083 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1084 }
1085
1086 // Create indices for the linear systems to be solved.
1087 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1088 std::vector<int> currIdx(1);
1089 currIdx[0] = 0;
1090
1091 // Inform the linear problem of the current linear system to solve.
1092 problem_->setLSIndex( currIdx );
1093
1094 // Check the number of blocks and change them if necessary.
1095 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1096 if (numBlocks_ > dim) {
1097 numBlocks_ = Teuchos::asSafe<int>(dim);
1098 params_->set("Num Blocks", numBlocks_);
1099 printer_->stream(Warnings) <<
1100 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1101 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1102 }
1103
1104 // Initialize storage for all state variables
1105 initializeStateStorage();
1106
1107 // Parameter list
1108 Teuchos::ParameterList plist;
1109 plist.set("Num Blocks",numBlocks_);
1110 plist.set("Recycled Blocks",recycleBlocks_);
1111
1112 // Reset the status test.
1113 outputTest_->reset();
1114
1115 // Assume convergence is achieved, then let any failed convergence set this to false.
1116 bool isConverged = true;
1117
1118 // Compute AU = A*U, UTAU = U'*AU, AUTAU = (AU)'*(AU)
1119 if (existU_) {
1120 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, rindex );
1121 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, rindex );
1122 // Initialize AU
1123 problem_->applyOp( *Utmp, *AUtmp );
1124 // Initialize UTAU
1125 MVT::MvTransMv( one, *Utmp, *AUtmp, *UTAU_ );
1126 // Initialize AUTAU ( AUTAU = AU'*(M\AU) )
1127 if ( precObj != Teuchos::null ) {
1128 Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, rindex ); // use U1 as temp storage
1129 OPT::Apply( *precObj, *AUtmp, *PCAU );
1130 MVT::MvTransMv( one, *AUtmp, *PCAU, *AUTAU_ );
1131 } else {
1132 MVT::MvTransMv( one, *AUtmp, *AUtmp, *AUTAU_ );
1133 }
1134 }
1135
1137 // RCG solver
1138
1139 Teuchos::RCP<RCGIter<ScalarType,MV,OP,DM> > rcg_iter;
1140 rcg_iter = Teuchos::rcp( new RCGIter<ScalarType,MV,OP,DM>(problem_,printer_,outputTest_,plist) );
1141
1142 // Enter solve() iterations
1143 {
1144#ifdef BELOS_TEUCHOS_TIME_MONITOR
1145 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1146#endif
1147
1148 while ( numRHS2Solve > 0 ) {
1149
1150 // Debugging output to tell use if recycle space exists and will be used
1151 if (printer_->isVerbosity( Debug ) ) {
1152 if (existU_) printer_->print( Debug, "Using recycle space generated from previous call to solve()." );
1153 else printer_->print( Debug, "No recycle space exists." );
1154 }
1155
1156 // Reset the number of iterations.
1157 rcg_iter->resetNumIters();
1158
1159 // Set the current number of recycle blocks and subspace dimension with the RCG iteration.
1160 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1161
1162 // Reset the number of calls that the status test output knows about.
1163 outputTest_->resetNumCalls();
1164
1165 // indicate that updated recycle space has not yet been generated for this linear system
1166 existU1_ = false;
1167
1168 // reset cycle count
1169 cycle = 0;
1170
1171 // Get the current residual
1172 problem_->computeCurrResVec( &*r_ );
1173
1174 // If U exists, find best soln over this space first
1175 if (existU_) {
1176 // Solve linear system UTAU * y = (U'*r)
1177 Teuchos::RCP<DM> Utr = DMT::Create(recycleBlocks_,1);
1178 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, rindex );
1179 MVT::MvTransMv( one, *Utmp, *r_, *Utr );
1180
1181 DMT::SyncHostToDevice(*LUUTAU_);
1182 DMT::Assign(*LUUTAU_,*UTAU_);
1183 DMT::SyncDeviceToHost( *LUUTAU_ );
1184 DMT::SyncDeviceToHost( *Utr );
1185 int info = 0;
1186 lapack.GESV(recycleBlocks_, 1, DMT::GetRawHostPtr(*LUUTAU_), DMT::GetStride(*LUUTAU_),
1187 &(*ipiv_)[0], DMT::GetRawHostPtr(*Utr), DMT::GetStride(*Utr), &info);
1188 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1189 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1190 DMT::SyncHostToDevice( *Utr );
1191
1192 // Update solution (x = x + U*y)
1193 MVT::MvTimesMatAddMv( one, *Utmp, *Utr, one, *problem_->getCurrLHSVec() );
1194
1195 // Update residual ( r = r - AU*y )
1196 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, rindex );
1197 MVT::MvTimesMatAddMv( -one, *AUtmp, *Utr, one, *r_ );
1198 }
1199
1200 if ( precObj != Teuchos::null ) {
1201 OPT::Apply( *precObj, *r_, *z_ );
1202 } else {
1203 z_ = r_;
1204 }
1205
1206 // rTz_old = r'*z
1207 MVT::MvDot( *r_, *z_, *rTz_old_ );
1208
1209 if ( existU_ ) {
1210 // mu = UTAU\‍(AU'*z);
1211 Teuchos::RCP<DM> mu = DMT::Subview(*Delta_, recycleBlocks_, 1);
1212 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, rindex );
1213 MVT::MvTransMv( one, *AUtmp, *z_, *mu );
1214
1215 DMT::SyncDeviceToHost( *Delta_ );
1216 char TRANS = 'N';
1217 int info;
1218 lapack.GETRS( TRANS, recycleBlocks_, 1, DMT::GetConstRawHostPtr(*LUUTAU_), DMT::GetStride(*LUUTAU_),
1219 &(*ipiv_)[0], DMT::GetRawHostPtr(*mu), DMT::GetStride(*mu), &info );
1220 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1221 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1222 DMT::SyncHostToDevice( *mu );
1223
1224 // p = z - U*mu;
1225 index.resize( 1 );
1226 index[0] = 0;
1227 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1228 MVT::Assign(*z_,*Ptmp);
1229 MVT::MvTimesMatAddMv( -one, *U_, *mu, one, *Ptmp );
1230 } else {
1231 // p = z;
1232 index.resize( 1 );
1233 index[0] = 0;
1234 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1235 MVT::Assign(*z_,*Ptmp);
1236 }
1237
1238 // Set the new state and initialize the solver.
1239 RCGIterState<ScalarType,MV,DM> newstate;
1240
1241 // Create RCP views here
1242 index.resize( numBlocks_+1 );
1243 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1244 newstate.P = MVT::CloneViewNonConst( *P_, index );
1245 newstate.U = MVT::CloneViewNonConst( *U_, rindex );
1246 newstate.AU = MVT::CloneViewNonConst( *AU_, rindex );
1247 newstate.Alpha = Alpha_;
1248 newstate.Beta = Beta_;
1249 newstate.Beta_i = 0;
1250 newstate.D = D_;
1251 newstate.Delta = Delta_;
1252 newstate.LUUTAU = LUUTAU_;
1253 // assign the rest of the values in the struct
1254 newstate.curDim = 1; // We have initialized the first search vector
1255 newstate.Ap = Ap_;
1256 newstate.r = r_;
1257 newstate.z = z_;
1258 newstate.existU = existU_;
1259 newstate.ipiv = ipiv_;
1260 newstate.rTz_old = rTz_old_;
1261
1262 rcg_iter->initialize(newstate);
1263
1264 while(1) {
1265
1266 // tell rcg_iter to iterate
1267 try {
1268 rcg_iter->iterate();
1269
1271 //
1272 // check convergence first
1273 //
1275 if ( convTest_->getStatus() == Passed ) {
1276 // We have convergence
1277 break; // break from while(1){rcg_iter->iterate()}
1278 }
1280 //
1281 // check for maximum iterations
1282 //
1284 else if ( maxIterTest_->getStatus() == Passed ) {
1285 // we don't have convergence
1286 retType = MaxItersReached;
1287 isConverged = false;
1288 break; // break from while(1){rcg_iter->iterate()}
1289 }
1291 //
1292 // check if cycle complete; update for next cycle
1293 //
1295 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1296 // index into P_ of last search vector generated this cycle
1297 int lastp = -1;
1298 // index into Beta_ of last entry generated this cycle
1299 int lastBeta = -1;
1300 if (recycleBlocks_ > 0) {
1301 if (!existU_) {
1302 if (cycle == 0) { // No U, no U1
1303
1304 Teuchos::RCP<DM> Ftmp = DMT::Subview( *F_, numBlocks_, numBlocks_ );
1305 Teuchos::RCP<DM> Gtmp = DMT::Subview( *G_, numBlocks_, numBlocks_ );
1306 DMT::PutScalar( *Ftmp, zero );
1307 DMT::PutScalar( *Gtmp, zero );
1308 DMT::SyncDeviceToHost( *F_ );
1309 DMT::SyncDeviceToHost( *G_ );
1310 for (int ii=0;ii<numBlocks_;ii++) {
1311 DMT::Value(*Gtmp,ii,ii) = ((*D_)[ii] / (*Alpha_)[ii])*(1 + (*Beta_)[ii]);
1312 if (ii > 0) {
1313 DMT::Value(*Gtmp,ii-1,ii) = -(*D_)[ii]/(*Alpha_)[ii-1];
1314 DMT::Value(*Gtmp,ii,ii-1) = -(*D_)[ii]/(*Alpha_)[ii-1];
1315 }
1316 DMT::Value(*Ftmp,ii,ii) = (*D_)[ii];
1317 }
1318 DMT::SyncHostToDevice( *F_ );
1319 DMT::SyncHostToDevice( *G_ );
1320
1321 // compute harmonic Ritz vectors
1322 DMT::SyncDeviceToHost( *Y_ );
1323 Teuchos::RCP<DM> Ytmp = DMT::Subview( *Y_, numBlocks_, recycleBlocks_ );
1324 getHarmonicVecs(*Ftmp,*Gtmp,*Ytmp);
1325 DMT::SyncHostToDevice( *Y_ );
1326
1327 // U1 = [P(:,1:end-1)*Y];
1328 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, nindex );
1329 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, rindex );
1330 MVT::MvTimesMatAddMv( one, *Ptmp, *Ytmp, zero, *U1tmp );
1331
1332 // Precompute some variables for next cycle
1333 DMT::SyncDeviceToHost(*GY_);
1334 DMT::SyncDeviceToHost(*AU1TAU1_);
1335 DMT::SyncDeviceToHost(*FY_);
1336 DMT::SyncDeviceToHost(*AU1TU1_);
1337
1338 // AU1TAU1 = Y'*G*Y;
1339 Teuchos::RCP<DM> GYtmp = DMT::Subview( *GY_, numBlocks_, recycleBlocks_ );
1340 //GYtmp->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,*Ytmp,zero);
1341 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, numBlocks_, recycleBlocks_, numBlocks_,
1342 one, DMT::GetConstRawHostPtr(*Gtmp), DMT::GetStride(*Gtmp),
1343 DMT::GetConstRawHostPtr(*Ytmp), DMT::GetStride(*Ytmp),
1344 zero, DMT::GetRawHostPtr(*GYtmp), DMT::GetStride(*GYtmp));
1345 //AU1TAU1_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Ytmp,*GYtmp,zero);
1346 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_, numBlocks_,
1347 one, DMT::GetConstRawHostPtr(*Ytmp), DMT::GetStride(*Ytmp),
1348 DMT::GetConstRawHostPtr(*GYtmp), DMT::GetStride(*GYtmp),
1349 zero, DMT::GetRawHostPtr(*AU1TAU1_), DMT::GetStride(*AU1TAU1_));
1350
1351
1352 // AU1TU1 = Y'*F*Y;
1353 Teuchos::RCP<DM> FYtmp = DMT::Subview( *FY_, numBlocks_, recycleBlocks_ );
1354 //FYtmp->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Ftmp,*Ytmp,zero);
1355 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, numBlocks_, recycleBlocks_, numBlocks_,
1356 one, DMT::GetConstRawHostPtr(*Ftmp), DMT::GetStride(*Ftmp),
1357 DMT::GetConstRawHostPtr(*Ytmp), DMT::GetStride(*Ytmp),
1358 zero, DMT::GetRawHostPtr(*FYtmp), DMT::GetStride(*FYtmp));
1359 //AU1TU1_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Ytmp,*FYtmp,zero);
1360 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_, numBlocks_,
1361 one, DMT::GetConstRawHostPtr(*Ytmp), DMT::GetStride(*Ytmp),
1362 DMT::GetConstRawHostPtr(*FYtmp), DMT::GetStride(*FYtmp),
1363 zero, DMT::GetRawHostPtr(*AU1TU1_), DMT::GetStride(*AU1TU1_));
1364
1365 DMT::SyncHostToDevice(*AU1TAU1_);
1366 DMT::SyncHostToDevice(*AU1TU1_);
1367 DMT::SyncHostToDevice(*AU1TAP_);
1368
1369 Teuchos::RCP<DM> AU1TAPtmp = DMT::Subview( *AU1TAP_, recycleBlocks_, 1 );
1370 // Must reinitialize AU1TAP; can become dense later
1371 DMT::PutScalar( *AU1TAPtmp, zero );
1372 // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
1373 DMT::SyncDeviceToHost( *AU1TAP_ );
1374 ScalarType alphatmp = -1.0 / (*Alpha_)[numBlocks_-1];
1375 for (int ii=0; ii<recycleBlocks_; ++ii) {
1376 DMT::Value(*AU1TAPtmp,ii,0) = DMT::ValueConst(*Ytmp,numBlocks_-1,ii) * alphatmp;
1377 }
1378 DMT::SyncHostToDevice(*AU1TAP_);
1379
1380 // indicate that updated recycle space now defined
1381 existU1_ = true;
1382
1383 // Indicate the size of the P, Beta structures generated this cycle
1384 lastp = numBlocks_;
1385 lastBeta = numBlocks_-1;
1386
1387 } // if (cycle == 0)
1388 else { // No U, but U1 guaranteed to exist now
1389
1390 // Finish computation of subblocks
1391 // AU1TAP = AU1TAP * D(1);
1392 DMT::Scale(*AU1TAP_,(*D_)[0]);
1393
1394 DMT::PutScalar(*APTAP_,zero);
1395 DMT::SyncDeviceToHost(*APTAP_);
1396 for (int ii=0; ii<numBlocks_; ii++) {
1397 DMT::Value(*APTAP_,ii,ii) = ((*D_)[ii] / (*Alpha_)[ii])*(1 + (*Beta_)[ii+1]);
1398 if (ii > 0) {
1399 DMT::Value(*APTAP_,ii-1,ii) = -(*D_)[ii]/(*Alpha_)[ii-1];
1400 DMT::Value(*APTAP_,ii,ii-1) = -(*D_)[ii]/(*Alpha_)[ii-1];
1401 }
1402 }
1403 DMT::SyncHostToDevice(*APTAP_);
1404
1405 // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
1406 DMT::PutScalar(*F_,zero);
1407 Teuchos::RCP<DM> F11 = DMT::Subview( *F_, recycleBlocks_, recycleBlocks_ );
1408 Teuchos::RCP<DM> F22 = DMT::Subview( *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1409 DMT::Assign(*F11,*AU1TU1_);
1410 DMT::SyncDeviceToHost(*F_);
1411 for(int ii=0;ii<numBlocks_;ii++) {
1412 DMT::Value(*F22,ii,ii) = (*D_)[ii];
1413 }
1414 DMT::SyncHostToDevice(*F_);
1415
1416 // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
1417 Teuchos::RCP<DM> G11 = DMT::Subview( *G_, recycleBlocks_, recycleBlocks_ );
1418 Teuchos::RCP<DM> G12 = DMT::Subview( *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1419 Teuchos::RCP<DM> G21 = DMT::Subview( *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1420 Teuchos::RCP<DM> G22 = DMT::Subview( *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1421 DMT::Assign(*G11, *AU1TAU1_);
1422 DMT::Assign(*G12, *AU1TAP_);
1423 DMT::Assign(*G22, *APTAP_);
1424 DMT::SyncDeviceToHost( *G_ );
1425 // G21 = G12'; (no transpose operator exists for DM; Do copy manually)
1426 for (int ii=0;ii<recycleBlocks_;++ii)
1427 for (int jj=0;jj<numBlocks_;++jj)
1428 DMT::Value(*G21,jj,ii) = DMT::ValueConst(*G12,ii,jj);
1429 DMT::SyncHostToDevice( *G_ );
1430
1431 // compute harmonic Ritz vectors
1432 getHarmonicVecs(*F_,*G_,*Y_);
1433 DMT::SyncHostToDevice( *Y_ );
1434
1435 // U1 = [U1 P(:,2:end-1)]*Y;
1436 index.resize( numBlocks_ );
1437 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1438 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1439 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, rindex );
1440 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, rindex );
1441 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, rindex );
1442 Teuchos::RCP<const DM> Y1 = DMT::SubviewConst( *Y_, recycleBlocks_, recycleBlocks_ );
1443 Teuchos::RCP<const DM> Y2 = DMT::SubviewConst( *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1444 MVT::MvTimesMatAddMv( one, *Ptmp, *Y2, zero, *PY2tmp );
1445 MVT::MvTimesMatAddMv( one, *U1tmp, *Y1, zero, *U1Y1tmp );
1446 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1447
1448 // Precompute some variables for next cycle
1449 DMT::SyncDeviceToHost(*GY_);
1450 DMT::SyncDeviceToHost(*AU1TAU1_);
1451
1452 // AU1TAU1 = Y'*G*Y;
1453 //GY_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*G_,*Y_,zero);
1454 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, numBlocks_+recycleBlocks_,
1455 recycleBlocks_, numBlocks_+recycleBlocks_,
1456 one, DMT::GetConstRawHostPtr(*G_), DMT::GetStride(*G_),
1457 DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1458 zero, DMT::GetRawHostPtr(*GY_), DMT::GetStride(*GY_));
1459 //AU1TAU1_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Y_,*GY_,zero);
1460 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_,
1461 numBlocks_+recycleBlocks_,
1462 one, DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1463 DMT::GetConstRawHostPtr(*GY_), DMT::GetStride(*GY_),
1464 zero, DMT::GetRawHostPtr(*AU1TAU1_), DMT::GetStride(*AU1TAU1_));
1465
1466 DMT::SyncHostToDevice(*GY_);
1467 DMT::SyncHostToDevice(*AU1TAU1_);
1468
1469 // AU1TAP = zeros(k,m);
1470 // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
1471 DMT::PutScalar(*AU1TAP_,zero);
1472 ScalarType alphatmp = -1.0 / (*Alpha_)[numBlocks_-1];
1473 DMT::SyncDeviceToHost(*AU1TAP_);
1474 for (int ii=0; ii<recycleBlocks_; ++ii) {
1475 DMT::Value(*AU1TAP_,ii,0) = DMT::ValueConst(*Y_,numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1476 }
1477 DMT::SyncHostToDevice(*AU1TAP_);
1478
1479 DMT::SyncDeviceToHost(*FY_);
1480 DMT::SyncDeviceToHost(*AU1TU1_);
1481
1482 // AU1TU1 = Y'*F*Y;
1483 //FY_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*F_,*Y_,zero);
1484 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, numBlocks_+recycleBlocks_,
1485 recycleBlocks_, numBlocks_+recycleBlocks_,
1486 one, DMT::GetConstRawHostPtr(*F_), DMT::GetStride(*F_),
1487 DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1488 zero, DMT::GetRawHostPtr(*FY_), DMT::GetStride(*FY_));
1489 //AU1TU1_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Y_,*FY_,zero);
1490 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_,
1491 numBlocks_+recycleBlocks_,
1492 one, DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1493 DMT::GetConstRawHostPtr(*FY_), DMT::GetStride(*FY_),
1494 zero, DMT::GetRawHostPtr(*AU1TU1_), DMT::GetStride(*AU1TU1_));
1495
1496 DMT::SyncHostToDevice(*FY_);
1497 DMT::SyncHostToDevice(*AU1TU1_);
1498
1499 // Indicate the size of the P, Beta structures generated this cycle
1500 lastp = numBlocks_+1;
1501 lastBeta = numBlocks_;
1502
1503 } // if (cycle != 1)
1504 } // if (!existU_)
1505 else { // U exists
1506 if (cycle == 0) { // No U1, but U exists
1507 DMT::PutScalar(*APTAP_,zero);
1508 DMT::SyncDeviceToHost(*APTAP_);
1509 for (int ii=0; ii<numBlocks_; ii++) {
1510 DMT::Value(*APTAP_,ii,ii) = ((*D_)[ii] / (*Alpha_)[ii])*(1 + (*Beta_)[ii]);
1511 if (ii > 0) {
1512 DMT::Value(*APTAP_,ii-1,ii) = -(*D_)[ii]/(*Alpha_)[ii-1];
1513 DMT::Value(*APTAP_,ii,ii-1) = -(*D_)[ii]/(*Alpha_)[ii-1];
1514 }
1515 }
1516 DMT::SyncHostToDevice(*APTAP_);
1517
1518 DMT::PutScalar(*L2_,zero);
1519 DMT::SyncDeviceToHost(*L2_);
1520 for(int ii=0;ii<numBlocks_;ii++) {
1521 DMT::Value(*L2_,ii,ii) = 1./(*Alpha_)[ii];
1522 DMT::Value(*L2_,ii+1,ii) = -1./(*Alpha_)[ii];
1523 }
1524 DMT::SyncHostToDevice(*L2_);
1525
1526 // AUTAP = UTAU*Delta*L2;
1527 DMT::SyncDeviceToHost(*Delta_);
1528 DMT::SyncDeviceToHost(*DeltaL2_);
1529 DMT::SyncDeviceToHost(*AUTAP_);
1530 DMT::SyncDeviceToHost(*UTAU_);
1531
1532 //DeltaL2_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Delta_,*L2_,zero);
1533 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, recycleBlocks_, numBlocks_, numBlocks_+1,
1534 one, DMT::GetConstRawHostPtr(*Delta_), DMT::GetStride(*Delta_),
1535 DMT::GetConstRawHostPtr(*L2_), DMT::GetStride(*L2_),
1536 zero, DMT::GetRawHostPtr(*DeltaL2_), DMT::GetStride(*DeltaL2_));
1537 //AUTAP_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*UTAU_,*DeltaL2_,zero);
1538 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, recycleBlocks_, numBlocks_, recycleBlocks_,
1539 one, DMT::GetConstRawHostPtr(*UTAU_), DMT::GetStride(*UTAU_),
1540 DMT::GetConstRawHostPtr(*DeltaL2_), DMT::GetStride(*DeltaL2_),
1541 zero, DMT::GetRawHostPtr(*AUTAP_), DMT::GetStride(*AUTAP_));
1542
1543 DMT::SyncHostToDevice(*DeltaL2_);
1544 DMT::SyncHostToDevice(*AUTAP_);
1545
1546 // F = [UTAU zeros(k,m); zeros(m,k) diag(D)];
1547 DMT::PutScalar(*F_,zero);
1548 Teuchos::RCP<DM> F11 = DMT::Subview( *F_, recycleBlocks_, recycleBlocks_ );
1549 Teuchos::RCP<DM> F22 = DMT::Subview( *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1550 DMT::Assign(*F11,*UTAU_);
1551 DMT::SyncDeviceToHost(*F_);
1552 for(int ii=0;ii<numBlocks_;ii++) {
1553 DMT::Value(*F22,ii,ii) = (*D_)[ii];
1554 }
1555 DMT::SyncHostToDevice(*F_);
1556
1557 // G = [AUTAU AUTAP; AUTAP' APTAP];
1558 Teuchos::RCP<DM> G11 = DMT::Subview( *G_, recycleBlocks_, recycleBlocks_ );
1559 Teuchos::RCP<DM> G12 = DMT::Subview( *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1560 Teuchos::RCP<DM> G21 = DMT::Subview( *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1561 Teuchos::RCP<DM> G22 = DMT::Subview( *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1562 DMT::Assign(*G11,*AUTAU_);
1563 DMT::Assign(*G12,*AUTAP_);
1564 DMT::Assign(*G22,*APTAP_);
1565 DMT::SyncDeviceToHost(*G_);
1566 // G21 = G12'; (no transpose operator exists for DM; Do copy manually)
1567 for (int ii=0;ii<recycleBlocks_;++ii)
1568 for (int jj=0;jj<numBlocks_;++jj)
1569 DMT::Value(*G21,jj,ii) = DMT::ValueConst(*G12,ii,jj);
1570 DMT::SyncHostToDevice(*G_);
1571
1572 // compute harmonic Ritz vectors
1573 getHarmonicVecs(*F_,*G_,*Y_);
1574 DMT::SyncHostToDevice(*Y_);
1575
1576 // U1 = [U P(:,1:end-1)]*Y;
1577 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, rindex );
1578 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, nindex );
1579 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, rindex );
1580 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, rindex );
1581 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, rindex );
1582 Teuchos::RCP<const DM> Y1 = DMT::SubviewConst( *Y_, recycleBlocks_, recycleBlocks_ );
1583 Teuchos::RCP<const DM> Y2 = DMT::SubviewConst( *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1584 MVT::MvTimesMatAddMv( one, *Ptmp, *Y2, zero, *PY2tmp );
1585 MVT::MvTimesMatAddMv( one, *Utmp, *Y1, zero, *UY1tmp );
1586 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1587
1588 // Precompute some variables for next cycle
1589 DMT::SyncDeviceToHost(*GY_);
1590 DMT::SyncDeviceToHost(*AU1TAU1_);
1591 DMT::SyncDeviceToHost(*FY_);
1592 DMT::SyncDeviceToHost(*AU1TU1_);
1593
1594 // AU1TAU1 = Y'*G*Y;
1595 //GY_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*G_,*Y_,zero);
1596 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, numBlocks_+recycleBlocks_,
1597 recycleBlocks_, numBlocks_+recycleBlocks_,
1598 one, DMT::GetConstRawHostPtr(*G_), DMT::GetStride(*G_),
1599 DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1600 zero, DMT::GetRawHostPtr(*GY_), DMT::GetStride(*GY_));
1601 //AU1TAU1_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Y_,*GY_,zero);
1602 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_,
1603 numBlocks_+recycleBlocks_,
1604 one, DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1605 DMT::GetConstRawHostPtr(*GY_), DMT::GetStride(*GY_),
1606 zero, DMT::GetRawHostPtr(*AU1TAU1_), DMT::GetStride(*AU1TAU1_));
1607
1608 // AU1TU1 = Y'*F*Y;
1609 //FY_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*F_,*Y_,zero);
1610 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, numBlocks_+recycleBlocks_,
1611 recycleBlocks_, numBlocks_+recycleBlocks_,
1612 one, DMT::GetConstRawHostPtr(*F_), DMT::GetStride(*F_),
1613 DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1614 zero, DMT::GetRawHostPtr(*FY_), DMT::GetStride(*FY_));
1615 //AU1TU1_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Y_,*FY_,zero);
1616 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_,
1617 numBlocks_+recycleBlocks_,
1618 one, DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1619 DMT::GetConstRawHostPtr(*FY_), DMT::GetStride(*FY_),
1620 zero, DMT::GetRawHostPtr(*AU1TU1_), DMT::GetStride(*AU1TU1_));
1621
1622 DMT::SyncHostToDevice(*GY_);
1623 DMT::SyncHostToDevice(*AU1TAU1_);
1624 DMT::SyncHostToDevice(*FY_);
1625 DMT::SyncHostToDevice(*AU1TU1_);
1626
1627 // AU1TU = UTAU;
1628 DMT::Assign(*AU1TU_,*UTAU_);
1629
1630 // dold = D(end);
1631 dold = (*D_)[numBlocks_-1];
1632
1633 // indicate that updated recycle space now defined
1634 existU1_ = true;
1635
1636 // Indicate the size of the P, Beta structures generated this cycle
1637 lastp = numBlocks_;
1638 lastBeta = numBlocks_-1;
1639 }
1640 else { // Have U and U1
1641
1642 DMT::SyncDeviceToHost(*APTAP_);
1643 for (int ii=0; ii<numBlocks_; ii++) {
1644 DMT::Value(*APTAP_,ii,ii) = ((*D_)[ii] / (*Alpha_)[ii])*(1 + (*Beta_)[ii+1]);
1645 if (ii > 0) {
1646 DMT::Value(*APTAP_,ii-1,ii) = -(*D_)[ii]/(*Alpha_)[ii-1];
1647 DMT::Value(*APTAP_,ii,ii-1) = -(*D_)[ii]/(*Alpha_)[ii-1];
1648 }
1649 }
1650 DMT::SyncHostToDevice(*APTAP_);
1651
1652 DMT::SyncDeviceToHost(*L2_);
1653 for(int ii=0;ii<numBlocks_;ii++) {
1654 DMT::Value(*L2_,ii,ii) = 1./(*Alpha_)[ii];
1655 DMT::Value(*L2_,ii+1,ii) = -1./(*Alpha_)[ii];
1656 }
1657 DMT::SyncHostToDevice(*L2_);
1658
1659 DMT::SyncDeviceToHost(*Delta_);
1660 DMT::SyncDeviceToHost(*DeltaL2_);
1661 DMT::SyncDeviceToHost(*AU1TUDeltaL2_);
1662 DMT::SyncDeviceToHost(*AU1TAP_);
1663
1664 // M(end,1) = dold*(-Beta(1)/Alpha(1));
1665 // AU1TAP = Y'*[AU1TU*Delta*L2; M];
1666 //DeltaL2_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Delta_,*L2_,zero);
1667 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, recycleBlocks_, numBlocks_, numBlocks_+1,
1668 one, DMT::GetConstRawHostPtr(*Delta_), DMT::GetStride(*Delta_),
1669 DMT::GetConstRawHostPtr(*L2_), DMT::GetStride(*L2_),
1670 zero, DMT::GetRawHostPtr(*DeltaL2_), DMT::GetStride(*DeltaL2_));
1671 //AU1TUDeltaL2_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*AU1TU_,*DeltaL2_,zero);
1672 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, recycleBlocks_, numBlocks_, recycleBlocks_,
1673 one, DMT::GetConstRawHostPtr(*AU1TU_), DMT::GetStride(*AU1TU_),
1674 DMT::GetConstRawHostPtr(*DeltaL2_), DMT::GetStride(*DeltaL2_),
1675 zero, DMT::GetRawHostPtr(*AU1TUDeltaL2_), DMT::GetStride(*AU1TUDeltaL2_));
1676
1677 DMT::SyncDeviceToHost( *Y_);
1678 Teuchos::RCP<const DM> Y1 = DMT::SubviewConst( *Y_, recycleBlocks_, recycleBlocks_ );
1679 Teuchos::RCP<const DM> Y2 = DMT::SubviewConst( *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1680
1681 //AU1TAP_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Y1,*AU1TUDeltaL2_,zero);
1682 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, numBlocks_, recycleBlocks_,
1683 one, DMT::GetConstRawHostPtr(*Y1), DMT::GetStride(*Y1),
1684 DMT::GetConstRawHostPtr(*AU1TUDeltaL2_), DMT::GetStride(*AU1TUDeltaL2_),
1685 zero, DMT::GetRawHostPtr(*AU1TAP_), DMT::GetStride(*AU1TAP_));
1686 ScalarType val = dold * (-(*Beta_)[0]/(*Alpha_)[0]);
1687 for(int ii=0;ii<recycleBlocks_;ii++) {
1688 DMT::Value(*AU1TAP_,ii,0) += DMT::ValueConst(*Y2,numBlocks_-1,ii)*val;
1689 }
1690 DMT::SyncHostToDevice(*AU1TAP_);
1691
1692 // AU1TU = Y1'*AU1TU
1693 Teuchos::RCP<DM> Y1TAU1TU = DMT::Subview( *GY_, recycleBlocks_, recycleBlocks_ );
1694 //Y1TAU1TU->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Y1,*AU1TU_,zero);
1695 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_, recycleBlocks_,
1696 one, DMT::GetConstRawHostPtr(*Y1), DMT::GetStride(*Y1),
1697 DMT::GetConstRawHostPtr(*AU1TU_), DMT::GetStride(*AU1TU_),
1698 zero, DMT::GetRawHostPtr(*Y1TAU1TU), DMT::GetStride(*Y1TAU1TU));
1699 DMT::SyncHostToDevice(*GY_);
1700 DMT::Assign(*AU1TU_,*Y1TAU1TU);
1701
1702 // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
1703 DMT::PutScalar(*F_,zero);
1704 Teuchos::RCP<DM> F11 = DMT::Subview( *F_, recycleBlocks_, recycleBlocks_ );
1705 Teuchos::RCP<DM> F22 = DMT::Subview( *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1706 DMT::Assign(*F11,*AU1TU1_);
1707 DMT::SyncDeviceToHost(*F_);
1708 for(int ii=0;ii<numBlocks_;ii++) {
1709 DMT::Value(*F22,ii,ii) = (*D_)[ii];
1710 }
1711 DMT::SyncHostToDevice(*F_);
1712
1713 // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
1714 Teuchos::RCP<DM> G11 = DMT::Subview( *G_, recycleBlocks_, recycleBlocks_ );
1715 Teuchos::RCP<DM> G12 = DMT::Subview( *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1716 Teuchos::RCP<DM> G21 = DMT::Subview( *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1717 Teuchos::RCP<DM> G22 = DMT::Subview( *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1718 DMT::Assign(*G11,*AU1TAU1_);
1719 DMT::Assign(*G12,*AU1TAP_);
1720 DMT::Assign(*G22,*APTAP_);
1721 DMT::SyncDeviceToHost(*G_);
1722 // G21 = G12'; (no transpose operator exists for DM; Do copy manually)
1723 for (int ii=0;ii<recycleBlocks_;++ii)
1724 for (int jj=0;jj<numBlocks_;++jj)
1725 DMT::Value(*G21,jj,ii) = DMT::ValueConst(*G12,ii,jj);
1726 DMT::SyncHostToDevice(*G_);
1727
1728 // compute harmonic Ritz vectors
1729 getHarmonicVecs(*F_,*G_,*Y_);
1730 DMT::SyncHostToDevice(*Y_);
1731
1732 // U1 = [U1 P(:,2:end-1)]*Y;
1733 index.resize( numBlocks_ );
1734 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1735 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1736 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, rindex );
1737 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, rindex );
1738 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, rindex );
1739 MVT::MvTimesMatAddMv( one, *Ptmp, *Y2, zero, *PY2tmp );
1740 MVT::MvTimesMatAddMv( one, *U1tmp, *Y1, zero, *U1Y1tmp );
1741 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1742
1743 // Precompute some variables for next cycle
1744 DMT::SyncDeviceToHost(*GY_);
1745 DMT::SyncDeviceToHost(*AU1TAU1_);
1746 DMT::SyncDeviceToHost(*FY_);
1747 DMT::SyncDeviceToHost(*AU1TU1_);
1748
1749 // AU1TAU1 = Y'*G*Y;
1750 //GY_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*G_,*Y_,zero);
1751 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, numBlocks_+recycleBlocks_,
1752 recycleBlocks_, numBlocks_+recycleBlocks_,
1753 one, DMT::GetConstRawHostPtr(*G_), DMT::GetStride(*G_),
1754 DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1755 zero, DMT::GetRawHostPtr(*GY_), DMT::GetStride(*GY_));
1756 //AU1TAU1_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Y_,*GY_,zero);
1757 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_,
1758 numBlocks_+recycleBlocks_,
1759 one, DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1760 DMT::GetConstRawHostPtr(*GY_), DMT::GetStride(*GY_),
1761 zero, DMT::GetRawHostPtr(*AU1TAU1_), DMT::GetStride(*AU1TAU1_));
1762
1763 // AU1TU1 = Y'*F*Y;
1764 //FY_->multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*F_,*Y_,zero);
1765 blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, numBlocks_+recycleBlocks_,
1766 recycleBlocks_, numBlocks_+recycleBlocks_,
1767 one, DMT::GetConstRawHostPtr(*F_), DMT::GetStride(*F_),
1768 DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1769 zero, DMT::GetRawHostPtr(*FY_), DMT::GetStride(*FY_));
1770 //AU1TU1_->multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,*Y_,*FY_,zero);
1771 blas.GEMM( Teuchos::TRANS, Teuchos::NO_TRANS, recycleBlocks_, recycleBlocks_,
1772 numBlocks_+recycleBlocks_,
1773 one, DMT::GetConstRawHostPtr(*Y_), DMT::GetStride(*Y_),
1774 DMT::GetConstRawHostPtr(*FY_), DMT::GetStride(*FY_),
1775 zero, DMT::GetRawHostPtr(*AU1TU1_), DMT::GetStride(*AU1TU1_));
1776
1777 DMT::SyncHostToDevice(*GY_);
1778 DMT::SyncHostToDevice(*AU1TAU1_);
1779 DMT::SyncHostToDevice(*FY_);
1780 DMT::SyncHostToDevice(*AU1TU1_);
1781
1782 // dold = D(end);
1783 dold = (*D_)[numBlocks_-1];
1784
1785 // Indicate the size of the P, Beta structures generated this cycle
1786 lastp = numBlocks_+1;
1787 lastBeta = numBlocks_;
1788
1789 }
1790 }
1791 } // if (recycleBlocks_ > 0)
1792
1793 // Cleanup after end of cycle
1794
1795 // P = P(:,end-1:end);
1796 index.resize( 2 );
1797 index[0] = lastp-1; index[1] = lastp;
1798 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1799 index[0] = 0; index[1] = 1;
1800 MVT::SetBlock(*Ptmp2,index,*P_);
1801
1802 // Beta = Beta(end);
1803 (*Beta_)[0] = (*Beta_)[lastBeta];
1804
1805 // Delta = Delta(:,end);
1806 if (existU_) { // Delta only initialized if U exists
1807 Teuchos::RCP<DM> mu1 = DMT::Subview( *Delta_, recycleBlocks_, 1, 0, 0 );
1808 Teuchos::RCP<DM> mu2 = DMT::Subview( *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1809 DMT::Assign(*mu1,*mu2);
1810 }
1811
1812 // Now reinitialize state variables for next cycle
1813 newstate.P = Teuchos::null;
1814 index.resize( numBlocks_+1 );
1815 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1816 newstate.P = MVT::CloneViewNonConst( *P_, index );
1817
1818 newstate.Beta = Beta_;
1819 newstate.Beta_i = 1;
1820
1821 newstate.Delta = Delta_;
1822
1823 newstate.curDim = 1; // We have initialized the first search vector
1824
1825 // Pass to iteration object
1826 rcg_iter->initialize(newstate);
1827
1828 // increment cycle count
1829 cycle = cycle + 1;
1830
1831 }
1833 //
1834 // a debug status test (if any) stopped the iteration
1835 //
1837 else if (Teuchos::nonnull(debugStatusTest_) &&
1838 debugStatusTest_->getStatus() == Passed) {
1839 // we don't have convergence, but a debug test asked us to stop
1840 retType = Unconverged;
1841 isConverged = false;
1842 break; // break from while(1){rcg_iter->iterate()}
1843 }
1845 //
1846 // we returned from iterate(), but none of our status tests Passed.
1847 // something is wrong, and it is probably our fault.
1848 //
1850 else {
1851 retType = InconsistentState;
1852 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1853 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1854 }
1855 }
1856 catch (const StatusTestNaNError& e) {
1857 // A NaN was detected in the solver. Set the solution to zero and return unconverged.
1858 retType = NaNDetected;
1859 achievedTol_ = MT::one();
1860 Teuchos::RCP<MV> X = problem_->getLHS();
1861 MVT::MvInit( *X, SCT::zero() );
1862 printer_->stream(Warnings) << "Belos::RCGSolMgr::solve(): Warning! NaN has been detected!"
1863 << std::endl;
1864 return retType;
1865 }
1866 catch (const std::exception &e) {
1867 retType = NonspecificException;
1868 printer_->stream(Errors) << "Error! Caught std::exception in RCGIter::iterate() at iteration "
1869 << rcg_iter->getNumIters() << std::endl
1870 << e.what() << std::endl;
1871 throw;
1872 }
1873 }
1874
1875 // Inform the linear problem that we are finished with this block linear system.
1876 problem_->setCurrLS();
1877
1878 // Update indices for the linear systems to be solved.
1879 numRHS2Solve -= 1;
1880 if ( numRHS2Solve > 0 ) {
1881 currIdx[0]++;
1882 // Set the next indices.
1883 problem_->setLSIndex( currIdx );
1884 }
1885 else {
1886 currIdx.resize( numRHS2Solve );
1887 }
1888
1889 // Update the recycle space for the next linear system
1890 if (existU1_) { // be sure updated recycle space was created
1891 // U = U1
1892 MVT::SetBlock(*U1_,rindex,*U_);
1893 // Set flag indicating recycle space is now defined
1894 existU_ = true;
1895 if (numRHS2Solve > 0) { // also update AU, UTAU, and AUTAU
1896 // Free pointers in newstate
1897 newstate.P = Teuchos::null;
1898 newstate.Ap = Teuchos::null;
1899 newstate.r = Teuchos::null;
1900 newstate.z = Teuchos::null;
1901 newstate.U = Teuchos::null;
1902 newstate.AU = Teuchos::null;
1903 newstate.Alpha = Teuchos::null;
1904 newstate.Beta = Teuchos::null;
1905 newstate.Beta_i = 0;
1906 newstate.D = Teuchos::null;
1907 newstate.Delta = Teuchos::null;
1908 newstate.LUUTAU = Teuchos::null;
1909 newstate.ipiv = Teuchos::null;
1910 newstate.rTz_old = Teuchos::null;
1911
1912 // Reinitialize AU, UTAU, AUTAU
1913 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, rindex );
1914 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, rindex );
1915 // Initialize AU
1916 problem_->applyOp( *Utmp, *AUtmp );
1917 // Initialize UTAU
1918 MVT::MvTransMv( one, *Utmp, *AUtmp, *UTAU_ );
1919 // Initialize AUTAU ( AUTAU = AU'*(M\AU) )
1920 if ( precObj != Teuchos::null ) {
1921 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, rindex ); // use U1 as temp storage
1922 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1923 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, *AUTAU_ );
1924 } else {
1925 MVT::MvTransMv( one, *AUtmp, *AUtmp, *AUTAU_ );
1926 }
1927 } // if (numRHS2Solve > 0)
1928
1929 } // if (existU1)
1930 }// while ( numRHS2Solve > 0 )
1931
1932 }
1933
1934 // print final summary
1935 sTest_->print( printer_->stream(FinalSummary) );
1936
1937 // print timing information
1938#ifdef BELOS_TEUCHOS_TIME_MONITOR
1939 // Calling summarize() can be expensive, so don't call unless the
1940 // user wants to print out timing details. summarize() will do all
1941 // the work even if it's passed a "black hole" output stream.
1942 if (verbosity_ & TimingDetails)
1943 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1944#endif
1945
1946 // get iteration information for this solve
1947 numIters_ = maxIterTest_->getNumIters();
1948
1949 // Save the convergence test value ("achieved tolerance") for this solve.
1950 {
1951 using Teuchos::rcp_dynamic_cast;
1952 typedef StatusTestGenResNorm<ScalarType,MV,OP,DM> conv_test_type;
1953 // testValues is nonnull and not persistent.
1954 const std::vector<MagnitudeType>* pTestValues =
1955 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1956
1957 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1958 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1959 "method returned NULL. Please report this bug to the Belos developers.");
1960
1961 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1962 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1963 "method returned a vector of length zero. Please report this bug to the "
1964 "Belos developers.");
1965
1966 // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1967 // achieved tolerances for all vectors in the current solve(), or
1968 // just for the vectors from the last deflation?
1969 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1970 }
1971
1972 if (!isConverged) {
1973 return retType; // return from RCGSolMgr::solve()
1974 }
1975 return Converged; // return from RCGSolMgr::solve()
1976}
1977
1978// Compute the harmonic eigenpairs of the projected, dense system.
1979template<class ScalarType, class MV, class OP, class DM>
1980void RCGSolMgr<ScalarType,MV,OP,DM,true>::getHarmonicVecs(const DM& F,
1981 const DM& G,
1982 DM& Y ) {
1983
1984 // order of F,G
1985 int n = DMT::GetNumCols(F);
1986
1987 // The LAPACK interface
1988 Teuchos::LAPACK<int,ScalarType> lapack;
1989
1990 // Magnitude of harmonic Ritz values
1991 std::vector<MagnitudeType> w(n);
1992
1993 // Sorted order of harmonic Ritz values
1994 std::vector<int> iperm(n);
1995
1996 // Compute k smallest harmonic Ritz pairs
1997 // SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO )
1998 int itype = 1; // solve A*x = (lambda)*B*x
1999 char jobz='V'; // compute eigenvalues and eigenvectors
2000 char uplo='U'; // since F,G symmetric, reference only their upper triangular data
2001 std::vector<ScalarType> work(1);
2002 int lwork = -1;
2003 int info = 0;
2004 // since SYGV destroys workspace, create copies of F,G
2005 Teuchos::RCP<DM> F2 = DMT::CreateCopy( F );
2006 Teuchos::RCP<DM> G2 = DMT::CreateCopy( G );
2007
2008 DMT::SyncDeviceToHost(*F2);
2009 DMT::SyncDeviceToHost(*G2);
2010
2011 // query for optimal workspace size
2012 lapack.SYGV(itype, jobz, uplo, n, DMT::GetRawHostPtr(*G2), DMT::GetStride(*G2),
2013 DMT::GetRawHostPtr(*F2), DMT::GetStride(*F2), &w[0], &work[0], lwork, &info);
2014 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
2015 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
2016 lwork = (int)work[0];
2017 work.resize(lwork);
2018 lapack.SYGV(itype, jobz, uplo, n, DMT::GetRawHostPtr(*G2), DMT::GetStride(*G2),
2019 DMT::GetRawHostPtr(*F2), DMT::GetStride(*F2), &w[0], &work[0], lwork, &info);
2020 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
2021 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to compute eigensolutions.");
2022
2023
2024 // Construct magnitude of each harmonic Ritz value
2025 this->sort(w,n,iperm);
2026
2027 // Select recycledBlocks_ smallest eigenvectors
2028 for( int i=0; i<recycleBlocks_; i++ ) {
2029 for( int j=0; j<n; j++ ) {
2030 DMT::Value(Y,j,i) = DMT::ValueConst(*G2,j,iperm[i]);
2031 }
2032 }
2033
2034}
2035
2036// This method sorts list of n floating-point numbers and return permutation vector
2037template<class ScalarType, class MV, class OP, class DM>
2038void RCGSolMgr<ScalarType,MV,OP,DM,true>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm)
2039{
2040 int l, r, j, i, flag;
2041 int RR2;
2042 double dRR, dK;
2043
2044 // Initialize the permutation vector.
2045 for(j=0;j<n;j++)
2046 iperm[j] = j;
2047
2048 if (n <= 1) return;
2049
2050 l = n / 2 + 1;
2051 r = n - 1;
2052 l = l - 1;
2053 dRR = dlist[l - 1];
2054 dK = dlist[l - 1];
2055
2056 RR2 = iperm[l - 1];
2057 while (r != 0) {
2058 j = l;
2059 flag = 1;
2060
2061 while (flag == 1) {
2062 i = j;
2063 j = j + j;
2064
2065 if (j > r + 1)
2066 flag = 0;
2067 else {
2068 if (j < r + 1)
2069 if (dlist[j] > dlist[j - 1]) j = j + 1;
2070
2071 if (dlist[j - 1] > dK) {
2072 dlist[i - 1] = dlist[j - 1];
2073 iperm[i - 1] = iperm[j - 1];
2074 }
2075 else {
2076 flag = 0;
2077 }
2078 }
2079 }
2080 dlist[i - 1] = dRR;
2081 iperm[i - 1] = RR2;
2082 if (l == 1) {
2083 dRR = dlist [r];
2084 RR2 = iperm[r];
2085 dK = dlist[r];
2086 dlist[r] = dlist[0];
2087 iperm[r] = iperm[0];
2088 r = r - 1;
2089 }
2090 else {
2091 l = l - 1;
2092 dRR = dlist[l - 1];
2093 RR2 = iperm[l - 1];
2094 dK = dlist[l - 1];
2095 }
2096 }
2097 dlist[0] = dRR;
2098 iperm[0] = RR2;
2099}
2100
2101// This method requires the solver manager to return a std::string that describes itself.
2102template<class ScalarType, class MV, class OP, class DM>
2103std::string RCGSolMgr<ScalarType,MV,OP,DM,true>::description() const
2104{
2105 std::ostringstream oss;
2106 oss << "Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
2107 return oss.str();
2108}
2109
2110} // end Belos namespace
2111
2112#ifdef HAVE_BELOS_TPETRA
2113#include "BelosTpetraETIHelpers.hpp"
2114
2115#define BELOS_TPETRA_RCGSOLMGR_NOEXTERN_CALL(SC, LO, GO, NT) \
2116 BELOS_TPETRA_CALL(Belos::RCGSolMgr, SC, LO, GO, NT)
2117
2118#define BELOS_TPETRA_RCGSOLMGR_EXTERN_CALL(SC, LO, GO, NT) \
2119 BELOS_TPETRA_EXTERN_CALL(Belos::RCGSolMgr, SC, LO, GO, NT)
2120
2121TPETRA_INSTANTIATE_SLGN_NO_ORDINAL_SCALAR(BELOS_TPETRA_RCGSOLMGR_EXTERN_CALL)
2122#endif
2123
2124#endif /* BELOS_RCG_SOLMGR_HPP */
BelosConfigDefs.hpp
Belos header file which uses auto-configuration information to include necessary C++ headers.
BelosLinearProblem.hpp
Class which describes the linear problem to be solved by the iterative solver.
BelosOutputManager.hpp
Class which manages the output and verbosity of the Belos solvers.
BelosRCGIter.hpp
Belos concrete class for performing the RCG iteration.
BelosSolverManager.hpp
Pure virtual base class which describes the basic interface for a solver manager.
BelosStatusTestCombo.hpp
Belos::StatusTest for logically combining several status tests.
BelosStatusTestGenResNorm.hpp
Belos::StatusTestResNorm for specifying general residual norm stopping criteria.
BelosStatusTestMaxIters.hpp
Belos::StatusTest class for specifying a maximum number of iterations.
BelosStatusTestOutputFactory.hpp
A factory class for generating StatusTestOutput objects.
BelosTypes.hpp
Collection of types and exceptions used within the Belos solvers.
Belos::BelosError::BelosError
BelosError(const std::string &what_arg)
Definition BelosTypes.hpp:40
Belos::DenseMatTraits
Virtual base class which defines basic traits for the multi-vector type.
Definition BelosDenseMatTraits.hpp:52
Belos::DenseMatTraits::Subview
static Teuchos::RCP< DM > Subview(DM &source, int numRows, int numCols, int startRow=0, int startCol=0)
Marks host data modified to avoid device sync errors.
Definition BelosDenseMatTraits.hpp:117
Belos::DenseMatTraits::SyncHostToDevice
static void SyncHostToDevice(DM &dm)
Definition BelosDenseMatTraits.hpp:174
Belos::DenseMatTraits::PutScalar
static void PutScalar(DM &dm, ScalarType value=Teuchos::ScalarTraits< ScalarType >::zero())
Fill all entries with value. Value is zero if not specified.
Definition BelosDenseMatTraits.hpp:185
Belos::DenseMatTraits::Scale
static void Scale(DM &dm, ScalarType value)
Multiply all entries by a scalar. DM = value.*DM.
Definition BelosDenseMatTraits.hpp:189
Belos::DenseMatTraits::GetRawHostPtr
static ScalarType * GetRawHostPtr(DM &dm)
Returns a raw pointer to the (non-const) data on the host.
Definition BelosDenseMatTraits.hpp:95
Belos::DenseMatTraits::ValueConst
static const ScalarType & ValueConst(const DM &dm, const int i, const int j)
Access a const reference to the (i,j) entry of dm, e_i^T dm e_j.
Definition BelosDenseMatTraits.hpp:168
Belos::DenseMatTraits::SyncDeviceToHost
static void SyncDeviceToHost(DM &dm)
Definition BelosDenseMatTraits.hpp:171
Belos::DenseMatTraits::Assign
static void Assign(DM &thisDM, const DM &sourceDM)
Copies entries of sourceDM to thisDM (deep copy).
Definition BelosDenseMatTraits.hpp:198
Belos::DenseMatTraits::SubviewConst
static Teuchos::RCP< const DM > SubviewConst(const DM &source, int numRows, int numCols, int startRow=0, int startCol=0)
Definition BelosDenseMatTraits.hpp:120
Belos::DenseMatTraits::GetConstRawHostPtr
static ScalarType const * GetConstRawHostPtr(const DM &dm)
Returns a raw pointer to const data on the host.
Definition BelosDenseMatTraits.hpp:99
Belos::DenseMatTraits::Value
static ScalarType & Value(DM &dm, const int i, const int j)
Access a reference to the (i,j) entry of dm, e_i^T dm e_j.
Definition BelosDenseMatTraits.hpp:164
Belos::DenseMatTraits::GetStride
static int GetStride(const DM &dm)
Obtain the stride between the columns of dm.
Definition BelosDenseMatTraits.hpp:139
Belos::DenseMatTraits::Create
static Teuchos::RCP< DM > Create()
Creates a new empty DM with no dimension.
Definition BelosDenseMatTraits.hpp:61
Belos::Details::SolverManagerRequiresRealLapack
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
Definition BelosSolverManager.hpp:411
Belos::MultiVecTraits
Traits class which defines basic operations on multivectors.
Definition BelosMultiVecTraits.hpp:102
Belos::MultiVecTraits< ScalarType, MV, DM >::MvInit
static void MvInit(MV &mv, const ScalarType alpha=Teuchos::ScalarTraits< ScalarType >::zero())
Definition BelosMultiVecTraits.hpp:315
Belos::MultiVecTraits< ScalarType, MV, DM >::MvAddMv
static void MvAddMv(const ScalarType alpha, const MV &A, const ScalarType beta, const MV &B, MV &mv)
Definition BelosMultiVecTraits.hpp:236
Belos::MultiVecTraits< ScalarType, MV, DM >::MvTimesMatAddMv
static void MvTimesMatAddMv(const ScalarType alpha, const MV &A, const DM &B, const ScalarType beta, MV &mv)
Definition BelosMultiVecTraits.hpp:229
Belos::MultiVecTraits< ScalarType, MV, DM >::CloneView
static Teuchos::RCP< const MV > CloneView(const MV &mv, const std::vector< int > &index)
Definition BelosMultiVecTraits.hpp:165
Belos::MultiVecTraits< ScalarType, MV, DM >::MvTransMv
static void MvTransMv(const ScalarType alpha, const MV &A, const MV &mv, DM &B)
Definition BelosMultiVecTraits.hpp:253
Belos::MultiVecTraits< ScalarType, MV, DM >::MvDot
static void MvDot(const MV &mv, const MV &A, std::vector< ScalarType > &b)
Definition BelosMultiVecTraits.hpp:258
Belos::MultiVecTraits< ScalarType, MV, DM >::SetBlock
static void SetBlock(const MV &A, const std::vector< int > &index, MV &mv)
Definition BelosMultiVecTraits.hpp:284
Belos::MultiVecTraits< ScalarType, MV, DM >::CloneViewNonConst
static Teuchos::RCP< MV > CloneViewNonConst(MV &mv, const std::vector< int > &index)
Definition BelosMultiVecTraits.hpp:146
Belos::OperatorTraits
Class which defines basic traits for the operator type.
Definition BelosOperatorTraits.hpp:78
Belos::OperatorTraits::Apply
static void Apply(const OP &Op, const MV &x, MV &y, ETrans trans=NOTRANS)
Apply Op to x, putting the result into y.
Definition BelosOperatorTraits.hpp:91
Belos::OutputManager
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
Definition BelosOutputManager.hpp:42
Belos::RCGIter
This class implements the RCG iteration, where a single-stdvector Krylov subspace is constructed.
Definition BelosRCGIter.hpp:110
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::setProblem
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP, DM > > &problem) override
Set the linear problem that needs to be solved.
Definition BelosRCGSolMgr.hpp:246
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::getTimers
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Definition BelosRCGSolMgr.hpp:220
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::setDebugStatusTest
void setDebugStatusTest(const Teuchos::RCP< StatusTest< ScalarType, MV, OP, DM > > &debugStatusTest) override
Set a debug status test that will be OR-combined into the top-level status test.
Definition BelosRCGSolMgr.hpp:252
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::RCGSolMgr
RCGSolMgr()
Empty constructor for RCGSolMgr. This constructor takes no arguments and sets the default values for ...
Definition BelosRCGSolMgr.hpp:453
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::setParameters
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params) override
Set the parameters the solver manager should use to solve the linear problem.
Definition BelosRCGSolMgr.hpp:531
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::reset
void reset(const ResetType type) override
Performs a reset of the solver manager specified by the ResetType. This informs the solver manager th...
Definition BelosRCGSolMgr.hpp:273
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::getCurrentParameters
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
Definition BelosRCGSolMgr.hpp:213
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::getNumIters
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Definition BelosRCGSolMgr.hpp:233
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::isLOADetected
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve.
Definition BelosRCGSolMgr.hpp:238
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::getValidParameters
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const override
Get a parameter list containing the valid parameters for this object.
Definition BelosRCGSolMgr.hpp:717
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::getProblem
const LinearProblem< ScalarType, MV, OP, DM > & getProblem() const override
Return a reference to the linear problem being solved by this solver manager.
Definition BelosRCGSolMgr.hpp:205
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::achievedTol
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
Definition BelosRCGSolMgr.hpp:228
Belos::RCGSolMgr< ScalarType, MV, OP, DM, true >::clone
Teuchos::RCP< SolverManager< ScalarType, MV, OP, DM > > clone() const override
clone for Inverted Injection (DII)
Definition BelosRCGSolMgr.hpp:197
Belos::RCGSolMgr
Implementation of the RCG (Recycling Conjugate Gradient) iterative linear solver.
Definition BelosRCGSolMgr.hpp:120
Belos::RCGSolMgr::RCGSolMgr
RCGSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP, DM > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
Definition BelosRCGSolMgr.hpp:131
Belos::RCGSolMgr::clone
Teuchos::RCP< SolverManager< ScalarType, MV, OP, DM > > clone() const override
clone for Inverted Injection (DII)
Definition BelosRCGSolMgr.hpp:138
Belos::RCGSolMgrLAPACKFailure
RCGSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
Definition BelosRCGSolMgr.hpp:103
Belos::RCGSolMgrLAPACKFailure::RCGSolMgrLAPACKFailure
RCGSolMgrLAPACKFailure(const std::string &what_arg)
Definition BelosRCGSolMgr.hpp:104
Belos::RCGSolMgrLinearProblemFailure
RCGSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i.e.
Definition BelosRCGSolMgr.hpp:93
Belos::RCGSolMgrLinearProblemFailure::RCGSolMgrLinearProblemFailure
RCGSolMgrLinearProblemFailure(const std::string &what_arg)
Definition BelosRCGSolMgr.hpp:94
Belos::StatusTestCombo
A class for extending the status testing capabilities of Belos via logical combinations.
Definition BelosStatusTestCombo.hpp:59
Belos::StatusTestGenResNorm
An implementation of StatusTestResNorm using a family of residual norms.
Definition BelosStatusTestGenResNorm.hpp:47
Belos::StatusTestGenResNorm::getTestValue
const std::vector< MagnitudeType > * getTestValue() const
Returns the test value, , computed in most recent call to CheckStatus.
Definition BelosStatusTestGenResNorm.hpp:199
Belos::StatusTest
A pure virtual class for defining the status tests for the Belos iterative solvers.
Definition BelosStatusTest.hpp:50
Belos::StatusTestMaxIters
A Belos::StatusTest class for specifying a maximum number of iterations.
Definition BelosStatusTestMaxIters.hpp:31
Belos::StatusTestOutputFactory
A factory class for generating StatusTestOutput objects.
Definition BelosStatusTestOutputFactory.hpp:37
Belos::StatusTestOutputFactory::create
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP, DM > > create(const Teuchos::RCP< OutputManager< ScalarType > > &printer, Teuchos::RCP< StatusTest< ScalarType, MV, OP, DM > > test, int mod, int printStates)
Create the StatusTestOutput object specified by the outputStyle.
Definition BelosStatusTestOutputFactory.hpp:78
Belos::StatusTestDetails
@ StatusTestDetails
Definition BelosTypes.hpp:292
Belos::FinalSummary
@ FinalSummary
Definition BelosTypes.hpp:290
Belos::TimingDetails
@ TimingDetails
Definition BelosTypes.hpp:291
Belos::ReturnType
ReturnType
Whether the Belos solve converged for all linear systems.
Definition BelosTypes.hpp:133
Belos::NaNDetected
@ NaNDetected
Definition BelosTypes.hpp:139
Belos::Unconverged
@ Unconverged
Definition BelosTypes.hpp:135
Belos::MaxItersReached
@ MaxItersReached
Definition BelosTypes.hpp:136
Belos::NonspecificException
@ NonspecificException
Definition BelosTypes.hpp:142
Belos::InconsistentState
@ InconsistentState
Definition BelosTypes.hpp:143
Belos::Undetermined
@ Undetermined
Definition BelosTypes.hpp:144
Belos::ResetType
ResetType
How to reset the solver.
Definition BelosTypes.hpp:237
Belos::RecycleSubspace
@ RecycleSubspace
Definition BelosTypes.hpp:238
Belos::DefaultSolverParameters::convTol
static const double convTol
Default convergence tolerance.
Definition BelosTypes.hpp:324
Belos::RCGIterState
Structure to contain pointers to RCGIter state variables.
Definition BelosRCGIter.hpp:55
Belos::RCGIterState::curDim
int curDim
The current dimension of the reduction.
Definition BelosRCGIter.hpp:60
Belos::RCGIterState::AU
Teuchos::RCP< MV > AU
Definition BelosRCGIter.hpp:78
Belos::RCGIterState::r
Teuchos::RCP< MV > r
The current residual.
Definition BelosRCGIter.hpp:69
Belos::RCGIterState::existU
bool existU
Flag to indicate the recycle space should be used.
Definition BelosRCGIter.hpp:75
Belos::RCGIterState::P
Teuchos::RCP< MV > P
The current Krylov basis.
Definition BelosRCGIter.hpp:63
Belos::RCGIterState::U
Teuchos::RCP< MV > U
The recycled subspace and its image.
Definition BelosRCGIter.hpp:78
Belos::RCGIterState::D
Teuchos::RCP< std::vector< ScalarType > > D
Definition BelosRCGIter.hpp:85
Belos::RCGIterState::Beta
Teuchos::RCP< std::vector< ScalarType > > Beta
Definition BelosRCGIter.hpp:83
Belos::RCGIterState::rTz_old
Teuchos::RCP< std::vector< ScalarType > > rTz_old
Definition BelosRCGIter.hpp:86
Belos::RCGIterState::Alpha
Teuchos::RCP< std::vector< ScalarType > > Alpha
Coefficients arising in RCG iteration.
Definition BelosRCGIter.hpp:82
Belos::RCGIterState::Delta
Teuchos::RCP< DM > Delta
Solutions to local least-squares problems.
Definition BelosRCGIter.hpp:89
Belos::RCGIterState::ipiv
Teuchos::RCP< std::vector< int > > ipiv
Data from LU factorization of U^T A U.
Definition BelosRCGIter.hpp:94
Belos::RCGIterState::z
Teuchos::RCP< MV > z
The current preconditioned residual.
Definition BelosRCGIter.hpp:72
Belos::RCGIterState::LUUTAU
Teuchos::RCP< DM > LUUTAU
The LU factorization of the matrix U^T A U.
Definition BelosRCGIter.hpp:92
Belos::RCGIterState::Ap
Teuchos::RCP< MV > Ap
A times current search vector.
Definition BelosRCGIter.hpp:66
Generated for Belos by doxygen 1.15.0