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
1029 Teuchos::LAPACK<int,ScalarType> lapack;
1030 std::vector<int> index(recycleBlocks_);
1031 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1032 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1033
1034 // Count of number of cycles performed on current rhs
1035 int cycle = 0;
1036
1037 // Set the current parameters if they were not set before.
1038 // NOTE: This may occur if the user generated the solver manager with the default constructor and
1039 // then didn't set any parameters using setParameters().
1040 if (!params_Set_) {
1041 setParameters(Teuchos::parameterList(*getValidParameters()));
1042 }
1043
1044 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,RCGSolMgrLinearProblemFailure,
1045 "Belos::RCGSolMgr::solve(): Linear problem is not a valid object.");
1046 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),RCGSolMgrLinearProblemFailure,
1047 "Belos::RCGSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1048 TEUCHOS_TEST_FOR_EXCEPTION((problem_->getLeftPrec() != Teuchos::null)&&(problem_->getRightPrec() != Teuchos::null),
1049 RCGSolMgrLinearProblemFailure,
1050 "Belos::RCGSolMgr::solve(): RCG does not support split preconditioning, only set left or right preconditioner.");
1051
1052 // Grab the preconditioning object
1053 Teuchos::RCP<OP> precObj;
1054 if (problem_->getLeftPrec() != Teuchos::null) {
1055 precObj = Teuchos::rcp_const_cast<OP>(problem_->getLeftPrec());
1056 }
1057 else if (problem_->getRightPrec() != Teuchos::null) {
1058 precObj = Teuchos::rcp_const_cast<OP>(problem_->getRightPrec());
1059 }
1060
1061 // Create indices for the linear systems to be solved.
1062 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1063 std::vector<int> currIdx(1);
1064 currIdx[0] = 0;
1065
1066 // Inform the linear problem of the current linear system to solve.
1067 problem_->setLSIndex( currIdx );
1068
1069 // Check the number of blocks and change them if necessary.
1070 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1071 if (numBlocks_ > dim) {
1072 numBlocks_ = Teuchos::asSafe<int>(dim);
1073 params_->set("Num Blocks", numBlocks_);
1074 printer_->stream(Warnings) <<
1075 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1076 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1077 }
1078
1079 // Initialize storage for all state variables
1080 initializeStateStorage();
1081
1082 // Parameter list
1083 Teuchos::ParameterList plist;
1084 plist.set("Num Blocks",numBlocks_);
1085 plist.set("Recycled Blocks",recycleBlocks_);
1086
1087 // Reset the status test.
1088 outputTest_->reset();
1089
1090 // Assume convergence is achieved, then let any failed convergence set this to false.
1091 bool isConverged = true;
1092
1093 // Compute AU = A*U, UTAU = U'*AU, AUTAU = (AU)'*(AU)
1094 if (existU_) {
1095 index.resize(recycleBlocks_);
1096 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1097 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1098 index.resize(recycleBlocks_);
1099 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1100 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1101 // Initialize AU
1102 problem_->applyOp( *Utmp, *AUtmp );
1103 // Initialize UTAU
1104 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1105 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1106 // Initialize AUTAU ( AUTAU = AU'*(M\AU) )
1107 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1108 if ( precObj != Teuchos::null ) {
1109 index.resize(recycleBlocks_);
1110 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1111 index.resize(recycleBlocks_);
1112 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1113 Teuchos::RCP<MV> PCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage
1114 OPT::Apply( *precObj, *AUtmp, *PCAU );
1115 MVT::MvTransMv( one, *AUtmp, *PCAU, AUTAUtmp );
1116 } else {
1117 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1118 }
1119 }
1120
1122 // RCG solver
1123
1124 Teuchos::RCP<RCGIter<ScalarType,MV,OP> > rcg_iter;
1125 rcg_iter = Teuchos::rcp( new RCGIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,plist) );
1126
1127 // Enter solve() iterations
1128 {
1129#ifdef BELOS_TEUCHOS_TIME_MONITOR
1130 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1131#endif
1132
1133 while ( numRHS2Solve > 0 ) {
1134
1135 // Debugging output to tell use if recycle space exists and will be used
1136 if (printer_->isVerbosity( Debug ) ) {
1137 if (existU_) printer_->print( Debug, "Using recycle space generated from previous call to solve()." );
1138 else printer_->print( Debug, "No recycle space exists." );
1139 }
1140
1141 // Reset the number of iterations.
1142 rcg_iter->resetNumIters();
1143
1144 // Set the current number of recycle blocks and subspace dimension with the RCG iteration.
1145 rcg_iter->setSize( recycleBlocks_, numBlocks_ );
1146
1147 // Reset the number of calls that the status test output knows about.
1148 outputTest_->resetNumCalls();
1149
1150 // indicate that updated recycle space has not yet been generated for this linear system
1151 existU1_ = false;
1152
1153 // reset cycle count
1154 cycle = 0;
1155
1156 // Get the current residual
1157 problem_->computeCurrResVec( &*r_ );
1158
1159 // If U exists, find best soln over this space first
1160 if (existU_) {
1161 // Solve linear system UTAU * y = (U'*r)
1162 Teuchos::SerialDenseMatrix<int,ScalarType> Utr(recycleBlocks_,1);
1163 index.resize(recycleBlocks_);
1164 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1165 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1166 MVT::MvTransMv( one, *Utmp, *r_, Utr );
1167 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1168 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1169 LUUTAUtmp.assign(UTAUtmp);
1170 int info = 0;
1171 lapack.GESV(recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], Utr.values(), Utr.stride(), &info);
1172 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1173 "Belos::RCGSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1174
1175 // Update solution (x = x + U*y)
1176 MVT::MvTimesMatAddMv( one, *Utmp, Utr, one, *problem_->getCurrLHSVec() );
1177
1178 // Update residual ( r = r - AU*y )
1179 index.resize(recycleBlocks_);
1180 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1181 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1182 MVT::MvTimesMatAddMv( -one, *AUtmp, Utr, one, *r_ );
1183 }
1184
1185 if ( precObj != Teuchos::null ) {
1186 OPT::Apply( *precObj, *r_, *z_ );
1187 } else {
1188 z_ = r_;
1189 }
1190
1191 // rTz_old = r'*z
1192 MVT::MvTransMv( one, *r_, *z_, *rTz_old_ );
1193
1194 if ( existU_ ) {
1195 // mu = UTAU\‍(AU'*z);
1196 Teuchos::SerialDenseMatrix<int,ScalarType> mu( Teuchos::View, *Delta_, recycleBlocks_, 1);
1197 index.resize(recycleBlocks_);
1198 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1199 Teuchos::RCP<const MV> AUtmp = MVT::CloneView( *AU_, index );
1200 MVT::MvTransMv( one, *AUtmp, *z_, mu );
1201 Teuchos::SerialDenseMatrix<int,ScalarType> LUUTAUtmp( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ );
1202 char TRANS = 'N';
1203 int info;
1204 lapack.GETRS( TRANS, recycleBlocks_, 1, LUUTAUtmp.values(), LUUTAUtmp.stride(), &(*ipiv_)[0], mu.values(), mu.stride(), &info );
1205 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1206 "Belos::RCGSolMgr::solve(): LAPACK GETRS failed to compute a solution.");
1207 // p = z - U*mu;
1208 index.resize( 1 );
1209 index[0] = 0;
1210 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1211 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1212 MVT::MvTimesMatAddMv( -one, *U_, mu, one, *Ptmp );
1213 } else {
1214 // p = z;
1215 index.resize( 1 );
1216 index[0] = 0;
1217 Teuchos::RCP<MV> Ptmp = MVT::CloneViewNonConst( *P_, index );
1218 MVT::MvAddMv(one,*z_,zero,*z_,*Ptmp);
1219 }
1220
1221 // Set the new state and initialize the solver.
1222 RCGIterState<ScalarType,MV> newstate;
1223
1224 // Create RCP views here
1225 index.resize( numBlocks_+1 );
1226 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1227 newstate.P = MVT::CloneViewNonConst( *P_, index );
1228 index.resize( recycleBlocks_ );
1229 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1230 newstate.U = MVT::CloneViewNonConst( *U_, index );
1231 index.resize( recycleBlocks_ );
1232 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1233 newstate.AU = MVT::CloneViewNonConst( *AU_, index );
1234 newstate.Alpha = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Alpha_, numBlocks_, 1 ) );
1235 newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1 ) );
1236 newstate.D = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *D_, numBlocks_, 1 ) );
1237 newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1238 newstate.LUUTAU = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *LUUTAU_, recycleBlocks_, recycleBlocks_ ) );
1239 // assign the rest of the values in the struct
1240 newstate.curDim = 1; // We have initialized the first search vector
1241 newstate.Ap = Ap_;
1242 newstate.r = r_;
1243 newstate.z = z_;
1244 newstate.existU = existU_;
1245 newstate.ipiv = ipiv_;
1246 newstate.rTz_old = rTz_old_;
1247
1248 rcg_iter->initialize(newstate);
1249
1250 while(1) {
1251
1252 // tell rcg_iter to iterate
1253 try {
1254 rcg_iter->iterate();
1255
1257 //
1258 // check convergence first
1259 //
1261 if ( convTest_->getStatus() == Passed ) {
1262 // We have convergence
1263 break; // break from while(1){rcg_iter->iterate()}
1264 }
1266 //
1267 // check for maximum iterations
1268 //
1270 else if ( maxIterTest_->getStatus() == Passed ) {
1271 // we don't have convergence
1272 isConverged = false;
1273 break; // break from while(1){rcg_iter->iterate()}
1274 }
1276 //
1277 // check if cycle complete; update for next cycle
1278 //
1280 else if ( rcg_iter->getCurSubspaceDim() == rcg_iter->getMaxSubspaceDim() ) {
1281 // index into P_ of last search vector generated this cycle
1282 int lastp = -1;
1283 // index into Beta_ of last entry generated this cycle
1284 int lastBeta = -1;
1285 if (recycleBlocks_ > 0) {
1286 if (!existU_) {
1287 if (cycle == 0) { // No U, no U1
1288
1289 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, numBlocks_, numBlocks_ );
1290 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, numBlocks_, numBlocks_ );
1291 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1292 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1293 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1294 Ftmp.putScalar(zero);
1295 Gtmp.putScalar(zero);
1296 for (int ii=0;ii<numBlocks_;ii++) {
1297 Gtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1298 if (ii > 0) {
1299 Gtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1300 Gtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1301 }
1302 Ftmp(ii,ii) = Dtmp(ii,0);
1303 }
1304
1305 // compute harmonic Ritz vectors
1306 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, numBlocks_, recycleBlocks_ );
1307 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1308
1309 // U1 = [P(:,1:end-1)*Y];
1310 index.resize( numBlocks_ );
1311 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1312 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1313 index.resize( recycleBlocks_ );
1314 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1315 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1316 MVT::MvTimesMatAddMv( one, *Ptmp, Ytmp, zero, *U1tmp );
1317
1318 // Precompute some variables for next cycle
1319
1320 // AU1TAU1 = Y'*G*Y;
1321 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, numBlocks_, recycleBlocks_ );
1322 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1323 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1324 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1325
1326 // AU1TU1 = Y'*F*Y;
1327 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, numBlocks_, recycleBlocks_ );
1328 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1329 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1330 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1331
1332 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, 1 );
1333 // Must reinitialize AU1TAP; can become dense later
1334 AU1TAPtmp.putScalar(zero);
1335 // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
1336 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1337 for (int ii=0; ii<recycleBlocks_; ++ii) {
1338 AU1TAPtmp(ii,0) = Ytmp(numBlocks_-1,ii) * alphatmp;
1339 }
1340
1341 // indicate that updated recycle space now defined
1342 existU1_ = true;
1343
1344 // Indicate the size of the P, Beta structures generated this cycle
1345 lastp = numBlocks_;
1346 lastBeta = numBlocks_-1;
1347
1348 } // if (cycle == 0)
1349 else { // No U, but U1 guaranteed to exist now
1350
1351 // Finish computation of subblocks
1352 // AU1TAP = AU1TAP * D(1);
1353 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1354 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1355 AU1TAPtmp.scale(Dtmp(0,0));
1356
1357 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1358 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1359 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1360 APTAPtmp.putScalar(zero);
1361 for (int ii=0; ii<numBlocks_; ii++) {
1362 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1363 if (ii > 0) {
1364 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1365 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1366 }
1367 }
1368
1369 // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
1370 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1371 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1372 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1373 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1374 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1375 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1376 F11.assign(AU1TU1tmp);
1377 F12.putScalar(zero);
1378 F21.putScalar(zero);
1379 F22.putScalar(zero);
1380 for(int ii=0;ii<numBlocks_;ii++) {
1381 F22(ii,ii) = Dtmp(ii,0);
1382 }
1383
1384 // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
1385 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1386 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1387 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1388 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1389 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1390 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1391 G11.assign(AU1TAU1tmp);
1392 G12.assign(AU1TAPtmp);
1393 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1394 for (int ii=0;ii<recycleBlocks_;++ii)
1395 for (int jj=0;jj<numBlocks_;++jj)
1396 G21(jj,ii) = G12(ii,jj);
1397 G22.assign(APTAPtmp);
1398
1399 // compute harmonic Ritz vectors
1400 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1401 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1402
1403 // U1 = [U1 P(:,2:end-1)]*Y;
1404 index.resize( numBlocks_ );
1405 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1406 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1407 index.resize( recycleBlocks_ );
1408 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1409 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1410 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1411 index.resize( recycleBlocks_ );
1412 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1413 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1414 index.resize( recycleBlocks_ );
1415 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1416 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1417 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1418 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1419 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1420 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1421
1422 // Precompute some variables for next cycle
1423
1424 // AU1TAU1 = Y'*G*Y;
1425 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1426 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1427 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1428
1429 // AU1TAP = zeros(k,m);
1430 // AU1TAP(:,1) = Y(end,:)' * (-1/Alpha(end));
1431 AU1TAPtmp.putScalar(zero);
1432 ScalarType alphatmp = -1.0 / Alphatmp(numBlocks_-1,0);
1433 for (int ii=0; ii<recycleBlocks_; ++ii) {
1434 AU1TAPtmp(ii,0) = Ytmp(numBlocks_+recycleBlocks_-1,ii) * alphatmp;
1435 }
1436
1437 // AU1TU1 = Y'*F*Y;
1438 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1439 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1440 //Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1441 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1442
1443 // Indicate the size of the P, Beta structures generated this cycle
1444 lastp = numBlocks_+1;
1445 lastBeta = numBlocks_;
1446
1447 } // if (cycle != 1)
1448 } // if (!existU_)
1449 else { // U exists
1450 if (cycle == 0) { // No U1, but U exists
1451 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1452 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_, 1 );
1453 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1454 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1455 APTAPtmp.putScalar(zero);
1456 for (int ii=0; ii<numBlocks_; ii++) {
1457 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii,0));
1458 if (ii > 0) {
1459 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1460 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1461 }
1462 }
1463
1464 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1465 L2tmp.putScalar(zero);
1466 for(int ii=0;ii<numBlocks_;ii++) {
1467 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1468 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1469 }
1470
1471 // AUTAP = UTAU*Delta*L2;
1472 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAPtmp( Teuchos::View, *AUTAP_, recycleBlocks_, numBlocks_ );
1473 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1474 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1475 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2tmp( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1476 DeltaL2tmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1477 AUTAPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,UTAUtmp,DeltaL2tmp,zero);
1478
1479 // F = [UTAU zeros(k,m); zeros(m,k) diag(D)];
1480 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1481 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1482 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1483 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1484 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1485 F11.assign(UTAUtmp);
1486 F12.putScalar(zero);
1487 F21.putScalar(zero);
1488 F22.putScalar(zero);
1489 for(int ii=0;ii<numBlocks_;ii++) {
1490 F22(ii,ii) = Dtmp(ii,0);
1491 }
1492
1493 // G = [AUTAU AUTAP; AUTAP' APTAP];
1494 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1495 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1496 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1497 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1498 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1499 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1500 G11.assign(AUTAUtmp);
1501 G12.assign(AUTAPtmp);
1502 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1503 for (int ii=0;ii<recycleBlocks_;++ii)
1504 for (int jj=0;jj<numBlocks_;++jj)
1505 G21(jj,ii) = G12(ii,jj);
1506 G22.assign(APTAPtmp);
1507
1508 // compute harmonic Ritz vectors
1509 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1510 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1511
1512 // U1 = [U P(:,1:end-1)]*Y;
1513 index.resize( recycleBlocks_ );
1514 for (int ii=0; ii<(recycleBlocks_); ++ii) { index[ii] = ii; }
1515 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1516 index.resize( numBlocks_ );
1517 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii; }
1518 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1519 index.resize( recycleBlocks_ );
1520 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1521 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1522 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1523 index.resize( recycleBlocks_ );
1524 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1525 Teuchos::RCP<MV> UY1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1526 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1527 index.resize( recycleBlocks_ );
1528 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1529 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1530 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1531 MVT::MvTimesMatAddMv( one, *Utmp, Y1, zero, *UY1tmp );
1532 MVT::MvAddMv(one,*UY1tmp, one, *PY2tmp, *U1tmp);
1533
1534 // Precompute some variables for next cycle
1535
1536 // AU1TAU1 = Y'*G*Y;
1537 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1538 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1539 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1540 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1541
1542 // AU1TU1 = Y'*F*Y;
1543 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1544 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1545 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1546 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1547
1548 // AU1TU = UTAU;
1549 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1550 AU1TUtmp.assign(UTAUtmp);
1551
1552 // dold = D(end);
1553 dold = Dtmp(numBlocks_-1,0);
1554
1555 // indicate that updated recycle space now defined
1556 existU1_ = true;
1557
1558 // Indicate the size of the P, Beta structures generated this cycle
1559 lastp = numBlocks_;
1560 lastBeta = numBlocks_-1;
1561 }
1562 else { // Have U and U1
1563 Teuchos::SerialDenseMatrix<int,ScalarType> Alphatmp( Teuchos::View, *Alpha_, numBlocks_, 1 );
1564 Teuchos::SerialDenseMatrix<int,ScalarType> Betatmp( Teuchos::View, *Beta_, numBlocks_+1, 1 );
1565 Teuchos::SerialDenseMatrix<int,ScalarType> Dtmp( Teuchos::View, *D_, numBlocks_, 1 );
1566 Teuchos::SerialDenseMatrix<int,ScalarType> APTAPtmp( Teuchos::View, *APTAP_, numBlocks_, numBlocks_ );
1567 for (int ii=0; ii<numBlocks_; ii++) {
1568 APTAPtmp(ii,ii) = (Dtmp(ii,0) / Alphatmp(ii,0))*(1 + Betatmp(ii+1,0));
1569 if (ii > 0) {
1570 APTAPtmp(ii-1,ii) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1571 APTAPtmp(ii,ii-1) = -Dtmp(ii,0)/Alphatmp(ii-1,0);
1572 }
1573 }
1574
1575 Teuchos::SerialDenseMatrix<int,ScalarType> L2tmp( Teuchos::View, *L2_, numBlocks_+1, numBlocks_ );
1576 for(int ii=0;ii<numBlocks_;ii++) {
1577 L2tmp(ii,ii) = 1./Alphatmp(ii,0);
1578 L2tmp(ii+1,ii) = -1./Alphatmp(ii,0);
1579 }
1580
1581 // M(end,1) = dold*(-Beta(1)/Alpha(1));
1582 // AU1TAP = Y'*[AU1TU*Delta*L2; M];
1583 Teuchos::SerialDenseMatrix<int,ScalarType> DeltaL2( Teuchos::View, *DeltaL2_, recycleBlocks_, numBlocks_ );
1584 Teuchos::SerialDenseMatrix<int,ScalarType> Deltatmp( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_+1 );
1585 DeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Deltatmp,L2tmp,zero);
1586 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUDeltaL2( Teuchos::View, *AU1TUDeltaL2_, recycleBlocks_, numBlocks_ );
1587 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TUtmp( Teuchos::View, *AU1TU_, recycleBlocks_, recycleBlocks_ );
1588 AU1TUDeltaL2.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,AU1TUtmp,DeltaL2,zero);
1589 Teuchos::SerialDenseMatrix<int,ScalarType> Y1( Teuchos::View, *Y_, recycleBlocks_, recycleBlocks_ );
1590 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAPtmp( Teuchos::View, *AU1TAP_, recycleBlocks_, numBlocks_ );
1591 AU1TAPtmp.putScalar(zero);
1592 AU1TAPtmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUDeltaL2,zero);
1593 Teuchos::SerialDenseMatrix<int,ScalarType> Y2( Teuchos::View, *Y_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1594 ScalarType val = dold * (-Betatmp(0,0)/Alphatmp(0,0));
1595 for(int ii=0;ii<recycleBlocks_;ii++) {
1596 AU1TAPtmp(ii,0) += Y2(numBlocks_-1,ii)*val;
1597 }
1598
1599 // AU1TU = Y1'*AU1TU
1600 Teuchos::SerialDenseMatrix<int,ScalarType> Y1TAU1TU( Teuchos::View, *GY_, recycleBlocks_, recycleBlocks_ );
1601 Y1TAU1TU.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Y1,AU1TUtmp,zero);
1602 AU1TUtmp.assign(Y1TAU1TU);
1603
1604 // F = [AU1TU1 zeros(k,m); zeros(m,k) diag(D)];
1605 Teuchos::SerialDenseMatrix<int,ScalarType> Ftmp( Teuchos::View, *F_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1606 Teuchos::SerialDenseMatrix<int,ScalarType> F11( Teuchos::View, *F_, recycleBlocks_, recycleBlocks_ );
1607 Teuchos::SerialDenseMatrix<int,ScalarType> F12( Teuchos::View, *F_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1608 Teuchos::SerialDenseMatrix<int,ScalarType> F21( Teuchos::View, *F_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1609 Teuchos::SerialDenseMatrix<int,ScalarType> F22( Teuchos::View, *F_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1610 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TU1tmp( Teuchos::View, *AU1TU1_, recycleBlocks_, recycleBlocks_ );
1611 F11.assign(AU1TU1tmp);
1612 for(int ii=0;ii<numBlocks_;ii++) {
1613 F22(ii,ii) = Dtmp(ii,0);
1614 }
1615
1616 // G = [AU1TAU1 AU1TAP; AU1TAP' APTAP];
1617 Teuchos::SerialDenseMatrix<int,ScalarType> Gtmp( Teuchos::View, *G_, (numBlocks_+recycleBlocks_), (numBlocks_+recycleBlocks_) );
1618 Teuchos::SerialDenseMatrix<int,ScalarType> G11( Teuchos::View, *G_, recycleBlocks_, recycleBlocks_ );
1619 Teuchos::SerialDenseMatrix<int,ScalarType> G12( Teuchos::View, *G_, recycleBlocks_, numBlocks_, 0, recycleBlocks_ );
1620 Teuchos::SerialDenseMatrix<int,ScalarType> G21( Teuchos::View, *G_, numBlocks_, recycleBlocks_, recycleBlocks_, 0 );
1621 Teuchos::SerialDenseMatrix<int,ScalarType> G22( Teuchos::View, *G_, numBlocks_, numBlocks_, recycleBlocks_, recycleBlocks_ );
1622 Teuchos::SerialDenseMatrix<int,ScalarType> AU1TAU1tmp( Teuchos::View, *AU1TAU1_, recycleBlocks_, recycleBlocks_ );
1623 G11.assign(AU1TAU1tmp);
1624 G12.assign(AU1TAPtmp);
1625 // G21 = G12'; (no transpose operator exists for SerialDenseMatrix; Do copy manually)
1626 for (int ii=0;ii<recycleBlocks_;++ii)
1627 for (int jj=0;jj<numBlocks_;++jj)
1628 G21(jj,ii) = G12(ii,jj);
1629 G22.assign(APTAPtmp);
1630
1631 // compute harmonic Ritz vectors
1632 Teuchos::SerialDenseMatrix<int,ScalarType> Ytmp( Teuchos::View, *Y_, (recycleBlocks_+numBlocks_), recycleBlocks_ );
1633 getHarmonicVecs(Ftmp,Gtmp,Ytmp);
1634
1635 // U1 = [U1 P(:,2:end-1)]*Y;
1636 index.resize( numBlocks_ );
1637 for (int ii=0; ii<numBlocks_; ++ii) { index[ii] = ii+1; }
1638 Teuchos::RCP<const MV> Ptmp = MVT::CloneView( *P_, index );
1639 index.resize( recycleBlocks_ );
1640 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1641 Teuchos::RCP<MV> PY2tmp = MVT::CloneViewNonConst( *PY2_, index );
1642 index.resize( recycleBlocks_ );
1643 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1644 Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1645 index.resize( recycleBlocks_ );
1646 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1647 Teuchos::RCP<MV> U1Y1tmp = MVT::CloneViewNonConst( *U1Y1_, index );
1648 MVT::MvTimesMatAddMv( one, *Ptmp, Y2, zero, *PY2tmp );
1649 MVT::MvTimesMatAddMv( one, *U1tmp, Y1, zero, *U1Y1tmp );
1650 MVT::MvAddMv(one,*U1Y1tmp, one, *PY2tmp, *U1tmp);
1651
1652 // Precompute some variables for next cycle
1653
1654 // AU1TAU1 = Y'*G*Y;
1655 Teuchos::SerialDenseMatrix<int,ScalarType> GYtmp( Teuchos::View, *GY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1656 GYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Gtmp,Ytmp,zero);
1657 AU1TAU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,GYtmp,zero);
1658
1659 // AU1TU1 = Y'*F*Y;
1660 Teuchos::SerialDenseMatrix<int,ScalarType> FYtmp( Teuchos::View, *FY_, (numBlocks_+recycleBlocks_), recycleBlocks_ );
1661 FYtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,Ftmp,Ytmp,zero);
1662 AU1TU1tmp.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,Ytmp,FYtmp,zero);
1663
1664 // dold = D(end);
1665 dold = Dtmp(numBlocks_-1,0);
1666
1667 // Indicate the size of the P, Beta structures generated this cycle
1668 lastp = numBlocks_+1;
1669 lastBeta = numBlocks_;
1670
1671 }
1672 }
1673 } // if (recycleBlocks_ > 0)
1674
1675 // Cleanup after end of cycle
1676
1677 // P = P(:,end-1:end);
1678 index.resize( 2 );
1679 index[0] = lastp-1; index[1] = lastp;
1680 Teuchos::RCP<const MV> Ptmp2 = MVT::CloneView( *P_, index );
1681 index[0] = 0; index[1] = 1;
1682 MVT::SetBlock(*Ptmp2,index,*P_);
1683
1684 // Beta = Beta(end);
1685 (*Beta_)(0,0) = (*Beta_)(lastBeta,0);
1686
1687 // Delta = Delta(:,end);
1688 if (existU_) { // Delta only initialized if U exists
1689 Teuchos::SerialDenseMatrix<int,ScalarType> mu1( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, 0 );
1690 Teuchos::SerialDenseMatrix<int,ScalarType> mu2( Teuchos::View, *Delta_, recycleBlocks_, 1, 0, numBlocks_ );
1691 mu1.assign(mu2);
1692 }
1693
1694 // Now reinitialize state variables for next cycle
1695 newstate.P = Teuchos::null;
1696 index.resize( numBlocks_+1 );
1697 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii+1; }
1698 newstate.P = MVT::CloneViewNonConst( *P_, index );
1699
1700 newstate.Beta = Teuchos::null;
1701 newstate.Beta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Beta_, numBlocks_, 1, 1, 0 ) );
1702
1703 newstate.Delta = Teuchos::null;
1704 newstate.Delta = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *Delta_, recycleBlocks_, numBlocks_, 0, 1 ) );
1705
1706 newstate.curDim = 1; // We have initialized the first search vector
1707
1708 // Pass to iteration object
1709 rcg_iter->initialize(newstate);
1710
1711 // increment cycle count
1712 cycle = cycle + 1;
1713
1714 }
1716 //
1717 // we returned from iterate(), but none of our status tests Passed.
1718 // something is wrong, and it is probably our fault.
1719 //
1721 else {
1722 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,
1723 "Belos::RCGSolMgr::solve(): Invalid return from RCGIter::iterate().");
1724 }
1725 }
1726 catch (const StatusTestNaNError& e) {
1727 // A NaN was detected in the solver. Set the solution to zero and return unconverged.
1728 achievedTol_ = MT::one();
1729 Teuchos::RCP<MV> X = problem_->getLHS();
1730 MVT::MvInit( *X, SCT::zero() );
1731 printer_->stream(Warnings) << "Belos::RCGSolMgr::solve(): Warning! NaN has been detected!"
1732 << std::endl;
1733 return Unconverged;
1734 }
1735 catch (const std::exception &e) {
1736 printer_->stream(Errors) << "Error! Caught std::exception in RCGIter::iterate() at iteration "
1737 << rcg_iter->getNumIters() << std::endl
1738 << e.what() << std::endl;
1739 throw;
1740 }
1741 }
1742
1743 // Inform the linear problem that we are finished with this block linear system.
1744 problem_->setCurrLS();
1745
1746 // Update indices for the linear systems to be solved.
1747 numRHS2Solve -= 1;
1748 if ( numRHS2Solve > 0 ) {
1749 currIdx[0]++;
1750 // Set the next indices.
1751 problem_->setLSIndex( currIdx );
1752 }
1753 else {
1754 currIdx.resize( numRHS2Solve );
1755 }
1756
1757 // Update the recycle space for the next linear system
1758 if (existU1_) { // be sure updated recycle space was created
1759 // U = U1
1760 index.resize(recycleBlocks_);
1761 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1762 MVT::SetBlock(*U1_,index,*U_);
1763 // Set flag indicating recycle space is now defined
1764 existU_ = true;
1765 if (numRHS2Solve > 0) { // also update AU, UTAU, and AUTAU
1766 // Free pointers in newstate
1767 newstate.P = Teuchos::null;
1768 newstate.Ap = Teuchos::null;
1769 newstate.r = Teuchos::null;
1770 newstate.z = Teuchos::null;
1771 newstate.U = Teuchos::null;
1772 newstate.AU = Teuchos::null;
1773 newstate.Alpha = Teuchos::null;
1774 newstate.Beta = Teuchos::null;
1775 newstate.D = Teuchos::null;
1776 newstate.Delta = Teuchos::null;
1777 newstate.LUUTAU = Teuchos::null;
1778 newstate.ipiv = Teuchos::null;
1779 newstate.rTz_old = Teuchos::null;
1780
1781 // Reinitialize AU, UTAU, AUTAU
1782 index.resize(recycleBlocks_);
1783 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1784 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1785 index.resize(recycleBlocks_);
1786 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1787 Teuchos::RCP<MV> AUtmp = MVT::CloneViewNonConst( *AU_, index );
1788 // Initialize AU
1789 problem_->applyOp( *Utmp, *AUtmp );
1790 // Initialize UTAU
1791 Teuchos::SerialDenseMatrix<int,ScalarType> UTAUtmp( Teuchos::View, *UTAU_, recycleBlocks_, recycleBlocks_ );
1792 MVT::MvTransMv( one, *Utmp, *AUtmp, UTAUtmp );
1793 // Initialize AUTAU ( AUTAU = AU'*(M\AU) )
1794 Teuchos::SerialDenseMatrix<int,ScalarType> AUTAUtmp( Teuchos::View, *AUTAU_, recycleBlocks_, recycleBlocks_ );
1795 if ( precObj != Teuchos::null ) {
1796 index.resize(recycleBlocks_);
1797 for (int i=0; i<recycleBlocks_; ++i) { index[i] = i; }
1798 index.resize(recycleBlocks_);
1799 for (int ii=0; ii<recycleBlocks_; ++ii) { index[ii] = ii; }
1800 Teuchos::RCP<MV> LeftPCAU = MVT::CloneViewNonConst( *U1_, index ); // use U1 as temp storage
1801 OPT::Apply( *precObj, *AUtmp, *LeftPCAU );
1802 MVT::MvTransMv( one, *AUtmp, *LeftPCAU, AUTAUtmp );
1803 } else {
1804 MVT::MvTransMv( one, *AUtmp, *AUtmp, AUTAUtmp );
1805 }
1806 } // if (numRHS2Solve > 0)
1807
1808 } // if (existU1)
1809 }// while ( numRHS2Solve > 0 )
1810
1811 }
1812
1813 // print final summary
1814 sTest_->print( printer_->stream(FinalSummary) );
1815
1816 // print timing information
1817#ifdef BELOS_TEUCHOS_TIME_MONITOR
1818 // Calling summarize() can be expensive, so don't call unless the
1819 // user wants to print out timing details. summarize() will do all
1820 // the work even if it's passed a "black hole" output stream.
1821 if (verbosity_ & TimingDetails)
1822 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1823#endif
1824
1825 // get iteration information for this solve
1826 numIters_ = maxIterTest_->getNumIters();
1827
1828 // Save the convergence test value ("achieved tolerance") for this solve.
1829 {
1830 using Teuchos::rcp_dynamic_cast;
1831 typedef StatusTestGenResNorm<ScalarType,MV,OP> conv_test_type;
1832 // testValues is nonnull and not persistent.
1833 const std::vector<MagnitudeType>* pTestValues =
1834 rcp_dynamic_cast<conv_test_type>(convTest_)->getTestValue();
1835
1836 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1837 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1838 "method returned NULL. Please report this bug to the Belos developers.");
1839
1840 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1841 "Belos::RCGSolMgr::solve(): The convergence test's getTestValue() "
1842 "method returned a vector of length zero. Please report this bug to the "
1843 "Belos developers.");
1844
1845 // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1846 // achieved tolerances for all vectors in the current solve(), or
1847 // just for the vectors from the last deflation?
1848 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1849 }
1850
1851 if (!isConverged) {
1852 return Unconverged; // return from RCGSolMgr::solve()
1853 }
1854 return Converged; // return from RCGSolMgr::solve()
1855}
1856
1857// Compute the harmonic eigenpairs of the projected, dense system.
1858template<class ScalarType, class MV, class OP>
1859void RCGSolMgr<ScalarType,MV,OP,true>::getHarmonicVecs(const Teuchos::SerialDenseMatrix<int,ScalarType>& F,
1860 const Teuchos::SerialDenseMatrix<int,ScalarType>& /* G */,
1861 Teuchos::SerialDenseMatrix<int,ScalarType>& Y) {
1862 // order of F,G
1863 int n = F.numCols();
1864
1865 // The LAPACK interface
1866 Teuchos::LAPACK<int,ScalarType> lapack;
1867
1868 // Magnitude of harmonic Ritz values
1869 std::vector<MagnitudeType> w(n);
1870
1871 // Sorted order of harmonic Ritz values
1872 std::vector<int> iperm(n);
1873
1874 // Compute k smallest harmonic Ritz pairs
1875 // SUBROUTINE DSYGV( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, LWORK, INFO )
1876 int itype = 1; // solve A*x = (lambda)*B*x
1877 char jobz='V'; // compute eigenvalues and eigenvectors
1878 char uplo='U'; // since F,G symmetric, reference only their upper triangular data
1879 std::vector<ScalarType> work(1);
1880 int lwork = -1;
1881 int info = 0;
1882 // since SYGV destroys workspace, create copies of F,G
1883 Teuchos::SerialDenseMatrix<int,ScalarType> F2( Teuchos::Copy, *F_ );
1884 Teuchos::SerialDenseMatrix<int,ScalarType> G2( Teuchos::Copy, *G_ );
1885
1886 // query for optimal workspace size
1887 lapack.SYGV(itype, jobz, uplo, n, G2.values(), G2.stride(), F2.values(), F2.stride(), &w[0], &work[0], lwork, &info);
1888 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, RCGSolMgrLAPACKFailure,
1889 "Belos::RCGSolMgr::solve(): LAPACK SYGV failed to query optimal work size.");
1890 lwork = (int)work[0];
1891 work.resize(lwork);
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 compute eigensolutions.");
1895
1896
1897 // Construct magnitude of each harmonic Ritz value
1898 this->sort(w,n,iperm);
1899
1900 // Select recycledBlocks_ smallest eigenvectors
1901 for( int i=0; i<recycleBlocks_; i++ ) {
1902 for( int j=0; j<n; j++ ) {
1903 Y(j,i) = G2(j,iperm[i]);
1904 }
1905 }
1906
1907}
1908
1909// This method sorts list of n floating-point numbers and return permutation vector
1910template<class ScalarType, class MV, class OP>
1911void RCGSolMgr<ScalarType,MV,OP,true>::sort(std::vector<ScalarType>& dlist, int n, std::vector<int>& iperm)
1912{
1913 int l, r, j, i, flag;
1914 int RR2;
1915 double dRR, dK;
1916
1917 // Initialize the permutation vector.
1918 for(j=0;j<n;j++)
1919 iperm[j] = j;
1920
1921 if (n <= 1) return;
1922
1923 l = n / 2 + 1;
1924 r = n - 1;
1925 l = l - 1;
1926 dRR = dlist[l - 1];
1927 dK = dlist[l - 1];
1928
1929 RR2 = iperm[l - 1];
1930 while (r != 0) {
1931 j = l;
1932 flag = 1;
1933
1934 while (flag == 1) {
1935 i = j;
1936 j = j + j;
1937
1938 if (j > r + 1)
1939 flag = 0;
1940 else {
1941 if (j < r + 1)
1942 if (dlist[j] > dlist[j - 1]) j = j + 1;
1943
1944 if (dlist[j - 1] > dK) {
1945 dlist[i - 1] = dlist[j - 1];
1946 iperm[i - 1] = iperm[j - 1];
1947 }
1948 else {
1949 flag = 0;
1950 }
1951 }
1952 }
1953 dlist[i - 1] = dRR;
1954 iperm[i - 1] = RR2;
1955 if (l == 1) {
1956 dRR = dlist [r];
1957 RR2 = iperm[r];
1958 dK = dlist[r];
1959 dlist[r] = dlist[0];
1960 iperm[r] = iperm[0];
1961 r = r - 1;
1962 }
1963 else {
1964 l = l - 1;
1965 dRR = dlist[l - 1];
1966 RR2 = iperm[l - 1];
1967 dK = dlist[l - 1];
1968 }
1969 }
1970 dlist[0] = dRR;
1971 iperm[0] = RR2;
1972}
1973
1974// This method requires the solver manager to return a std::string that describes itself.
1975template<class ScalarType, class MV, class OP>
1976std::string RCGSolMgr<ScalarType,MV,OP,true>::description() const
1977{
1978 std::ostringstream oss;
1979 oss << "Belos::RCGSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
1980 return oss.str();
1981}
1982
1983} // end Belos namespace
1984
1985#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:262
Belos::Details::SolverManagerRequiresRealLapack
Base class for Belos::SolverManager subclasses which normally can only compile with real ScalarType t...
Definition BelosSolverManager.hpp:409
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::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::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::StatusTestDetails
@ StatusTestDetails
Definition BelosTypes.hpp:229
Belos::FinalSummary
@ FinalSummary
Definition BelosTypes.hpp:227
Belos::TimingDetails
@ TimingDetails
Definition BelosTypes.hpp:228
Belos::ReturnType
ReturnType
Whether the Belos solve converged for all linear systems.
Definition BelosTypes.hpp:123
Belos::Unconverged
@ Unconverged
Definition BelosTypes.hpp:125
Belos::ResetType
ResetType
How to reset the solver.
Definition BelosTypes.hpp:174
Belos::RecycleSubspace
@ RecycleSubspace
Definition BelosTypes.hpp:175
Belos::DefaultSolverParameters::convTol
static const double convTol
Default convergence tolerance.
Definition BelosTypes.hpp:261
Generated for Belos by doxygen 1.9.8