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