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