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