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