Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosGCRODRSolMgr.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_GCRODR_SOLMGR_HPP
43#define BELOS_GCRODR_SOLMGR_HPP
44
48
49#include "BelosConfigDefs.hpp"
53#include "BelosTypes.hpp"
54
55#include "BelosGCRODRIter.hpp"
62#include "Teuchos_BLAS.hpp" // includes Teuchos_ConfigDefs.hpp
63#include "Teuchos_LAPACK.hpp"
64#include "Teuchos_as.hpp"
65#ifdef BELOS_TEUCHOS_TIME_MONITOR
66# include "Teuchos_TimeMonitor.hpp"
67#endif // BELOS_TEUCHOS_TIME_MONITOR
68#if defined(HAVE_TEUCHOSCORE_CXX11)
69# include <type_traits>
70#endif // defined(HAVE_TEUCHOSCORE_CXX11)
71
79namespace Belos {
80
82
83
91 public:
92 GCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
93 };
94
102 public:
103 GCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
104 };
105
113 public:
114 GCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
115 };
116
124 public:
125 GCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
126 };
127
129
154 template<class ScalarType, class MV, class OP,
155 const bool lapackSupportsScalarType =
158 public Details::SolverManagerRequiresLapack<ScalarType,MV,OP>
159 {
160 static const bool requiresLapack =
162 typedef Details::SolverManagerRequiresLapack<ScalarType, MV, OP,
164
165 public:
167 base_type ()
168 {}
169 GCRODRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem,
170 const Teuchos::RCP<Teuchos::ParameterList>& pl) :
171 base_type ()
172 {}
173 virtual ~GCRODRSolMgr () {}
174 };
175
179 template<class ScalarType, class MV, class OP>
180 class GCRODRSolMgr<ScalarType, MV, OP, true> :
181 public Details::SolverManagerRequiresLapack<ScalarType, MV, OP, true>
182 {
183
184#if defined(HAVE_TEUCHOSCORE_CXX11)
185# if defined(HAVE_TEUCHOS_COMPLEX)
186 #if defined(HAVE_TEUCHOS_LONG_DOUBLE)
187 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
188 std::is_same<ScalarType, std::complex<double> >::value ||
189 std::is_same<ScalarType, float>::value ||
190 std::is_same<ScalarType, double>::value ||
191 std::is_same<ScalarType, long double>::value,
192 "Belos::GCRODRSolMgr: ScalarType must be one of the four "
193 "types (S,D,C,Z) supported by LAPACK or long double (largely not impl'd).");
194 #else
195 static_assert (std::is_same<ScalarType, std::complex<float> >::value ||
196 std::is_same<ScalarType, std::complex<double> >::value ||
197 std::is_same<ScalarType, float>::value ||
198 std::is_same<ScalarType, double>::value,
199 "Belos::GCRODRSolMgr: ScalarType must be one of the four "
200 "types (S,D,C,Z) supported by LAPACK.");
201 #endif
202# else
203 #if defined(HAVE_TEUCHOS_LONG_DOUBLE)
204 static_assert (std::is_same<ScalarType, float>::value ||
205 std::is_same<ScalarType, double>::value ||
206 std::is_same<ScalarType, long double>::value,
207 "Belos::GCRODRSolMgr: ScalarType must be float, double or long double. "
208 "Complex arithmetic support is currently disabled. To "
209 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
210 #else
211 static_assert (std::is_same<ScalarType, float>::value ||
212 std::is_same<ScalarType, double>::value,
213 "Belos::GCRODRSolMgr: ScalarType must be float or double. "
214 "Complex arithmetic support is currently disabled. To "
215 "enable it, set Teuchos_ENABLE_COMPLEX=ON.");
216 #endif
217# endif // defined(HAVE_TEUCHOS_COMPLEX)
218#endif // defined(HAVE_TEUCHOSCORE_CXX11)
219
220 private:
223 typedef Teuchos::ScalarTraits<ScalarType> SCT;
224 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
225 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
227
228 public:
230
231
237 GCRODRSolMgr();
238
291 GCRODRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
292 const Teuchos::RCP<Teuchos::ParameterList> &pl);
293
295 virtual ~GCRODRSolMgr() {};
296
298 Teuchos::RCP<SolverManager<ScalarType, MV, OP> > clone () const override {
299 return Teuchos::rcp(new GCRODRSolMgr<ScalarType,MV,OP,true>);
300 }
302
304
305
309 return *problem_;
310 }
311
314 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const override;
315
318 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const override {
319 return params_;
320 }
321
327 Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
328 return Teuchos::tuple(timerSolve_);
329 }
330
336 MagnitudeType achievedTol() const override {
337 return achievedTol_;
338 }
339
341 int getNumIters() const override {
342 return numIters_;
343 }
344
347 bool isLOADetected() const override { return false; }
348
350
352
353
355 void setProblem( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem ) override {
356 problem_ = problem;
357 }
358
360 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params ) override;
361
363
365
366
370 void reset( const ResetType type ) override {
371 if ((type & Belos::Problem) && !Teuchos::is_null(problem_)) {
372 bool set = problem_->setProblem();
373 if (!set)
374 throw "Could not set problem.";
375 }
376 else if (type & Belos::RecycleSubspace) {
377 keff = 0;
378 }
379 }
381
383
384
411 ReturnType solve() override;
412
414
416
418 std::string description() const override;
419
421
422 private:
423
424 // Called by all constructors; Contains init instructions common to all constructors
425 void init();
426
427 // Initialize solver state storage
428 void initializeStateStorage();
429
430 // Compute updated recycle space given existing recycle space and newly generated Krylov space
431 void buildRecycleSpace2(Teuchos::RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter);
432
433 // Computes harmonic eigenpairs of projected matrix created during the priming solve.
434 // HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension m+1 x m.
435 // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
436 // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
437 int getHarmonicVecs1(int m,
438 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
439 Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
440
441 // Computes harmonic eigenpairs of projected matrix created during one cycle.
442 // HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+m+1 x keff+m.
443 // VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) m+1 vectors.
444 // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
445 // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
446 int getHarmonicVecs2(int keff, int m,
447 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
448 const Teuchos::RCP<const MV>& VV,
449 Teuchos::SerialDenseMatrix<int,ScalarType>& PP);
450
451 // Sort list of n floating-point numbers and return permutation vector
452 void sort(std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm);
453
454 // Lapack interface
455 Teuchos::LAPACK<int,ScalarType> lapack;
456
457 // Linear problem.
458 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
459
460 // Output manager.
461 Teuchos::RCP<OutputManager<ScalarType> > printer_;
462 Teuchos::RCP<std::ostream> outputStream_;
463
464 // Status test.
465 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
466 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
467 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
468 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
469 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
470
474 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
475
476 // Current parameter list.
477 Teuchos::RCP<Teuchos::ParameterList> params_;
478
479 // Default solver values.
480 static constexpr double orthoKappa_default_ = 0.0;
481 static constexpr int maxRestarts_default_ = 100;
482 static constexpr int maxIters_default_ = 1000;
483 static constexpr int numBlocks_default_ = 50;
484 static constexpr int blockSize_default_ = 1;
485 static constexpr int recycledBlocks_default_ = 5;
486 static constexpr int verbosity_default_ = Belos::Errors;
487 static constexpr int outputStyle_default_ = Belos::General;
488 static constexpr int outputFreq_default_ = -1;
489 static constexpr const char * impResScale_default_ = "Norm of Preconditioned Initial Residual";
490 static constexpr const char * expResScale_default_ = "Norm of Initial Residual";
491 static constexpr const char * label_default_ = "Belos";
492 static constexpr const char * orthoType_default_ = "ICGS";
493
494 // Current solver values.
495 MagnitudeType convTol_, orthoKappa_, achievedTol_;
496 int maxRestarts_, maxIters_, numIters_;
497 int verbosity_, outputStyle_, outputFreq_;
498 std::string orthoType_;
499 std::string impResScale_, expResScale_;
500
502 // Solver State Storage
504 //
505 // The number of blocks and recycle blocks (m and k, respectively)
506 int numBlocks_, recycledBlocks_;
507 // Current size of recycled subspace
508 int keff;
509 //
510 // Residual vector
511 Teuchos::RCP<MV> r_;
512 //
513 // Search space
514 Teuchos::RCP<MV> V_;
515 //
516 // Recycled subspace and its image
517 Teuchos::RCP<MV> U_, C_;
518 //
519 // Updated recycle space and its image
520 Teuchos::RCP<MV> U1_, C1_;
521 //
522 // Storage used in constructing harmonic Ritz values/vectors
523 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2_;
524 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H_;
525 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B_;
526 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PP_;
527 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > HP_;
528 std::vector<ScalarType> tau_;
529 std::vector<ScalarType> work_;
530 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > R_;
531 std::vector<int> ipiv_;
533
534 // Timers.
535 std::string label_;
536 Teuchos::RCP<Teuchos::Time> timerSolve_;
537
538 // Internal state variables.
539 bool isSet_;
540
541 // Have we generated or regenerated a recycle space yet this solve?
543 };
544
545
546// Empty Constructor
547template<class ScalarType, class MV, class OP>
549 achievedTol_(0.0),
550 numIters_(0)
551{
552 init ();
553}
554
555
556// Basic Constructor
557template<class ScalarType, class MV, class OP>
559GCRODRSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem,
560 const Teuchos::RCP<Teuchos::ParameterList>& pl):
561 achievedTol_(0.0),
562 numIters_(0)
563{
564 // Initialize local pointers to null, and initialize local variables
565 // to default values.
566 init ();
567
568 TEUCHOS_TEST_FOR_EXCEPTION(
569 problem == Teuchos::null, std::invalid_argument,
570 "Belos::GCRODRSolMgr constructor: The solver manager's "
571 "constructor needs the linear problem argument 'problem' "
572 "to be non-null.");
573 problem_ = problem;
574
575 // Set the parameters using the list that was passed in. If null,
576 // we defer initialization until a non-null list is set (by the
577 // client calling setParameters(), or by calling solve() -- in
578 // either case, a null parameter list indicates that default
579 // parameters should be used).
580 if (! pl.is_null ()) {
581 setParameters (pl);
582 }
583}
584
585// Common instructions executed in all constructors
586template<class ScalarType, class MV, class OP>
588 outputStream_ = Teuchos::rcpFromRef(std::cout);
590 orthoKappa_ = orthoKappa_default_;
591 maxRestarts_ = maxRestarts_default_;
592 maxIters_ = maxIters_default_;
593 numBlocks_ = numBlocks_default_;
594 recycledBlocks_ = recycledBlocks_default_;
595 verbosity_ = verbosity_default_;
596 outputStyle_ = outputStyle_default_;
597 outputFreq_ = outputFreq_default_;
598 orthoType_ = orthoType_default_;
599 impResScale_ = impResScale_default_;
600 expResScale_ = expResScale_default_;
601 label_ = label_default_;
602 isSet_ = false;
603 builtRecycleSpace_ = false;
604 keff = 0;
605 r_ = Teuchos::null;
606 V_ = Teuchos::null;
607 U_ = Teuchos::null;
608 C_ = Teuchos::null;
609 U1_ = Teuchos::null;
610 C1_ = Teuchos::null;
611 PP_ = Teuchos::null;
612 HP_ = Teuchos::null;
613 H2_ = Teuchos::null;
614 R_ = Teuchos::null;
615 H_ = Teuchos::null;
616 B_ = Teuchos::null;
617}
618
619template<class ScalarType, class MV, class OP>
620void
622setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params)
623{
624 using Teuchos::isParameterType;
625 using Teuchos::getParameter;
626 using Teuchos::null;
627 using Teuchos::ParameterList;
628 using Teuchos::parameterList;
629 using Teuchos::RCP;
630 using Teuchos::rcp;
631 using Teuchos::rcp_dynamic_cast;
632 using Teuchos::rcpFromRef;
633 using Teuchos::Exceptions::InvalidParameter;
634 using Teuchos::Exceptions::InvalidParameterName;
635 using Teuchos::Exceptions::InvalidParameterType;
636
637 // The default parameter list contains all parameters that
638 // GCRODRSolMgr understands, and none that it doesn't understand.
639 RCP<const ParameterList> defaultParams = getValidParameters();
640
641 // Create the internal parameter list if one doesn't already exist.
642 //
643 // (mfh 28 Feb 2011, 10 Mar 2011) At the time this code was written,
644 // ParameterList did not have validators or validateParameters().
645 // This is why the code below carefully validates the parameters one
646 // by one and fills in defaults. This code could be made a lot
647 // shorter by using validators. To do so, we would have to define
648 // appropriate validators for all the parameters. (This would more
649 // or less just move all that validation code out of this routine
650 // into to getValidParameters().)
651 //
652 // For an analogous reason, GCRODRSolMgr defines default parameter
653 // values as class data, as well as in the default ParameterList.
654 // This redundancy could be removed by defining the default
655 // parameter values only in the default ParameterList (which
656 // documents each parameter as well -- handy!).
657 if (params_.is_null()) {
658 params_ = parameterList (*defaultParams);
659 } else {
660 // A common case for setParameters() is for it to be called at the
661 // beginning of the solve() routine. This follows the Belos
662 // pattern of delaying initialization until the last possible
663 // moment (when the user asks Belos to perform the solve). In
664 // this common case, we save ourselves a deep copy of the input
665 // parameter list.
666 if (params_ != params) {
667 // Make a deep copy of the input parameter list. This allows
668 // the caller to modify or change params later, without
669 // affecting the behavior of this solver. This solver will then
670 // only change its internal parameters if setParameters() is
671 // called again.
672 params_ = parameterList (*params);
673 }
674
675 // Fill in any missing parameters and their default values. Also,
676 // throw an exception if the parameter list has any misspelled or
677 // "extra" parameters. If you don't like this behavior, you'll
678 // want to replace the line of code below with your desired
679 // validation scheme. Note that Teuchos currently only implements
680 // two options:
681 //
682 // 1. validateParameters() requires that params_ has all the
683 // parameters that the default list has, and none that it
684 // doesn't have.
685 //
686 // 2. validateParametersAndSetDefaults() fills in missing
687 // parameters in params_ using the default list, but requires
688 // that any parameter provided in params_ is also in the
689 // default list.
690 //
691 // Here is an easy way to ignore any "extra" or misspelled
692 // parameters: Make a deep copy of the default list, fill in any
693 // "missing" parameters from the _input_ list, and then validate
694 // the input list using the deep copy of the default list. We
695 // show this in code:
696 //
697 // RCP<ParameterList> defaultCopy = parameterList (*getValidParameters ());
698 // defaultCopy->validateParametersAndSetDefaults (params);
699 // params->validateParametersAndSetDefaults (defaultCopy);
700 //
701 // This method is not entirely robust, because the input list may
702 // have incorrect validators set for existing parameters in the
703 // default list. This would then cause "validation" of the
704 // default list to throw an exception. As a result, we've chosen
705 // for now to be intolerant of misspellings and "extra" parameters
706 // in the input list.
707 params_->validateParametersAndSetDefaults (*defaultParams);
708 }
709
710 // Check for maximum number of restarts.
711 if (params->isParameter ("Maximum Restarts")) {
712 maxRestarts_ = params->get("Maximum Restarts", maxRestarts_default_);
713
714 // Update parameter in our list.
715 params_->set ("Maximum Restarts", maxRestarts_);
716 }
717
718 // Check for maximum number of iterations
719 if (params->isParameter ("Maximum Iterations")) {
720 maxIters_ = params->get ("Maximum Iterations", maxIters_default_);
721
722 // Update parameter in our list and in status test.
723 params_->set ("Maximum Iterations", maxIters_);
724 if (! maxIterTest_.is_null())
725 maxIterTest_->setMaxIters (maxIters_);
726 }
727
728 // Check for the maximum number of blocks.
729 if (params->isParameter ("Num Blocks")) {
730 numBlocks_ = params->get ("Num Blocks", numBlocks_default_);
731 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks_ <= 0, std::invalid_argument,
732 "Belos::GCRODRSolMgr: The \"Num Blocks\" parameter must "
733 "be strictly positive, but you specified a value of "
734 << numBlocks_ << ".");
735 // Update parameter in our list.
736 params_->set ("Num Blocks", numBlocks_);
737 }
738
739 // Check for the maximum number of blocks.
740 if (params->isParameter ("Num Recycled Blocks")) {
741 recycledBlocks_ = params->get ("Num Recycled Blocks",
742 recycledBlocks_default_);
743 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ <= 0, std::invalid_argument,
744 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
745 "parameter must be strictly positive, but you specified "
746 "a value of " << recycledBlocks_ << ".");
747 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks_ >= numBlocks_, std::invalid_argument,
748 "Belos::GCRODRSolMgr: The \"Num Recycled Blocks\" "
749 "parameter must be less than the \"Num Blocks\" "
750 "parameter, but you specified \"Num Recycled Blocks\" "
751 "= " << recycledBlocks_ << " and \"Num Blocks\" = "
752 << numBlocks_ << ".");
753 // Update parameter in our list.
754 params_->set("Num Recycled Blocks", recycledBlocks_);
755 }
756
757 // Check to see if the timer label changed. If it did, update it in
758 // the parameter list, and create a new timer with that label (if
759 // Belos was compiled with timers enabled).
760 if (params->isParameter ("Timer Label")) {
761 std::string tempLabel = params->get ("Timer Label", label_default_);
762
763 // Update parameter in our list and solver timer
764 if (tempLabel != label_) {
765 label_ = tempLabel;
766 params_->set ("Timer Label", label_);
767 std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
768#ifdef BELOS_TEUCHOS_TIME_MONITOR
769 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
770#endif
771 if (ortho_ != Teuchos::null) {
772 ortho_->setLabel( label_ );
773 }
774 }
775 }
776
777 // Check for a change in verbosity level
778 if (params->isParameter ("Verbosity")) {
779 if (isParameterType<int> (*params, "Verbosity")) {
780 verbosity_ = params->get ("Verbosity", verbosity_default_);
781 } else {
782 verbosity_ = (int) getParameter<Belos::MsgType> (*params, "Verbosity");
783 }
784 // Update parameter in our list.
785 params_->set ("Verbosity", verbosity_);
786 // If the output manager (printer_) is null, then we will
787 // instantiate it later with the correct verbosity.
788 if (! printer_.is_null())
789 printer_->setVerbosity (verbosity_);
790 }
791
792 // Check for a change in output style
793 if (params->isParameter ("Output Style")) {
794 if (isParameterType<int> (*params, "Output Style")) {
795 outputStyle_ = params->get ("Output Style", outputStyle_default_);
796 } else {
797 outputStyle_ = (int) getParameter<OutputType> (*params, "Output Style");
798 }
799
800 // Update parameter in our list.
801 params_->set ("Output Style", outputStyle_);
802 // We will (re)instantiate the output status test afresh below.
803 outputTest_ = null;
804 }
805
806 // Get the output stream for the output manager.
807 //
808 // While storing the output stream in the parameter list (either as
809 // an RCP or as a nonconst reference) is convenient and safe for
810 // programming, it makes it impossible to serialize the parameter
811 // list, read it back in from the serialized representation, and get
812 // the same output stream as before. This is because output streams
813 // may be arbitrary constructed objects.
814 //
815 // In case the user tried reading in the parameter list from a
816 // serialized representation and the output stream can't be read
817 // back in, we set the output stream to point to std::cout. This
818 // ensures reasonable behavior.
819 if (params->isParameter ("Output Stream")) {
820 try {
821 outputStream_ = getParameter<RCP<std::ostream> > (*params, "Output Stream");
822 } catch (InvalidParameter&) {
823 outputStream_ = rcpFromRef (std::cout);
824 }
825 // We assume that a null output stream indicates that the user
826 // doesn't want to print anything, so we replace it with a "black
827 // hole" stream that prints nothing sent to it. (We can't use a
828 // null output stream, since the output manager always sends
829 // things it wants to print to the output stream.)
830 if (outputStream_.is_null()) {
831 outputStream_ = rcp (new Teuchos::oblackholestream);
832 }
833 // Update parameter in our list.
834 params_->set ("Output Stream", outputStream_);
835 // If the output manager (printer_) is null, then we will
836 // instantiate it later with the correct output stream.
837 if (! printer_.is_null()) {
838 printer_->setOStream (outputStream_);
839 }
840 }
841
842 // frequency level
843 if (verbosity_ & Belos::StatusTestDetails) {
844 if (params->isParameter ("Output Frequency")) {
845 outputFreq_ = params->get ("Output Frequency", outputFreq_default_);
846 }
847
848 // Update parameter in out list and output status test.
849 params_->set("Output Frequency", outputFreq_);
850 if (! outputTest_.is_null())
851 outputTest_->setOutputFrequency (outputFreq_);
852 }
853
854 // Create output manager if we need to, using the verbosity level
855 // and output stream that we fetched above. We do this here because
856 // instantiating an OrthoManager using OrthoManagerFactory requires
857 // a valid OutputManager.
858 if (printer_.is_null()) {
859 printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
860 }
861
862 // Get the orthogonalization manager name ("Orthogonalization").
863 //
864 // Getting default values for the orthogonalization manager
865 // parameters ("Orthogonalization Parameters") requires knowing the
866 // orthogonalization manager name. Save it for later, and also
867 // record whether it's different than before.
869 bool changedOrthoType = false;
870 if (params->isParameter ("Orthogonalization")) {
871 const std::string& tempOrthoType =
872 params->get ("Orthogonalization", orthoType_default_);
873 // Ensure that the specified orthogonalization type is valid.
874 if (! factory.isValidName (tempOrthoType)) {
875 std::ostringstream os;
876 os << "Belos::GCRODRSolMgr: Invalid orthogonalization name \""
877 << tempOrthoType << "\". The following are valid options "
878 << "for the \"Orthogonalization\" name parameter: ";
879 factory.printValidNames (os);
880 throw std::invalid_argument (os.str());
881 }
882 if (tempOrthoType != orthoType_) {
883 changedOrthoType = true;
884 orthoType_ = tempOrthoType;
885 // Update parameter in our list.
886 params_->set ("Orthogonalization", orthoType_);
887 }
888 }
889
890 // Get any parameters for the orthogonalization ("Orthogonalization
891 // Parameters"). If not supplied, the orthogonalization manager
892 // factory will supply default values.
893 //
894 // NOTE (mfh 12 Jan 2011) For the sake of backwards compatibility,
895 // if params has an "Orthogonalization Constant" parameter and the
896 // DGKS orthogonalization manager is to be used, the value of this
897 // parameter will override DGKS's "depTol" parameter.
898 //
899 // Users must supply the orthogonalization manager parameters as a
900 // sublist (supplying it as an RCP<ParameterList> would make the
901 // resulting parameter list not serializable).
902 RCP<ParameterList> orthoParams;
903 { // The nonmember function sublist() returns an RCP<ParameterList>,
904 // which is what we want here.
905 using Teuchos::sublist;
906 // Abbreviation to avoid typos.
907 const std::string paramName ("Orthogonalization Parameters");
908
909 try {
910 orthoParams = sublist (params_, paramName, true);
911 } catch (InvalidParameter&) {
912 // We didn't get the parameter list from params, so get a
913 // default parameter list from the OrthoManagerFactory. Modify
914 // params_ so that it has the default parameter list, and set
915 // orthoParams to ensure it's a sublist of params_ (and not just
916 // a copy of one).
917 params_->set (paramName, factory.getDefaultParameters (orthoType_));
918 orthoParams = sublist (params_, paramName, true);
919 }
920 }
921 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
922 "Failed to get orthogonalization parameters. "
923 "Please report this bug to the Belos developers.");
924
925 // Instantiate a new MatOrthoManager subclass instance if necessary.
926 // If not necessary, then tell the existing instance about the new
927 // parameters.
928 if (ortho_.is_null() || changedOrthoType) {
929 // We definitely need to make a new MatOrthoManager, since either
930 // we haven't made one yet, or we've changed orthogonalization
931 // methods. Creating the orthogonalization manager requires that
932 // the OutputManager (printer_) already be initialized.
933 ortho_ = factory.makeMatOrthoManager (orthoType_, null, printer_,
934 label_, orthoParams);
935 } else {
936 // If the MatOrthoManager implements the ParameterListAcceptor
937 // mix-in interface, we can propagate changes to its parameters
938 // without reinstantiating the MatOrthoManager.
939 //
940 // We recommend that all MatOrthoManager subclasses implement
941 // Teuchos::ParameterListAcceptor, but do not (yet) require this.
942 typedef Teuchos::ParameterListAcceptor PLA;
943 RCP<PLA> pla = rcp_dynamic_cast<PLA> (ortho_);
944 if (pla.is_null()) {
945 // Oops, it's not a ParameterListAcceptor. We have to
946 // reinstantiate the MatOrthoManager in order to pass in the
947 // possibly new parameters.
948 ortho_ = factory.makeMatOrthoManager (orthoType_, null, printer_,
949 label_, orthoParams);
950 } else {
951 pla->setParameterList (orthoParams);
952 }
953 }
954
955 // The DGKS orthogonalization accepts a "Orthogonalization Constant"
956 // parameter (also called kappa in the code, but not in the
957 // parameter list). If its value is provided in the given parameter
958 // list, and its value is positive, use it. Ignore negative values.
959 //
960 // NOTE (mfh 12 Jan 2011) This overrides the "depTol" parameter that
961 // may have been specified in "Orthogonalization Parameters". We
962 // retain this behavior for backwards compatibility.
963 if (params->isParameter ("Orthogonalization Constant")) {
964 MagnitudeType orthoKappa = orthoKappa_default_;
965 if (params->isType<MagnitudeType> ("Orthogonalization Constant")) {
966 orthoKappa = params->get ("Orthogonalization Constant", orthoKappa);
967 }
968 else {
969 orthoKappa = params->get ("Orthogonalization Constant", orthoKappa_default_);
970 }
971
972 if (orthoKappa > 0) {
973 orthoKappa_ = orthoKappa;
974 // Update parameter in our list.
975 params_->set("Orthogonalization Constant", orthoKappa_);
976 // Only DGKS currently accepts this parameter.
977 if (orthoType_ == "DGKS" && ! ortho_.is_null()) {
978 typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
979 // This cast should always succeed; it's a bug
980 // otherwise. (If the cast fails, then orthoType_
981 // doesn't correspond to the OrthoManager subclass
982 // instance that we think we have, so we initialized the
983 // wrong subclass somehow.)
984 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
985 }
986 }
987 }
988
989 // Convergence
990 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
991 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
992
993 // Check for convergence tolerance
994 if (params->isParameter("Convergence Tolerance")) {
995 if (params->isType<MagnitudeType> ("Convergence Tolerance")) {
996 convTol_ = params->get ("Convergence Tolerance",
998 }
999 else {
1000 convTol_ = params->get ("Convergence Tolerance", DefaultSolverParameters::convTol);
1001 }
1002
1003 // Update parameter in our list and residual tests.
1004 params_->set ("Convergence Tolerance", convTol_);
1005 if (! impConvTest_.is_null())
1006 impConvTest_->setTolerance (convTol_);
1007 if (! expConvTest_.is_null())
1008 expConvTest_->setTolerance (convTol_);
1009 }
1010
1011 // Check for a change in scaling, if so we need to build new residual tests.
1012 if (params->isParameter ("Implicit Residual Scaling")) {
1013 std::string tempImpResScale =
1014 getParameter<std::string> (*params, "Implicit Residual Scaling");
1015
1016 // Only update the scaling if it's different.
1017 if (impResScale_ != tempImpResScale) {
1018 ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
1019 impResScale_ = tempImpResScale;
1020
1021 // Update parameter in our list and residual tests
1022 params_->set("Implicit Residual Scaling", impResScale_);
1023 // NOTE (mfh 28 Feb 2011) StatusTestImpResNorm only lets you
1024 // call defineScaleForm() once. The code below attempts to call
1025 // defineScaleForm(); if the scale form has already been
1026 // defined, it constructs a new StatusTestImpResNorm instance.
1027 // StatusTestImpResNorm really should not expose the
1028 // defineScaleForm() method, since it's serving an
1029 // initialization purpose; all initialization should happen in
1030 // the constructor whenever possible. In that case, the code
1031 // below could be simplified into a single (re)instantiation.
1032 if (! impConvTest_.is_null()) {
1033 try {
1034 impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
1035 }
1036 catch (StatusTestError&) {
1037 // Delete the convergence test so it gets constructed again.
1038 impConvTest_ = null;
1039 convTest_ = null;
1040 }
1041 }
1042 }
1043 }
1044
1045 if (params->isParameter("Explicit Residual Scaling")) {
1046 std::string tempExpResScale =
1047 getParameter<std::string> (*params, "Explicit Residual Scaling");
1048
1049 // Only update the scaling if it's different.
1050 if (expResScale_ != tempExpResScale) {
1051 ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
1052 expResScale_ = tempExpResScale;
1053
1054 // Update parameter in our list and residual tests
1055 params_->set("Explicit Residual Scaling", expResScale_);
1056 // NOTE (mfh 28 Feb 2011) See note above on the (re)construction
1057 // of StatusTestImpResNorm.
1058 if (! expConvTest_.is_null()) {
1059 try {
1060 expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
1061 }
1062 catch (StatusTestError&) {
1063 // Delete the convergence test so it gets constructed again.
1064 expConvTest_ = null;
1065 convTest_ = null;
1066 }
1067 }
1068 }
1069 }
1070 //
1071 // Create iteration stopping criteria ("status tests") if we need
1072 // to, by combining three different stopping criteria.
1073 //
1074 // First, construct maximum-number-of-iterations stopping criterion.
1075 if (maxIterTest_.is_null())
1076 maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
1077
1078 // Implicit residual test, using the native residual to determine if
1079 // convergence was achieved.
1080 if (impConvTest_.is_null()) {
1081 impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1082 impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_),
1084 }
1085
1086 // Explicit residual test once the native residual is below the tolerance
1087 if (expConvTest_.is_null()) {
1088 expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1089 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
1090 expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_),
1092 }
1093 // Convergence test first tests the implicit residual, then the
1094 // explicit residual if the implicit residual test passes.
1095 if (convTest_.is_null()) {
1096 convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1097 impConvTest_,
1098 expConvTest_));
1099 }
1100 // Construct the complete iteration stopping criterion:
1101 //
1102 // "Stop iterating if the maximum number of iterations has been
1103 // reached, or if the convergence test passes."
1104 sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR,
1105 maxIterTest_,
1106 convTest_));
1107 // Create the status test output class.
1108 // This class manages and formats the output from the status test.
1109 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
1110 outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
1112
1113 // Set the solver string for the output test
1114 std::string solverDesc = " GCRODR ";
1115 outputTest_->setSolverDesc( solverDesc );
1116
1117 // Create the timer if we need to.
1118 if (timerSolve_.is_null()) {
1119 std::string solveLabel = label_ + ": GCRODRSolMgr total solve time";
1120#ifdef BELOS_TEUCHOS_TIME_MONITOR
1121 timerSolve_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
1122#endif
1123 }
1124
1125 // Inform the solver manager that the current parameters were set.
1126 isSet_ = true;
1127}
1128
1129
1130template<class ScalarType, class MV, class OP>
1131Teuchos::RCP<const Teuchos::ParameterList>
1133{
1134 using Teuchos::ParameterList;
1135 using Teuchos::parameterList;
1136 using Teuchos::RCP;
1137
1138 static RCP<const ParameterList> validPL;
1139 if (is_null(validPL)) {
1140 RCP<ParameterList> pl = parameterList ();
1141
1142 // Set all the valid parameters and their default values.
1143 pl->set("Convergence Tolerance", static_cast<MagnitudeType>(DefaultSolverParameters::convTol),
1144 "The relative residual tolerance that needs to be achieved by the\n"
1145 "iterative solver in order for the linear system to be declared converged.");
1146 pl->set("Maximum Restarts", static_cast<int>(maxRestarts_default_),
1147 "The maximum number of cycles allowed for each\n"
1148 "set of RHS solved.");
1149 pl->set("Maximum Iterations", static_cast<int>(maxIters_default_),
1150 "The maximum number of iterations allowed for each\n"
1151 "set of RHS solved.");
1152 // mfh 25 Oct 2010: "Block Size" must be 1 because GCRODR is
1153 // currently not a block method: i.e., it does not work on
1154 // multiple right-hand sides at once.
1155 pl->set("Block Size", static_cast<int>(blockSize_default_),
1156 "Block Size Parameter -- currently must be 1 for GCRODR");
1157 pl->set("Num Blocks", static_cast<int>(numBlocks_default_),
1158 "The maximum number of vectors allowed in the Krylov subspace\n"
1159 "for each set of RHS solved.");
1160 pl->set("Num Recycled Blocks", static_cast<int>(recycledBlocks_default_),
1161 "The maximum number of vectors in the recycled subspace." );
1162 pl->set("Verbosity", static_cast<int>(verbosity_default_),
1163 "What type(s) of solver information should be outputted\n"
1164 "to the output stream.");
1165 pl->set("Output Style", static_cast<int>(outputStyle_default_),
1166 "What style is used for the solver information outputted\n"
1167 "to the output stream.");
1168 pl->set("Output Frequency", static_cast<int>(outputFreq_default_),
1169 "How often convergence information should be outputted\n"
1170 "to the output stream.");
1171 pl->set("Output Stream", Teuchos::rcpFromRef(std::cout),
1172 "A reference-counted pointer to the output stream where all\n"
1173 "solver output is sent.");
1174 pl->set("Implicit Residual Scaling", static_cast<const char *>(impResScale_default_),
1175 "The type of scaling used in the implicit residual convergence test.");
1176 pl->set("Explicit Residual Scaling", static_cast<const char *>(expResScale_default_),
1177 "The type of scaling used in the explicit residual convergence test.");
1178 pl->set("Timer Label", static_cast<const char *>(label_default_),
1179 "The string to use as a prefix for the timer labels.");
1180 {
1182 pl->set("Orthogonalization", static_cast<const char *>(orthoType_default_),
1183 "The type of orthogonalization to use. Valid options: " +
1184 factory.validNamesString());
1185 RCP<const ParameterList> orthoParams =
1186 factory.getDefaultParameters (orthoType_default_);
1187 pl->set ("Orthogonalization Parameters", *orthoParams,
1188 "Parameters specific to the type of orthogonalization used.");
1189 }
1190 pl->set("Orthogonalization Constant",static_cast<MagnitudeType>(orthoKappa_default_),
1191 "When using DGKS orthogonalization: the \"depTol\" constant, used "
1192 "to determine whether another step of classical Gram-Schmidt is "
1193 "necessary. Otherwise ignored.");
1194 validPL = pl;
1195 }
1196 return validPL;
1197}
1198
1199// initializeStateStorage
1200template<class ScalarType, class MV, class OP>
1202
1203 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1204
1205 // Check if there is any multivector to clone from.
1206 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1207 if (rhsMV == Teuchos::null) {
1208 // Nothing to do
1209 return;
1210 }
1211 else {
1212
1213 // Initialize the state storage
1214 TEUCHOS_TEST_FOR_EXCEPTION(static_cast<ptrdiff_t>(numBlocks_) > MVT::GetGlobalLength(*rhsMV),std::invalid_argument,
1215 "Belos::GCRODRSolMgr::initializeStateStorage(): Cannot generate a Krylov basis with dimension larger the operator!");
1216
1217 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1218 if (U_ == Teuchos::null) {
1219 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1220 }
1221 else {
1222 // Generate U_ by cloning itself ONLY if more space is needed.
1223 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1224 Teuchos::RCP<const MV> tmp = U_;
1225 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1226 }
1227 }
1228
1229 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1230 if (C_ == Teuchos::null) {
1231 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1232 }
1233 else {
1234 // Generate C_ by cloning itself ONLY if more space is needed.
1235 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1236 Teuchos::RCP<const MV> tmp = C_;
1237 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1238 }
1239 }
1240
1241 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1242 if (V_ == Teuchos::null) {
1243 V_ = MVT::Clone( *rhsMV, numBlocks_+1 );
1244 }
1245 else {
1246 // Generate V_ by cloning itself ONLY if more space is needed.
1247 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1248 Teuchos::RCP<const MV> tmp = V_;
1249 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1250 }
1251 }
1252
1253 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1254 if (U1_ == Teuchos::null) {
1255 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1256 }
1257 else {
1258 // Generate U1_ by cloning itself ONLY if more space is needed.
1259 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1260 Teuchos::RCP<const MV> tmp = U1_;
1261 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1262 }
1263 }
1264
1265 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1266 if (C1_ == Teuchos::null) {
1267 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1268 }
1269 else {
1270 // Generate C1_ by cloning itself ONLY if more space is needed.
1271 if (MVT::GetNumberVecs(*C1_) < recycledBlocks_+1) {
1272 Teuchos::RCP<const MV> tmp = C1_;
1273 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1274 }
1275 }
1276
1277 // Generate r_ only if it doesn't exist
1278 if (r_ == Teuchos::null)
1279 r_ = MVT::Clone( *rhsMV, 1 );
1280
1281 // Size of tau_ will change during computation, so just be sure it starts with appropriate size
1282 tau_.resize(recycledBlocks_+1);
1283
1284 // Size of work_ will change during computation, so just be sure it starts with appropriate size
1285 work_.resize(recycledBlocks_+1);
1286
1287 // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
1288 ipiv_.resize(recycledBlocks_+1);
1289
1290 // Generate H2_ only if it doesn't exist, otherwise resize it.
1291 if (H2_ == Teuchos::null)
1292 H2_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1293 else {
1294 if ( (H2_->numRows() != numBlocks_+recycledBlocks_+2) || (H2_->numCols() != numBlocks_+recycledBlocks_+1) )
1295 H2_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1296 }
1297 H2_->putScalar(zero);
1298
1299 // Generate R_ only if it doesn't exist, otherwise resize it.
1300 if (R_ == Teuchos::null)
1301 R_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( recycledBlocks_+1, recycledBlocks_+1 ) );
1302 else {
1303 if ( (R_->numRows() != recycledBlocks_+1) || (R_->numCols() != recycledBlocks_+1) )
1304 R_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1305 }
1306 R_->putScalar(zero);
1307
1308 // Generate PP_ only if it doesn't exist, otherwise resize it.
1309 if (PP_ == Teuchos::null)
1310 PP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 ) );
1311 else {
1312 if ( (PP_->numRows() != numBlocks_+recycledBlocks_+2) || (PP_->numCols() != recycledBlocks_+1) )
1313 PP_->reshape( numBlocks_+recycledBlocks_+2, recycledBlocks_+1 );
1314 }
1315
1316 // Generate HP_ only if it doesn't exist, otherwise resize it.
1317 if (HP_ == Teuchos::null)
1318 HP_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 ) );
1319 else {
1320 if ( (HP_->numRows() != numBlocks_+recycledBlocks_+2) || (HP_->numCols() != numBlocks_+recycledBlocks_+1) )
1321 HP_->reshape( numBlocks_+recycledBlocks_+2, numBlocks_+recycledBlocks_+1 );
1322 }
1323
1324 } // end else
1325}
1326
1327
1328// solve()
1329template<class ScalarType, class MV, class OP>
1331 using Teuchos::RCP;
1332 using Teuchos::rcp;
1333
1334 // Set the current parameters if they were not set before.
1335 // NOTE: This may occur if the user generated the solver manager with the default constructor and
1336 // then didn't set any parameters using setParameters().
1337 if (!isSet_) { setParameters( params_ ); }
1338
1339 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1340 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1341 std::vector<int> index(numBlocks_+1);
1342
1343 TEUCHOS_TEST_FOR_EXCEPTION(problem_ == Teuchos::null,GCRODRSolMgrLinearProblemFailure, "Belos::GCRODRSolMgr::solve(): Linear problem is not a valid object.");
1344
1345 TEUCHOS_TEST_FOR_EXCEPTION(!problem_->isProblemSet(),GCRODRSolMgrLinearProblemFailure,"Belos::GCRODRSolMgr::solve(): Linear problem is not ready, setProblem() has not been called.");
1346
1347 // Create indices for the linear systems to be solved.
1348 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1349 std::vector<int> currIdx(1);
1350 currIdx[0] = 0;
1351
1352 // Inform the linear problem of the current linear system to solve.
1353 problem_->setLSIndex( currIdx );
1354
1355 // Check the number of blocks and change them is necessary.
1356 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1357 if (static_cast<ptrdiff_t>(numBlocks_) > dim) {
1358 numBlocks_ = Teuchos::as<int>(dim);
1359 printer_->stream(Warnings) <<
1360 "Warning! Requested Krylov subspace dimension is larger than operator dimension!" << std::endl <<
1361 " The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << numBlocks_ << std::endl;
1362 params_->set("Num Blocks", numBlocks_);
1363 }
1364
1365 // Assume convergence is achieved, then let any failed convergence set this to false.
1366 bool isConverged = true;
1367
1368 // Initialize storage for all state variables
1369 initializeStateStorage();
1370
1372 // Parameter list
1373 Teuchos::ParameterList plist;
1374
1375 plist.set("Num Blocks",numBlocks_);
1376 plist.set("Recycled Blocks",recycledBlocks_);
1377
1379 // GCRODR solver
1380
1381 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_iter;
1382 gcrodr_iter = rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,plist) );
1383 // Number of iterations required to generate initial recycle space (if needed)
1384 int prime_iterations = 0;
1385
1386 // Enter solve() iterations
1387 {
1388#ifdef BELOS_TEUCHOS_TIME_MONITOR
1389 Teuchos::TimeMonitor slvtimer(*timerSolve_);
1390#endif
1391
1392 while ( numRHS2Solve > 0 ) {
1393
1394 // Set flag indicating recycle space has not been generated this solve
1395 builtRecycleSpace_ = false;
1396
1397 // Reset the status test.
1398 outputTest_->reset();
1399
1401 // Initialize recycled subspace for GCRODR
1402
1403 // If there is a subspace to recycle, recycle it, otherwise generate the initial recycled subspace.
1404 if (keff > 0) {
1405 TEUCHOS_TEST_FOR_EXCEPTION(keff < recycledBlocks_,GCRODRSolMgrRecyclingFailure,
1406 "Belos::GCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1407
1408 printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
1409 // Compute image of U_ under the new operator
1410 index.resize(keff);
1411 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1412 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1413 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1414 problem_->apply( *Utmp, *Ctmp );
1415
1416 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1417
1418 // Orthogonalize this block
1419 // Get a matrix to hold the orthonormalization coefficients.
1420 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1421 int rank = ortho_->normalize(*Ctmp, rcp(&Rtmp,false));
1422 // Throw an error if we could not orthogonalize this block
1423 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,GCRODRSolMgrOrthoFailure,"Belos::GCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1424
1425 // U_ = U_*R^{-1}
1426 // First, compute LU factorization of R
1427 int info = 0;
1428 ipiv_.resize(Rtmp.numRows());
1429 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1430 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1431
1432 // Now, form inv(R)
1433 int lwork = Rtmp.numRows();
1434 work_.resize(lwork);
1435 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1436 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1437
1438 // U_ = U1_; (via a swap)
1439 MVT::MvTimesMatAddMv( one, *Utmp, Rtmp, zero, *U1tmp );
1440 std::swap(U_, U1_);
1441
1442 // Must reinitialize after swap
1443 index.resize(keff);
1444 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1445 Ctmp = MVT::CloneViewNonConst( *C_, index );
1446 Utmp = MVT::CloneView( *U_, index );
1447
1448 // Compute C_'*r_
1449 Teuchos::SerialDenseMatrix<int,ScalarType> Ctr(keff,1);
1450 problem_->computeCurrPrecResVec( &*r_ );
1451 MVT::MvTransMv( one, *Ctmp, *r_, Ctr );
1452
1453 // Update solution ( x += U_*C_'*r_ )
1454 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), 1 );
1455 MVT::MvInit( *update, 0.0 );
1456 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1457 problem_->updateSolution( update, true );
1458
1459 // Update residual norm ( r -= C_*C_'*r_ )
1460 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *r_ );
1461
1462 // We recycled space from previous call
1463 prime_iterations = 0;
1464
1465 }
1466 else {
1467
1468 // Do one cycle of Gmres to "prime the pump" if there is no subspace to recycle
1469 printer_->stream(Debug) << " No recycled subspace available for RHS index " << currIdx[0] << std::endl << std::endl;
1470
1471 Teuchos::ParameterList primeList;
1472
1473 // Tell the block solver that the block size is one.
1474 primeList.set("Num Blocks",numBlocks_);
1475 primeList.set("Recycled Blocks",0);
1476
1477 // Create GCRODR iterator object to perform one cycle of GMRES.
1478 RCP<GCRODRIter<ScalarType,MV,OP> > gcrodr_prime_iter;
1479 gcrodr_prime_iter = rcp( new GCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
1480
1481 // Create the first block in the current Krylov basis (residual).
1482 problem_->computeCurrPrecResVec( &*r_ );
1483 index.resize( 1 ); index[0] = 0;
1484 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1485 MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
1486
1487 // Set the new state and initialize the solver.
1489 index.resize( numBlocks_+1 );
1490 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1491 newstate.V = MVT::CloneViewNonConst( *V_, index );
1492 newstate.U = Teuchos::null;
1493 newstate.C = Teuchos::null;
1494 newstate.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, recycledBlocks_+1, recycledBlocks_+1 ) );
1495 newstate.B = Teuchos::null;
1496 newstate.curDim = 0;
1497 gcrodr_prime_iter->initialize(newstate);
1498
1499 // Perform one cycle of GMRES
1500 bool primeConverged = false;
1501 try {
1502 gcrodr_prime_iter->iterate();
1503
1504 // Check convergence first
1505 if ( convTest_->getStatus() == Passed ) {
1506 // we have convergence
1507 primeConverged = true;
1508 }
1509 }
1510 catch (const GCRODRIterOrthoFailure &e) {
1511 // Try to recover the most recent least-squares solution
1512 gcrodr_prime_iter->updateLSQR( gcrodr_prime_iter->getCurSubspaceDim() );
1513
1514 // Check to see if the most recent least-squares solution yielded convergence.
1515 sTest_->checkStatus( &*gcrodr_prime_iter );
1516 if (convTest_->getStatus() == Passed)
1517 primeConverged = true;
1518 }
1519 catch (const std::exception &e) {
1520 printer_->stream(Errors) << "Error! Caught exception in GCRODRIter::iterate() at iteration "
1521 << gcrodr_prime_iter->getNumIters() << std::endl
1522 << e.what() << std::endl;
1523 throw;
1524 }
1525 // Record number of iterations in generating initial recycle spacec
1526 prime_iterations = gcrodr_prime_iter->getNumIters();
1527
1528 // Update the linear problem.
1529 RCP<MV> update = gcrodr_prime_iter->getCurrentUpdate();
1530 problem_->updateSolution( update, true );
1531
1532 // Get the state.
1533 newstate = gcrodr_prime_iter->getState();
1534 int p = newstate.curDim;
1535
1536 // Compute harmonic Ritz vectors
1537 // NOTE: The storage for the harmonic Ritz vectors (PP) is made one column larger
1538 // just in case we split a complex conjugate pair.
1539 // NOTE: Generate a recycled subspace only if we have enough vectors. If we converged
1540 // too early, move on to the next linear system and try to generate a subspace again.
1541 if (recycledBlocks_ < p+1) {
1542 int info = 0;
1543 RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1544 // getHarmonicVecs1 assumes PP has recycledBlocks_+1 columns available
1545 keff = getHarmonicVecs1( p, *newstate.H, *PPtmp );
1546 // Hereafter, only keff columns of PP are needed
1547 PPtmp = rcp (new Teuchos::SerialDenseMatrix<int,ScalarType> ( Teuchos::View, *PP_, p, keff ) );
1548 // Now get views into C, U, V
1549 index.resize(keff);
1550 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1551 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1552 RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1553 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1554 index.resize(p);
1555 for (int ii=0; ii < p; ++ii) { index[ii] = ii; }
1556 RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1557
1558 // Form U (the subspace to recycle)
1559 // U = newstate.V(:,1:p) * PP;
1560 MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1561
1562 // Form orthonormalized C and adjust U so that C = A*U
1563
1564 // First, compute [Q, R] = qr(H*P);
1565
1566 // Step #1: Form HP = H*P
1567 Teuchos::SerialDenseMatrix<int,ScalarType> Htmp( Teuchos::View, *H2_, p+1, p, recycledBlocks_+1,recycledBlocks_+1);
1568 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+1, keff );
1569 HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, Htmp, *PPtmp, zero );
1570
1571 // Step #1.5: Perform workspace size query for QR
1572 // factorization of HP. On input, lwork must be -1.
1573 // _GEQRF will put the workspace size in work_[0].
1574 int lwork = -1;
1575 tau_.resize (keff);
1576 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1577 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1578 TEUCHOS_TEST_FOR_EXCEPTION(
1579 info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve:"
1580 " LAPACK's _GEQRF failed to compute a workspace size.");
1581
1582 // Step #2: Compute QR factorization of HP
1583 //
1584 // NOTE (mfh 17 Apr 2014) LAPACK promises that the value of
1585 // work_[0] after the workspace query will fit in int. This
1586 // justifies the cast. We call real() first because
1587 // static_cast from std::complex to int doesn't work.
1588 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1589 work_.resize (lwork); // Allocate workspace for the QR factorization
1590 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1591 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1592 TEUCHOS_TEST_FOR_EXCEPTION(
1593 info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve:"
1594 " LAPACK's _GEQRF failed to compute a QR factorization.");
1595
1596 // Step #3: Explicitly construct Q and R factors
1597 // NOTE: The upper triangular part of HP is copied into R and HP becomes Q.
1598 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff, keff );
1599 for (int ii = 0; ii < keff; ++ii) {
1600 for (int jj = ii; jj < keff; ++jj) {
1601 Rtmp(ii,jj) = HPtmp(ii,jj);
1602 }
1603 }
1604 // NOTE (mfh 17 Apr 2014): Teuchos::LAPACK's wrapper for
1605 // UNGQR dispatches to the correct Scalar-specific routine.
1606 // It calls {S,D}ORGQR if Scalar is real, and {C,Z}UNGQR if
1607 // Scalar is complex.
1608 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
1609 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
1610 lwork, &info);
1611 TEUCHOS_TEST_FOR_EXCEPTION(
1612 info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1613 "LAPACK's _UNGQR failed to construct the Q factor.");
1614
1615 // Now we have [Q,R] = qr(H*P)
1616
1617 // Now compute C = V(:,1:p+1) * Q
1618 index.resize (p + 1);
1619 for (int ii = 0; ii < (p+1); ++ii) {
1620 index[ii] = ii;
1621 }
1622 Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+1 vectors now; needed p above)
1623 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1624
1625 // Finally, compute U = U*R^{-1}.
1626 // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as
1627 // backsolve capabilities don't exist in the Belos::MultiVec class
1628
1629 // Step #1: First, compute LU factorization of R
1630 ipiv_.resize(Rtmp.numRows());
1631 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1632 TEUCHOS_TEST_FOR_EXCEPTION(
1633 info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1634 "LAPACK's _GETRF failed to compute an LU factorization.");
1635
1636 // FIXME (mfh 17 Apr 2014) We have to compute the explicit
1637 // inverse of R here because Belos::MultiVecTraits doesn't
1638 // have a triangular solve (where the triangular matrix is
1639 // globally replicated and the "right-hand side" is the
1640 // distributed MultiVector).
1641
1642 // Step #2: Form inv(R)
1643 lwork = Rtmp.numRows();
1644 work_.resize(lwork);
1645 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1646 TEUCHOS_TEST_FOR_EXCEPTION(
1647 info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1648 "LAPACK's _GETRI failed to invert triangular matrix.");
1649
1650 // Step #3: Let U = U * R^{-1}
1651 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1652
1653 printer_->stream(Debug)
1654 << " Generated recycled subspace using RHS index " << currIdx[0]
1655 << " of dimension " << keff << std::endl << std::endl;
1656
1657 } // if (recycledBlocks_ < p+1)
1658
1659 // Return to outer loop if the priming solve converged, set the next linear system.
1660 if (primeConverged) {
1661 // Inform the linear problem that we are finished with this block linear system.
1662 problem_->setCurrLS();
1663
1664 // Update indices for the linear systems to be solved.
1665 numRHS2Solve -= 1;
1666 if (numRHS2Solve > 0) {
1667 currIdx[0]++;
1668 problem_->setLSIndex (currIdx); // Set the next indices
1669 }
1670 else {
1671 currIdx.resize (numRHS2Solve);
1672 }
1673
1674 continue;
1675 }
1676 } // if (keff > 0) ...
1677
1678 // Prepare for the Gmres iterations with the recycled subspace.
1679
1680 // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
1681 gcrodr_iter->setSize( keff, numBlocks_ );
1682
1683 // Reset the number of iterations.
1684 gcrodr_iter->resetNumIters(prime_iterations);
1685
1686 // Reset the number of calls that the status test output knows about.
1687 outputTest_->resetNumCalls();
1688
1689 // Compute the residual after the priming solve, it will be the first block in the current Krylov basis.
1690 problem_->computeCurrPrecResVec( &*r_ );
1691 index.resize( 1 ); index[0] = 0;
1692 RCP<MV> v0 = MVT::CloneViewNonConst( *V_, index );
1693 MVT::SetBlock(*r_,index,*v0); // V(:,0) = r
1694
1695 // Set the new state and initialize the solver.
1697 index.resize( numBlocks_+1 );
1698 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1699 newstate.V = MVT::CloneViewNonConst( *V_, index );
1700 index.resize( keff );
1701 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1702 newstate.C = MVT::CloneViewNonConst( *C_, index );
1703 newstate.U = MVT::CloneViewNonConst( *U_, index );
1704 newstate.B = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1705 newstate.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1706 newstate.curDim = 0;
1707 gcrodr_iter->initialize(newstate);
1708
1709 // variables needed for inner loop
1710 int numRestarts = 0;
1711 while(1) {
1712
1713 // tell gcrodr_iter to iterate
1714 try {
1715 gcrodr_iter->iterate();
1716
1718 //
1719 // check convergence first
1720 //
1722 if ( convTest_->getStatus() == Passed ) {
1723 // we have convergence
1724 break; // break from while(1){gcrodr_iter->iterate()}
1725 }
1727 //
1728 // check for maximum iterations
1729 //
1731 else if ( maxIterTest_->getStatus() == Passed ) {
1732 // we don't have convergence
1733 isConverged = false;
1734 break; // break from while(1){gcrodr_iter->iterate()}
1735 }
1737 //
1738 // check for restarting, i.e. the subspace is full
1739 //
1741 else if ( gcrodr_iter->getCurSubspaceDim() == gcrodr_iter->getMaxSubspaceDim() ) {
1742
1743 // Update the recycled subspace even if we have hit the maximum number of restarts.
1744
1745 // Update the linear problem.
1746 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1747 problem_->updateSolution( update, true );
1748
1749 buildRecycleSpace2(gcrodr_iter);
1750
1751 printer_->stream(Debug)
1752 << " Generated new recycled subspace using RHS index "
1753 << currIdx[0] << " of dimension " << keff << std::endl
1754 << std::endl;
1755
1756 // NOTE: If we have hit the maximum number of restarts then quit
1757 if (numRestarts >= maxRestarts_) {
1758 isConverged = false;
1759 break; // break from while(1){gcrodr_iter->iterate()}
1760 }
1761 numRestarts++;
1762
1763 printer_->stream(Debug)
1764 << " Performing restart number " << numRestarts << " of "
1765 << maxRestarts_ << std::endl << std::endl;
1766
1767 // Create the restart vector (first block in the current Krylov basis)
1768 problem_->computeCurrPrecResVec( &*r_ );
1769 index.resize( 1 ); index[0] = 0;
1770 RCP<MV> v00 = MVT::CloneViewNonConst( *V_, index );
1771 MVT::SetBlock(*r_,index,*v00); // V(:,0) = r
1772
1773 // Set the new state and initialize the solver.
1774 GCRODRIterState<ScalarType,MV> restartState;
1775 index.resize( numBlocks_+1 );
1776 for (int ii=0; ii<(numBlocks_+1); ++ii) { index[ii] = ii; }
1777 restartState.V = MVT::CloneViewNonConst( *V_, index );
1778 index.resize( keff );
1779 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1780 restartState.U = MVT::CloneViewNonConst( *U_, index );
1781 restartState.C = MVT::CloneViewNonConst( *C_, index );
1782 restartState.B = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, keff, numBlocks_, 0, keff ) );
1783 restartState.H = rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, numBlocks_+1, numBlocks_, keff, keff ) );
1784 restartState.curDim = 0;
1785 gcrodr_iter->initialize(restartState);
1786
1787
1788 } // end of restarting
1789
1791 //
1792 // we returned from iterate(), but none of our status tests Passed.
1793 // something is wrong, and it is probably our fault.
1794 //
1796
1797 else {
1798 TEUCHOS_TEST_FOR_EXCEPTION(
1799 true, std::logic_error, "Belos::GCRODRSolMgr::solve: "
1800 "Invalid return from GCRODRIter::iterate().");
1801 }
1802 }
1803 catch (const GCRODRIterOrthoFailure &e) {
1804 // Try to recover the most recent least-squares solution
1805 gcrodr_iter->updateLSQR( gcrodr_iter->getCurSubspaceDim() );
1806
1807 // Check to see if the most recent least-squares solution yielded convergence.
1808 sTest_->checkStatus( &*gcrodr_iter );
1809 if (convTest_->getStatus() != Passed)
1810 isConverged = false;
1811 break;
1812 }
1813 catch (const std::exception& e) {
1814 printer_->stream(Errors)
1815 << "Error! Caught exception in GCRODRIter::iterate() at iteration "
1816 << gcrodr_iter->getNumIters() << std::endl << e.what() << std::endl;
1817 throw;
1818 }
1819 }
1820
1821 // Compute the current solution.
1822 // Update the linear problem.
1823 RCP<MV> update = gcrodr_iter->getCurrentUpdate();
1824 problem_->updateSolution( update, true );
1825
1826 // Inform the linear problem that we are finished with this block linear system.
1827 problem_->setCurrLS();
1828
1829 // If we didn't build a recycle space this solve but ran at least k iterations,
1830 // force build of new recycle space
1831
1832 if (!builtRecycleSpace_) {
1833 buildRecycleSpace2(gcrodr_iter);
1834 printer_->stream(Debug)
1835 << " Generated new recycled subspace using RHS index " << currIdx[0]
1836 << " of dimension " << keff << std::endl << std::endl;
1837 }
1838
1839 // Update indices for the linear systems to be solved.
1840 numRHS2Solve -= 1;
1841 if (numRHS2Solve > 0) {
1842 currIdx[0]++;
1843 problem_->setLSIndex (currIdx); // Set the next indices
1844 }
1845 else {
1846 currIdx.resize (numRHS2Solve);
1847 }
1848 } // while (numRHS2Solve > 0)
1849 }
1850
1851 sTest_->print (printer_->stream (FinalSummary)); // print final summary
1852
1853 // print timing information
1854#ifdef BELOS_TEUCHOS_TIME_MONITOR
1855 // Calling summarize() can be expensive, so don't call unless the
1856 // user wants to print out timing details. summarize() will do all
1857 // the work even if it's passed a "black hole" output stream.
1858 if (verbosity_ & TimingDetails)
1859 Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
1860#endif // BELOS_TEUCHOS_TIME_MONITOR
1861
1862 // get iteration information for this solve
1863 numIters_ = maxIterTest_->getNumIters ();
1864
1865 // Save the convergence test value ("achieved tolerance") for this
1866 // solve. This solver (unlike BlockGmresSolMgr) always has two
1867 // residual norm status tests: an explicit and an implicit test.
1868 // The master convergence test convTest_ is a SEQ combo of the
1869 // implicit resp. explicit tests. If the implicit test never
1870 // passes, then the explicit test won't ever be executed. This
1871 // manifests as expConvTest_->getTestValue()->size() < 1. We deal
1872 // with this case by using the values returned by
1873 // impConvTest_->getTestValue().
1874 {
1875 const std::vector<MagnitudeType>* pTestValues = expConvTest_->getTestValue();
1876 if (pTestValues == NULL || pTestValues->size() < 1) {
1877 pTestValues = impConvTest_->getTestValue();
1878 }
1879 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
1880 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1881 "method returned NULL. Please report this bug to the Belos developers.");
1882 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
1883 "Belos::GCRODRSolMgr::solve(): The implicit convergence test's getTestValue() "
1884 "method returned a vector of length zero. Please report this bug to the "
1885 "Belos developers.");
1886
1887 // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
1888 // achieved tolerances for all vectors in the current solve(), or
1889 // just for the vectors from the last deflation?
1890 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
1891 }
1892
1893 return isConverged ? Converged : Unconverged; // return from solve()
1894}
1895
1896// Given existing recycle space and Krylov space, build new recycle space
1897template<class ScalarType, class MV, class OP>
1899
1900 MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
1901 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1902
1903 std::vector<MagnitudeType> d(keff);
1904 std::vector<ScalarType> dscalar(keff);
1905 std::vector<int> index(numBlocks_+1);
1906
1907 // Get the state
1908 GCRODRIterState<ScalarType,MV> oldState = gcrodr_iter->getState();
1909 int p = oldState.curDim;
1910
1911 // insufficient new information to update recycle space
1912 if (p<1) return;
1913
1914 // Take the norm of the recycled vectors
1915 {
1916 index.resize(keff);
1917 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1918 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1919 d.resize(keff);
1920 dscalar.resize(keff);
1921 MVT::MvNorm( *Utmp, d );
1922 for (int i=0; i<keff; ++i) {
1923 d[i] = one / d[i];
1924 dscalar[i] = (ScalarType)d[i];
1925 }
1926 MVT::MvScale( *Utmp, dscalar );
1927 }
1928
1929 // Get view into current "full" upper Hessnburg matrix
1930 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H2tmp = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>( Teuchos::View, *H2_, p+keff+1, p+keff ) );
1931
1932 // Insert D into the leading keff x keff block of H2
1933 for (int i=0; i<keff; ++i) {
1934 (*H2tmp)(i,i) = d[i];
1935 }
1936
1937 // Compute the harmoic Ritz pairs for the generalized eigenproblem
1938 // getHarmonicVecs2 assumes PP has recycledBlocks_+1 columns available
1939 int keff_new;
1940 {
1941 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, recycledBlocks_+1 );
1942 keff_new = getHarmonicVecs2( keff, p, *H2tmp, oldState.V, PPtmp );
1943 }
1944
1945 // Code to form new U, C
1946 // U = [U V(:,1:p)] * P; (in two steps)
1947
1948 // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
1949 Teuchos::RCP<MV> U1tmp;
1950 {
1951 index.resize( keff );
1952 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1953 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1954 index.resize( keff_new );
1955 for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1956 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1957 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, keff, keff_new );
1958 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1959 }
1960
1961 // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
1962 {
1963 index.resize(p);
1964 for (int ii=0; ii < p; ii++) { index[ii] = ii; }
1965 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1966 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p, keff_new, keff );
1967 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1968 }
1969
1970 // Form HP = H*P
1971 Teuchos::SerialDenseMatrix<int,ScalarType> HPtmp( Teuchos::View, *HP_, p+keff+1, keff_new );
1972 {
1973 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *PP_, p+keff, keff_new );
1974 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*H2tmp,PPtmp,zero);
1975 }
1976
1977 // Workspace size query for QR factorization of HP (the worksize will be placed in work_[0])
1978 int info = 0, lwork = -1;
1979 tau_.resize (keff_new);
1980 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1981 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1982 TEUCHOS_TEST_FOR_EXCEPTION(
1983 info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1984 "LAPACK's _GEQRF failed to compute a workspace size.");
1985
1986 // NOTE (mfh 18 Apr 2014) LAPACK promises that the value of work_[0]
1987 // after the workspace query will fit in int. This justifies the
1988 // cast. We call real() first because static_cast from std::complex
1989 // to int doesn't work.
1990 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0])));
1991 work_.resize (lwork); // Allocate workspace for the QR factorization
1992 lapack.GEQRF (HPtmp.numRows (), HPtmp.numCols (), HPtmp.values (),
1993 HPtmp.stride (), &tau_[0], &work_[0], lwork, &info);
1994 TEUCHOS_TEST_FOR_EXCEPTION(
1995 info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
1996 "LAPACK's _GEQRF failed to compute a QR factorization.");
1997
1998 // Explicitly construct Q and R factors
1999 // NOTE: The upper triangular part of HP is copied into R and HP becomes Q.
2000 Teuchos::SerialDenseMatrix<int,ScalarType> Rtmp( Teuchos::View, *R_, keff_new, keff_new );
2001 for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Rtmp(i,j) = HPtmp(i,j); }
2002
2003 // NOTE (mfh 18 Apr 2014): Teuchos::LAPACK's wrapper for UNGQR
2004 // dispatches to the correct Scalar-specific routine. It calls
2005 // {S,D}ORGQR if Scalar is real, and {C,Z}UNGQR if Scalar is
2006 // complex.
2007 lapack.UNGQR (HPtmp.numRows (), HPtmp.numCols (), HPtmp.numCols (),
2008 HPtmp.values (), HPtmp.stride (), &tau_[0], &work_[0],
2009 lwork, &info);
2010 TEUCHOS_TEST_FOR_EXCEPTION(
2011 info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve: "
2012 "LAPACK's _UNGQR failed to construct the Q factor.");
2013
2014 // Form orthonormalized C and adjust U accordingly so that C = A*U
2015 // C = [C V] * Q;
2016
2017 // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
2018 {
2019 Teuchos::RCP<MV> C1tmp;
2020 {
2021 index.resize(keff);
2022 for (int i=0; i < keff; i++) { index[i] = i; }
2023 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2024 index.resize(keff_new);
2025 for (int i=0; i < keff_new; i++) { index[i] = i; }
2026 C1tmp = MVT::CloneViewNonConst( *C1_, index );
2027 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, keff, keff_new );
2028 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
2029 }
2030 // Now compute C += V(:,1:p+1) * Q
2031 {
2032 index.resize( p+1 );
2033 for (int i=0; i < p+1; ++i) { index[i] = i; }
2034 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
2035 Teuchos::SerialDenseMatrix<int,ScalarType> PPtmp( Teuchos::View, *HP_, p+1, keff_new, keff, 0 );
2036 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
2037 }
2038 }
2039
2040 // C_ = C1_; (via a swap)
2041 std::swap(C_, C1_);
2042
2043 // Finally, compute U_ = U_*R^{-1}
2044 // First, compute LU factorization of R
2045 ipiv_.resize(Rtmp.numRows());
2046 lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
2047 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
2048
2049 // Now, form inv(R)
2050 lwork = Rtmp.numRows();
2051 work_.resize(lwork);
2052 lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
2053 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
2054
2055 {
2056 index.resize(keff_new);
2057 for (int i=0; i < keff_new; i++) { index[i] = i; }
2058 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
2059 MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
2060 }
2061
2062 // Set the current number of recycled blocks and subspace dimension with the GCRO-DR iteration.
2063 if (keff != keff_new) {
2064 keff = keff_new;
2065 gcrodr_iter->setSize( keff, numBlocks_ );
2066 // Important to zero this out before next cyle
2067 Teuchos::SerialDenseMatrix<int,ScalarType> b1( Teuchos::View, *H2_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
2068 b1.putScalar(zero);
2069 }
2070
2071}
2072
2073
2074// Compute the harmonic eigenpairs of the projected, dense system.
2075template<class ScalarType, class MV, class OP>
2077 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2078 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2079 int i, j;
2080 bool xtraVec = false;
2081 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2082
2083 // Real and imaginary eigenvalue components
2084 std::vector<MagnitudeType> wr(m), wi(m);
2085
2086 // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
2087 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m,m,false);
2088
2089 // Magnitude of harmonic Ritz values
2090 std::vector<MagnitudeType> w(m);
2091
2092 // Sorted order of harmonic Ritz values, also used for GEEV
2093 std::vector<int> iperm(m);
2094
2095 // Output info
2096 int info = 0;
2097
2098 // Set flag indicating recycle space has been generated this solve
2099 builtRecycleSpace_ = true;
2100
2101 // Solve linear system: H_m^{-H}*e_m
2102 Teuchos::SerialDenseMatrix<int, ScalarType> HHt( HH, Teuchos::TRANS );
2103 Teuchos::SerialDenseVector<int, ScalarType> e_m( m );
2104 e_m[m-1] = one;
2105 lapack.GESV(m, 1, HHt.values(), HHt.stride(), &iperm[0], e_m.values(), e_m.stride(), &info);
2106 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
2107
2108 // Compute H_m + d*H_m^{-H}*e_m*e_m^H
2109 ScalarType d = HH(m, m-1) * HH(m, m-1);
2110 Teuchos::SerialDenseMatrix<int, ScalarType> harmHH( Teuchos::Copy, HH, m, m );
2111 for( i=0; i<m; ++i )
2112 harmHH(i, m-1) += d * e_m[i];
2113
2114 // Revise to do query for optimal workspace first
2115 // Create simple storage for the left eigenvectors, which we don't care about.
2116 const int ldvl = 1;
2117 ScalarType* vl = 0;
2118
2119 // Size of workspace and workspace for GEEV
2120 int lwork = -1;
2121 std::vector<ScalarType> work(1);
2122 std::vector<MagnitudeType> rwork(2*m);
2123
2124 // First query GEEV for the optimal workspace size
2125 lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2126 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2127
2128 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
2129 work.resize( lwork );
2130
2131 lapack.GEEV('N', 'V', m, harmHH.values(), harmHH.stride(), &wr[0], &wi[0],
2132 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
2133 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure,"Belos::GCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
2134
2135 // Construct magnitude of each harmonic Ritz value
2136 for( i=0; i<m; ++i )
2137 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
2138
2139 // Construct magnitude of each harmonic Ritz value
2140 this->sort(w, m, iperm);
2141
2142 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2143
2144 // Select recycledBlocks_ smallest eigenvectors
2145 for( i=0; i<recycledBlocks_; ++i ) {
2146 for( j=0; j<m; j++ ) {
2147 PP(j,i) = vr(j,iperm[i]);
2148 }
2149 }
2150
2151 if(!scalarTypeIsComplex) {
2152
2153 // Determine exact size for PP (i.e., determine if we need to store an additional vector)
2154 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
2155 int countImag = 0;
2156 for ( i=0; i<recycledBlocks_; ++i ) {
2157 if (wi[iperm[i]] != 0.0)
2158 countImag++;
2159 }
2160 // Check to see if this count is even or odd:
2161 if (countImag % 2)
2162 xtraVec = true;
2163 }
2164
2165 if (xtraVec) { // we need to store one more vector
2166 if (wi[iperm[recycledBlocks_-1]] > 0.0) { // I picked the "real" component
2167 for( j=0; j<m; ++j ) { // so get the "imag" component
2168 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
2169 }
2170 }
2171 else { // I picked the "imag" component
2172 for( j=0; j<m; ++j ) { // so get the "real" component
2173 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
2174 }
2175 }
2176 }
2177
2178 }
2179
2180 // Return whether we needed to store an additional vector
2181 if (xtraVec) {
2182 return recycledBlocks_+1;
2183 }
2184 else {
2185 return recycledBlocks_;
2186 }
2187
2188}
2189
2190// Compute the harmonic eigenpairs of the projected, dense system.
2191template<class ScalarType, class MV, class OP>
2193 const Teuchos::SerialDenseMatrix<int,ScalarType>& HH,
2194 const Teuchos::RCP<const MV>& VV,
2195 Teuchos::SerialDenseMatrix<int,ScalarType>& PP) {
2196 int i, j;
2197 int m2 = HH.numCols();
2198 bool xtraVec = false;
2199 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
2200 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
2201 std::vector<int> index;
2202
2203 // Real and imaginary eigenvalue components
2204 std::vector<MagnitudeType> wr(m2), wi(m2);
2205
2206 // Magnitude of harmonic Ritz values
2207 std::vector<MagnitudeType> w(m2);
2208
2209 // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
2210 Teuchos::SerialDenseMatrix<int,ScalarType> vr(m2,m2,false);
2211
2212 // Sorted order of harmonic Ritz values
2213 std::vector<int> iperm(m2);
2214
2215 // Set flag indicating recycle space has been generated this solve
2216 builtRecycleSpace_ = true;
2217
2218 // Form matrices for generalized eigenproblem
2219
2220 // B = H2' * H2; Don't zero out matrix when constructing
2221 Teuchos::SerialDenseMatrix<int,ScalarType> B(m2,m2,false);
2222 B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,HH,HH,zero);
2223
2224 // A_tmp = | C'*U 0 |
2225 // | V_{m+1}'*U I |
2226 Teuchos::SerialDenseMatrix<int,ScalarType> A_tmp( keffloc+m+1, keffloc+m );
2227
2228 // A_tmp(1:keffloc,1:keffloc) = C' * U;
2229 index.resize(keffloc);
2230 for (i=0; i<keffloc; ++i) { index[i] = i; }
2231 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
2232 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
2233 Teuchos::SerialDenseMatrix<int,ScalarType> A11( Teuchos::View, A_tmp, keffloc, keffloc );
2234 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
2235
2236 // A_tmp(keffloc+1:m-k+keffloc+1,1:keffloc) = V' * U;
2237 Teuchos::SerialDenseMatrix<int,ScalarType> A21( Teuchos::View, A_tmp, m+1, keffloc, keffloc );
2238 index.resize(m+1);
2239 for (i=0; i < m+1; i++) { index[i] = i; }
2240 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
2241 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
2242
2243 // A_tmp(keffloc+1:m-k+keffloc,keffloc+1:m-k+keffloc) = eye(m-k);
2244 for( i=keffloc; i<keffloc+m; i++ ) {
2245 A_tmp(i,i) = one;
2246 }
2247
2248 // A = H2' * A_tmp;
2249 Teuchos::SerialDenseMatrix<int,ScalarType> A( m2, A_tmp.numCols() );
2250 A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, HH, A_tmp, zero );
2251
2252 // Compute k smallest harmonic Ritz pairs
2253 // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
2254 // ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
2255 // IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
2256 // RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
2257 // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
2258 char balanc='P', jobvl='N', jobvr='V', sense='N';
2259 int ld = A.numRows();
2260 int lwork = 6*ld;
2261 int ldvl = ld, ldvr = ld;
2262 int info = 0,ilo = 0,ihi = 0;
2263 MagnitudeType abnrm = 0.0, bbnrm = 0.0;
2264 ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
2265 std::vector<ScalarType> beta(ld);
2266 std::vector<ScalarType> work(lwork);
2267 std::vector<MagnitudeType> rwork(lwork);
2268 std::vector<MagnitudeType> lscale(ld), rscale(ld);
2269 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
2270 std::vector<int> iwork(ld+6);
2271 int *bwork = 0; // If sense == 'N', bwork is never referenced
2272 //lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
2273 // &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2274 // &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
2275 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
2276 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
2277 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
2278 &iwork[0], bwork, &info);
2279 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, GCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
2280
2281 // Construct magnitude of each harmonic Ritz value
2282 // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
2283 for( i=0; i<ld; i++ ) {
2284 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot (wr[i]*wr[i] + wi[i]*wi[i]) /
2285 Teuchos::ScalarTraits<ScalarType>::magnitude (beta[i]);
2286 }
2287
2288 // Construct magnitude of each harmonic Ritz value
2289 this->sort(w,ld,iperm);
2290
2291 const bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
2292
2293 // Select recycledBlocks_ smallest eigenvectors
2294 for( i=0; i<recycledBlocks_; i++ ) {
2295 for( j=0; j<ld; j++ ) {
2296 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
2297 }
2298 }
2299
2300 if(!scalarTypeIsComplex) {
2301
2302 // Determine exact size for PP (i.e., determine if we need to store an additional vector)
2303 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
2304 int countImag = 0;
2305 for ( i=ld-recycledBlocks_; i<ld; i++ ) {
2306 if (wi[iperm[i]] != 0.0)
2307 countImag++;
2308 }
2309 // Check to see if this count is even or odd:
2310 if (countImag % 2)
2311 xtraVec = true;
2312 }
2313
2314 if (xtraVec) { // we need to store one more vector
2315 if (wi[iperm[ld-recycledBlocks_]] > 0.0) { // I picked the "real" component
2316 for( j=0; j<ld; j++ ) { // so get the "imag" component
2317 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
2318 }
2319 }
2320 else { // I picked the "imag" component
2321 for( j=0; j<ld; j++ ) { // so get the "real" component
2322 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
2323 }
2324 }
2325 }
2326
2327 }
2328
2329 // Return whether we needed to store an additional vector
2330 if (xtraVec) {
2331 return recycledBlocks_+1;
2332 }
2333 else {
2334 return recycledBlocks_;
2335 }
2336
2337}
2338
2339
2340// This method sorts list of n floating-point numbers and return permutation vector
2341template<class ScalarType, class MV, class OP>
2342void GCRODRSolMgr<ScalarType,MV,OP,true>::sort(std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm) {
2343 int l, r, j, i, flag;
2344 int RR2;
2345 MagnitudeType dRR, dK;
2346
2347 // Initialize the permutation vector.
2348 for(j=0;j<n;j++)
2349 iperm[j] = j;
2350
2351 if (n <= 1) return;
2352
2353 l = n / 2 + 1;
2354 r = n - 1;
2355 l = l - 1;
2356 dRR = dlist[l - 1];
2357 dK = dlist[l - 1];
2358
2359 RR2 = iperm[l - 1];
2360 while (r != 0) {
2361 j = l;
2362 flag = 1;
2363
2364 while (flag == 1) {
2365 i = j;
2366 j = j + j;
2367
2368 if (j > r + 1)
2369 flag = 0;
2370 else {
2371 if (j < r + 1)
2372 if (dlist[j] > dlist[j - 1]) j = j + 1;
2373
2374 if (dlist[j - 1] > dK) {
2375 dlist[i - 1] = dlist[j - 1];
2376 iperm[i - 1] = iperm[j - 1];
2377 }
2378 else {
2379 flag = 0;
2380 }
2381 }
2382 }
2383 dlist[i - 1] = dRR;
2384 iperm[i - 1] = RR2;
2385
2386 if (l == 1) {
2387 dRR = dlist [r];
2388 RR2 = iperm[r];
2389 dK = dlist[r];
2390 dlist[r] = dlist[0];
2391 iperm[r] = iperm[0];
2392 r = r - 1;
2393 }
2394 else {
2395 l = l - 1;
2396 dRR = dlist[l - 1];
2397 RR2 = iperm[l - 1];
2398 dK = dlist[l - 1];
2399 }
2400 }
2401 dlist[0] = dRR;
2402 iperm[0] = RR2;
2403}
2404
2405
2406template<class ScalarType, class MV, class OP>
2408 std::ostringstream out;
2409 out << "Belos::GCRODRSolMgr<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
2410 out << "{";
2411 out << "Ortho Type: \"" << orthoType_ << "\"";
2412 out << ", Num Blocks: " <<numBlocks_;
2413 out << ", Num Recycle Blocks: " << recycledBlocks_;
2414 out << ", Max Restarts: " << maxRestarts_;
2415 out << "}";
2416 return out.str ();
2417}
2418
2419} // namespace Belos
2420
2421#endif /* BELOS_GCRODR_SOLMGR_HPP */
Belos concrete class for performing the block, flexible GMRES iteration.
Belos header file which uses auto-configuration information to include necessary C++ headers.
Belos concrete class for performing the GCRO-DR iteration.
Class which describes the linear problem to be solved by the iterative solver.
Class which manages the output and verbosity of the Belos solvers.
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
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
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 ScalarType types ...
GCRODRIterOrthoFailure is thrown when the GCRODRIter object is unable to compute independent directio...
This class implements the GCRODR iteration, where a single-std::vector Krylov subspace is constructed...
GCRODRSolMgrLAPACKFailure is thrown when a nonzero value is retuned from an LAPACK call.
GCRODRSolMgrLAPACKFailure(const std::string &what_arg)
GCRODRSolMgrLinearProblemFailure is thrown when the linear problem is not setup (i....
GCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
GCRODRSolMgrOrthoFailure is thrown when the orthogonalization manager is unable to generate orthonorm...
GCRODRSolMgrOrthoFailure(const std::string &what_arg)
GCRODRSolMgrRecyclingFailure is thrown when any problem occurs in using/creating the recycling subspa...
GCRODRSolMgrRecyclingFailure(const std::string &what_arg)
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
const LinearProblem< ScalarType, MV, OP > & getProblem() const override
Get current linear problem being solved for in this object.
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< SolverManager< ScalarType, MV, OP > > clone() const override
clone for Inverted Injection (DII)
Teuchos::Array< Teuchos::RCP< Teuchos::Time > > getTimers() const
Return the timers for this object.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H_
Teuchos::RCP< OutputManager< ScalarType > > printer_
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
OrthoManagerFactory< ScalarType, MV, OP > ortho_factory_type
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > PP_
int getNumIters() const override
Get the iteration count for the most recent call to solve().
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::ScalarTraits< MagnitudeType > MT
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem) override
Set the linear problem that needs to be solved.
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > expConvTest_
Teuchos::RCP< Teuchos::ParameterList > params_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H2_
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const override
Get a parameter list containing the current parameters for this object.
MagnitudeType achievedTol() const override
Tolerance achieved by the last solve() invocation.
bool isLOADetected() const override
Return whether a loss of accuracy was detected by this solver during the most current solve.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > R_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B_
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > HP_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
Implementation of the GCRODR (Recycling GMRES) iterative linear solver.
static const bool requiresLapack
Details::SolverManagerRequiresLapack< ScalarType, MV, OP, requiresLapack > base_type
GCRODRSolMgr(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Teuchos::ParameterList > &pl)
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.
Enumeration of all valid Belos (Mat)OrthoManager classes.
std::ostream & printValidNames(std::ostream &out) const
Print all recognized MatOrthoManager names to the given ostream.
bool isValidName(const std::string &name) const
Whether this factory recognizes the MatOrthoManager with the given name.
std::string validNamesString() const
List (as a string) of recognized MatOrthoManager names.
Teuchos::RCP< Belos::MatOrthoManager< Scalar, MV, OP > > makeMatOrthoManager(const std::string &ortho, const Teuchos::RCP< const OP > &M, const Teuchos::RCP< OutputManager< Scalar > > &, const std::string &label, const Teuchos::RCP< Teuchos::ParameterList > &params)
Return an instance of the specified MatOrthoManager subclass.
Teuchos::RCP< const Teuchos::ParameterList > getDefaultParameters(const std::string &name) const
Default parameters for the given MatOrthoManager subclass.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A class for extending the status testing capabilities of Belos via logical combinations.
Exception thrown to signal error in a status test during Belos::StatusTest::checkStatus().
An implementation of StatusTestResNorm using a family of residual norms.
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.
ScaleType convertStringToScaleType(const std::string &scaleType)
Convert the given string to its ScaleType enum value.
Definition: BelosTypes.cpp:106
@ TwoNorm
Definition: BelosTypes.hpp:98
@ 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
ScaleType
The type of scaling to use on the residual norm value.
Definition: BelosTypes.hpp:120
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 GCRODRIter state variables.
Teuchos::RCP< MV > C
Teuchos::RCP< MV > U
The recycled subspace and its projection.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > B
The projection of the Krylov subspace against the recycled subspace.
Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.
Teuchos::RCP< MV > V
The current Krylov basis.
int curDim
The current dimension of the reduction.