Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosBlockGCRODRSolMgr.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
45#ifndef BELOS_BLOCK_GCRODR_SOLMGR_HPP
46#define BELOS_BLOCK_GCRODR_SOLMGR_HPP
47
48#include "BelosConfigDefs.hpp"
49#include "BelosTypes.hpp"
62#include "Teuchos_LAPACK.hpp"
63#include "Teuchos_as.hpp"
64#ifdef BELOS_TEUCHOS_TIME_MONITOR
65#include "Teuchos_TimeMonitor.hpp"
66#endif // BELOS_TEUCHOS_TIME_MONITOR
67
68// MLP: Add links to examples here
69
70namespace Belos{
71
73
74
80 public:
81 BlockGCRODRSolMgrLinearProblemFailure(const std::string& what_arg) : BelosError(what_arg) {}
82 };
83
89 public:
90 BlockGCRODRSolMgrOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
91 };
92
98 public:
99 BlockGCRODRSolMgrLAPACKFailure(const std::string& what_arg) : BelosError(what_arg) {}
100 };
101
107 public:
108 BlockGCRODRSolMgrRecyclingFailure(const std::string& what_arg) : BelosError(what_arg) {}
109 };
110
112
124
125template<class ScalarType, class MV, class OP>
126class BlockGCRODRSolMgr : public SolverManager<ScalarType, MV, OP> {
127private:
128
131 typedef Teuchos::ScalarTraits<ScalarType> SCT;
132 typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;
133 typedef Teuchos::ScalarTraits<MagnitudeType> MT;
134 typedef Teuchos::ScalarTraits<MagnitudeType> SMT;
136 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SDM;
137 typedef Teuchos::SerialDenseVector<int,ScalarType> SDV;
138
139public:
141
142
149
175 BlockGCRODRSolMgr (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
176 const Teuchos::RCP<Teuchos::ParameterList> &pl);
177
179 virtual ~BlockGCRODRSolMgr() {};
181
184
186 std::string description() const;
187
189
190
192
193
196 return *problem_;
197 }
198
200 Teuchos::RCP<const Teuchos::ParameterList> getValidParameters() const;
201
203 Teuchos::RCP<const Teuchos::ParameterList> getCurrentParameters() const { return params_; }
204
207 return achievedTol_;
208 }
209
211 int getNumIters() const { return numIters_; }
212
214 bool isLOADetected() const { return loaDetected_; }
215
217
219
220
222 void setProblem (const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >& problem) {
223 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
224 "Belos::BlockGCRODRSolMgr::setProblem: The input LinearProblem cannot be null.");
225
226 // Check state of problem object before proceeding
227 if (! problem->isProblemSet()) {
228 const bool success = problem->setProblem();
229 TEUCHOS_TEST_FOR_EXCEPTION(success, std::runtime_error,
230 "Belos::BlockGCRODRSolMgr::setProblem: Calling the input LinearProblem's setProblem() method failed. This likely means that the "
231 "LinearProblem has a missing (null) matrix A, solution vector X, or right-hand side vector B. Please set these items in the LinearProblem and try again.");
232 }
233
234 problem_ = problem;
235 }
236
238 void setParameters( const Teuchos::RCP<Teuchos::ParameterList> &params );
239
241
243
244
251 void reset (const ResetType type) {
252 if ((type & Belos::Problem) && ! problem_.is_null ())
253 problem_->setProblem();
254 else if (type & Belos::RecycleSubspace)
255 keff = 0;
256 }
257
259
261
262
278 // \returns ::ReturnType specifying:
282
284
285private:
286
287 // Called by all constructors; Contains init instructions common to all constructors
288 void init();
289
290 // Initialize solver state storage
292
293 // Recycling Methods
294 // Appending Function name by:
295 // "Kryl" indicates it is specialized for building a recycle space after an
296 // initial run of Block GMRES which generates an initial unaugmented block Krylov subspace
297 // "AugKryl" indicates it is specialized for building a recycle space from the augmented Krylov subspace
298
299 // Functions which control the building of a recycle space
300 void buildRecycleSpaceKryl(int& keff, Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter);
301 void buildRecycleSpaceAugKryl(Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > gcrodr_iter);
302
303 // Recycling with Harmonic Ritz Vectors
304 // Computes harmonic eigenpairs of projected matrix created during the priming solve.
305 // The return value is the number of vectors needed to be stored, recycledBlocks or recycledBlocks+1.
306
307 // HH is the projected problem from the initial cycle of Gmres, it is (at least) of dimension (numBlocks+1)*blockSize x numBlocks.
308 // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
309 int getHarmonicVecsKryl (int m, const SDM& HH, SDM& PP);
310
311 // HH is the total block projected problem from the GCRO-DR algorithm, it is (at least) of dimension keff+(numBlocks+1)*blockSize x keff+numBlocksm.
312 // VV is the Krylov vectors from the projected GMRES algorithm, which has (at least) (numBlocks+1)*blockSize vectors.
313 // PP contains the harmonic eigenvectors corresponding to the recycledBlocks eigenvalues of smallest magnitude.
315 int m,
316 const SDM& HH,
317 const Teuchos::RCP<const MV>& VV,
318 SDM& PP);
319
320 // Sort list of n floating-point numbers and return permutation vector
321 void sort (std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm);
322
323 // Lapack interface
324 Teuchos::LAPACK<int,ScalarType> lapack;
325
327 Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > problem_;
328
329 //Output Manager
330 Teuchos::RCP<OutputManager<ScalarType> > printer_;
331 Teuchos::RCP<std::ostream> outputStream_;
332
333 //Status Test
334 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > sTest_;
335 Teuchos::RCP<StatusTestMaxIters<ScalarType,MV,OP> > maxIterTest_;
336 Teuchos::RCP<StatusTest<ScalarType,MV,OP> > convTest_;
337 Teuchos::RCP<StatusTestGenResNorm<ScalarType,MV,OP> > expConvTest_, impConvTest_;
338 Teuchos::RCP<StatusTestOutput<ScalarType,MV,OP> > outputTest_;
339
342
348 Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > ortho_;
349
351 Teuchos::RCP<Teuchos::ParameterList> params_;
352
363 mutable Teuchos::RCP<const Teuchos::ParameterList> defaultParams_;
364
365 //Default Solver Values
366 static const bool adaptiveBlockSize_default_;
367 static const std::string recycleMethod_default_;
368
369 //Current Solver Values
376 std::string label_;
377
379 // Solver State Storage
381 //
382 // The number of blocks and recycle blocks (m and k, respectively)
384 // Current size of recycled subspace
385 int keff;
386 //
387 // Residual Vector
388 Teuchos::RCP<MV> R_;
389 //
390 // Search Space
391 Teuchos::RCP<MV> V_;
392 //
393 // Recycle subspace and its image under action of the operator
394 Teuchos::RCP<MV> U_, C_;
395 //
396 // Updated recycle Space and its image under action of the operator
397 Teuchos::RCP<MV> U1_, C1_;
398 //
399 // Storage used in constructing recycle space
400 Teuchos::RCP<SDM > G_;
401 Teuchos::RCP<SDM > H_;
402 Teuchos::RCP<SDM > B_;
403 Teuchos::RCP<SDM > PP_;
404 Teuchos::RCP<SDM > HP_;
405 std::vector<ScalarType> tau_;
406 std::vector<ScalarType> work_;
407 Teuchos::RCP<SDM > F_;
408 std::vector<int> ipiv_;
409
411 Teuchos::RCP<Teuchos::Time> timerSolve_;
412
414 bool isSet_;
415
418
421};
422
423 //
424 // Set default solver values
425 //
426 template<class ScalarType, class MV, class OP>
428
429 template<class ScalarType, class MV, class OP>
431
432 //
433 // Method definitions
434 //
435
436 template<class ScalarType, class MV, class OP>
440
441 //Basic Constructor
442 template<class ScalarType, class MV, class OP>
444 BlockGCRODRSolMgr(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
445 const Teuchos::RCP<Teuchos::ParameterList> &pl ) {
446 // Initialize local pointers to null, and initialize local
447 // variables to default values.
448 init ();
449
450 TEUCHOS_TEST_FOR_EXCEPTION(problem.is_null(), std::invalid_argument,
451 "Belos::BlockGCRODR constructor: The solver manager's constructor needs "
452 "the linear problem argument 'problem' to be nonnull.");
453
454 problem_ = problem;
455
456 // Set the parameters using the list that was passed in. If null, we defer initialization until a non-null list is set (by the
457 // client calling setParameters(), or by calling solve() -- in either case, a null parameter list indicates that default parameters should be used).
458 if (! pl.is_null())
459 setParameters (pl);
460 }
461
462 template<class ScalarType, class MV, class OP>
464 adaptiveBlockSize_ = adaptiveBlockSize_default_;
465 recycleMethod_ = recycleMethod_default_;
466 isSet_ = false;
467 loaDetected_ = false;
468 builtRecycleSpace_ = false;
469 keff = 0;//Effective Size of Recycle Space
470 //The following are all RCP smart pointers to indicated matrices/vectors.
471 //Some MATLAB notation used in comments.
472 R_ = Teuchos::null;//The Block Residual
473 V_ = Teuchos::null;//Block Arnoldi Vectors
474 U_ = Teuchos::null;//Recycle Space
475 C_ = Teuchos::null;//Image of U Under Action of Operator
476 U1_ = Teuchos::null;//Newly Computed Recycle Space
477 C1_ = Teuchos::null;//Image of Newly Computed Recycle Space
478 PP_ = Teuchos::null;//Coordinates of New Recyc. Vectors in Augmented Space
479 HP_ = Teuchos::null;//H_*PP_ or G_*PP_
480 G_ = Teuchos::null;//G_ such that A*[U V(:,1:numBlocks_*blockSize_)] = [C V_]*G_
481 F_ = Teuchos::null;//Upper Triangular portion of QR factorization for HP_
482 H_ = Teuchos::null;//H_ such that A*V(:,1:numBlocks_*blockSize_) = V_*H_ + C_*C_'*V_
483 B_ = Teuchos::null;//B_ = C_'*V_
484
485 //THIS BLOCK OF CODE IS COMMENTED OUT AND PLACED ELSEWHERE IN THE CODE
486/*
487 //WE TEMPORARILY INITIALIZE STATUS TESTS HERE FOR TESTING PURPOSES, BUT THEY SHOULD BE
488 //INITIALIZED SOMEWHERE ELSE, LIKE THE setParameters() FUNCTION
489
490 //INSTANTIATE AND INITIALIZE TEST OBJECTS AS NEEDED
491 if (maxIterTest_.is_null()){
492 maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
493 }
494 //maxIterTest_->setMaxIters (maxIters_);
495
496 //INSTANTIATE THE PRINTER
497 if (printer_.is_null()) {
498 printer_ = Teuchos::rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
499 }
500
501 if (ortho_.is_null()) // || changedOrthoType) %%%In other codes, this is also triggered if orthogonalization type changed {
502 // Create orthogonalization manager. This requires that the OutputManager (printer_) already be initialized.
503 Teuchos::RCP<const Teuchos::ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType_);
504 ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, Teuchos::null, printer_, label_, orthoParams);
505 }
506
507 // Convenience typedefs
508 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
509 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
510
511 if (impConvTest_.is_null()) {
512 impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
513 impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_), Belos::TwoNorm);
514 impConvTest_->setShowMaxResNormOnly( true );
515 }
516
517 if (expConvTest_.is_null()) {
518 expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
519 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
520 expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_), Belos::TwoNorm);
521 expConvTest_->setShowMaxResNormOnly( true );
522 }
523
524 if (convTest_.is_null()) {
525 convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ, impConvTest_, expConvTest_));
526 }
527
528 sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR, maxIterTest_, convTest_));
529
530 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
531 outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
532 Passed+Failed+Undefined);
533
534 */
535
536 }
537
538 // This method requires the solver manager to return a string that describes itself.
539 template<class ScalarType, class MV, class OP>
541 std::ostringstream oss;
542 oss << "Belos::BlockGCRODRSolMgr<" << SCT::name() << ", ...>";
543 oss << "{";
544 oss << "Ortho Type='"<<orthoType_ ;
545 oss << ", Num Blocks=" <<numBlocks_;
546 oss << ", Block Size=" <<blockSize_;
547 oss << ", Num Recycle Blocks=" << recycledBlocks_;
548 oss << ", Max Restarts=" << maxRestarts_;
549 oss << "}";
550 return oss.str();
551 }
552
553 template<class ScalarType, class MV, class OP>
554 Teuchos::RCP<const Teuchos::ParameterList>
556 using Teuchos::ParameterList;
557 using Teuchos::parameterList;
558 using Teuchos::RCP;
559 using Teuchos::rcpFromRef;
560
561 if (defaultParams_.is_null()) {
562 RCP<ParameterList> pl = parameterList ();
563
564 const MagnitudeType convTol = SMT::squareroot (SCT::magnitude (SCT::eps()));
565 const int maxRestarts = 1000;
566 const int maxIters = 5000;
567 const int blockSize = 2;
568 const int numBlocks = 100;
569 const int numRecycledBlocks = 25;
571 const Belos::OutputType outputStyle = Belos::General;
572 const int outputFreq = 1;
573 RCP<std::ostream> outputStream = rcpFromRef (std::cout);
574 const std::string impResScale ("Norm of Preconditioned Initial Residual");
575 const std::string expResScale ("Norm of Initial Residual");
576 const std::string timerLabel ("Belos");
577 const std::string orthoType ("ICGS");
578 RCP<const ParameterList> orthoParams = orthoFactory_.getDefaultParameters (orthoType);
579 //const MagnitudeType orthoKappa = SCT::magnitude (SCT::eps());
580 //
581 // mfh 31 Dec 2011: Negative values of orthoKappa are ignored.
582 // Setting this to a negative value by default ensures that
583 // this parameter is only _not_ ignored if the user
584 // specifically sets a valid value.
585 const MagnitudeType orthoKappa = -SMT::one();
586
587 // Set all the valid parameters and their default values.
588 pl->set ("Convergence Tolerance", convTol,
589 "The tolerance that the solver needs to achieve in order for "
590 "the linear system(s) to be declared converged. The meaning "
591 "of this tolerance depends on the convergence test details.");
592 pl->set("Maximum Restarts", maxRestarts,
593 "The maximum number of restart cycles (not counting the first) "
594 "allowed for each set of right-hand sides solved.");
595 pl->set("Maximum Iterations", maxIters,
596 "The maximum number of iterations allowed for each set of "
597 "right-hand sides solved.");
598 pl->set("Block Size", blockSize,
599 "How many linear systems to solve at once.");
600 pl->set("Num Blocks", numBlocks,
601 "The maximum number of blocks allowed in the Krylov subspace "
602 "for each set of right-hand sides solved. This includes the "
603 "initial block. Total Krylov basis storage, not counting the "
604 "recycled subspace, is \"Num Blocks\" times \"Block Size\".");
605 pl->set("Num Recycled Blocks", numRecycledBlocks,
606 "The maximum number of vectors in the recycled subspace." );
607 pl->set("Verbosity", verbosity,
608 "What type(s) of solver information should be written "
609 "to the output stream.");
610 pl->set("Output Style", outputStyle,
611 "What style is used for the solver information to write "
612 "to the output stream.");
613 pl->set("Output Frequency", outputFreq,
614 "How often convergence information should be written "
615 "to the output stream.");
616 pl->set("Output Stream", outputStream,
617 "A reference-counted pointer to the output stream where all "
618 "solver output is sent.");
619 pl->set("Implicit Residual Scaling", impResScale,
620 "The type of scaling used in the implicit residual convergence test.");
621 pl->set("Explicit Residual Scaling", expResScale,
622 "The type of scaling used in the explicit residual convergence test.");
623 pl->set("Timer Label", timerLabel,
624 "The string to use as a prefix for the timer labels.");
625 pl->set("Orthogonalization", orthoType,
626 "The orthogonalization method to use. Valid options: " +
627 orthoFactory_.validNamesString());
628 pl->set ("Orthogonalization Parameters", *orthoParams,
629 "Sublist of parameters specific to the orthogonalization method to use.");
630 pl->set("Orthogonalization Constant", orthoKappa,
631 "When using DGKS orthogonalization: the \"depTol\" constant, used "
632 "to determine whether another step of classical Gram-Schmidt is "
633 "necessary. Otherwise ignored. Nonpositive values are ignored.");
634 defaultParams_ = pl;
635 }
636 return defaultParams_;
637 }
638
639 template<class ScalarType, class MV, class OP>
640 void
642 setParameters (const Teuchos::RCP<Teuchos::ParameterList> &params) {
643 using Teuchos::isParameterType;
644 using Teuchos::getParameter;
645 using Teuchos::null;
646 using Teuchos::ParameterList;
647 using Teuchos::parameterList;
648 using Teuchos::RCP;
649 using Teuchos::rcp;
650 using Teuchos::rcp_dynamic_cast;
651 using Teuchos::rcpFromRef;
652 using Teuchos::Exceptions::InvalidParameter;
653 using Teuchos::Exceptions::InvalidParameterName;
654 using Teuchos::Exceptions::InvalidParameterType;
655
656 // The default parameter list contains all parameters that BlockGCRODRSolMgr understands, and none that it doesn't
657 RCP<const ParameterList> defaultParams = getValidParameters();
658
659 if (params.is_null()) {
660 // Users really should supply a non-null input ParameterList,
661 // so that we can fill it in. However, if they really did give
662 // us a null list, be nice and use default parameters. In this
663 // case, we don't have to do any validation.
664 params_ = parameterList (*defaultParams);
665 }
666 else {
667 // Validate the user's parameter list and fill in any missing parameters with default values.
668 //
669 // Currently, validation is quite strict. The following line
670 // will throw an exception for misspelled or extra parameters.
671 // However, it will _not_ throw an exception if there are
672 // missing parameters: these will be filled in with their
673 // default values. There is additional discussion of other
674 // validation strategies in the comments of this function for
675 // Belos::GCRODRSolMgr.
676 params->validateParametersAndSetDefaults (*defaultParams);
677 // No side effects on the solver until after validation passes.
678 params_ = params;
679 }
680
681 //
682 // Validate and set parameters relating to the Krylov subspace
683 // dimensions and the number of blocks.
684 //
685 // We try to read and validate all these parameters' values
686 // before setting them in the solver. This is because the
687 // validity of each may depend on the others' values. Our goal
688 // is to avoid incorrect side effects, so we don't set the values
689 // in the solver until we know they are right.
690 //
691 {
692 const int maxRestarts = params_->get<int> ("Maximum Restarts");
693 TEUCHOS_TEST_FOR_EXCEPTION(maxRestarts <= 0, std::invalid_argument,
694 "Belos::BlockGCRODRSolMgr: The \"Maximum Restarts\" parameter "
695 "must be nonnegative, but you specified a negative value of "
696 << maxRestarts << ".");
697
698 const int maxIters = params_->get<int> ("Maximum Iterations");
699 TEUCHOS_TEST_FOR_EXCEPTION(maxIters <= 0, std::invalid_argument,
700 "Belos::BlockGCRODRSolMgr: The \"Maximum Iterations\" parameter "
701 "must be positive, but you specified a nonpositive value of "
702 << maxIters << ".");
703
704 const int numBlocks = params_->get<int> ("Num Blocks");
705 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument,
706 "Belos::BlockGCRODRSolMgr: The \"Num Blocks\" parameter must be "
707 "positive, but you specified a nonpositive value of " << numBlocks
708 << ".");
709
710 const int blockSize = params_->get<int> ("Block Size");
711 TEUCHOS_TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument,
712 "Belos::BlockGCRODRSolMgr: The \"Block Size\" parameter must be "
713 "positive, but you specified a nonpositive value of " << blockSize
714 << ".");
715
716 const int recycledBlocks = params_->get<int> ("Num Recycled Blocks");
717 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks <= 0, std::invalid_argument,
718 "Belos::BlockGCRODRSolMgr: The \"Num Recycled Blocks\" parameter must "
719 "be positive, but you specified a nonpositive value of "
720 << recycledBlocks << ".");
721 TEUCHOS_TEST_FOR_EXCEPTION(recycledBlocks >= numBlocks,
722 std::invalid_argument, "Belos::BlockGCRODRSolMgr: The \"Num Recycled "
723 "Blocks\" parameter must be less than the \"Num Blocks\" parameter, "
724 "but you specified \"Num Recycled Blocks\" = " << recycledBlocks
725 << " and \"Num Blocks\" = " << numBlocks << ".");
726
727 // Now that we've validated the various dimension-related
728 // parameters, we can set them in the solvers.
729 maxRestarts_ = maxRestarts;
730 maxIters_ = maxIters;
731 numBlocks_ = numBlocks;
732 blockSize_ = blockSize;
733 recycledBlocks_ = recycledBlocks;
734 }
735
736 // Check to see if the timer label changed. If it did, update it in
737 // the parameter list, and create a new timer with that label (if
738 // Belos was compiled with timers enabled).
739 {
740 std::string tempLabel = params_->get<std::string> ("Timer Label");
741 const bool labelChanged = (tempLabel != label_);
742
743#ifdef BELOS_TEUCHOS_TIME_MONITOR
744 std::string solveLabel = label_ + ": BlockGCRODRSolMgr total solve time";
745 if (timerSolve_.is_null()) {
746 // We haven't created a timer yet.
747 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
748 } else if (labelChanged) {
749 // We created a timer before with a different label. In that
750 // case, clear the old timer and create a new timer with the
751 // new label. Clearing old timers prevents them from piling
752 // up, since they are stored statically, possibly without
753 // checking for duplicates.
754 Teuchos::TimeMonitor::clearCounter (solveLabel);
755 timerSolve_ = Teuchos::TimeMonitor::getNewCounter (solveLabel);
756 }
757#endif // BELOS_TEUCHOS_TIME_MONITOR
758 }
759
760 // Check for a change in verbosity level. Just in case, we allow
761 // the type of this parameter to be either int or MsgType (an
762 // enum).
763 if (params_->isParameter ("Verbosity")) {
764 if (isParameterType<int> (*params_, "Verbosity")) {
765 verbosity_ = params_->get<int> ("Verbosity");
766 }
767 else {
768 verbosity_ = (int) getParameter<MsgType> (*params_, "Verbosity");
769 }
770 }
771
772 // Check for a change in output style. Just in case, we allow
773 // the type of this parameter to be either int or OutputType (an
774 // enum).
775 if (params_->isParameter ("Output Style")) {
776 if (isParameterType<int> (*params_, "Output Style")) {
777 outputStyle_ = params_->get<int> ("Output Style");
778 }
779 else {
780 outputStyle_ = (int) getParameter<OutputType> (*params_, "Output Style");
781 }
782
783 // Currently, the only way we can update the output style is to
784 // create a new output test. This is because we must first
785 // instantiate a new StatusTestOutputFactory with the new
786 // style, and then use the factory to make a new output test.
787 // Thus, we set outputTest_ to null for now, so we remember
788 // later to create a new one.
789 outputTest_ = null;
790 }
791
792 // Get the output stream for the output manager.
793 //
794 // It has been pointed out (mfh 28 Feb 2011 in GCRODRSolMgr code)
795 // that including an RCP<std::ostream> parameter makes it
796 // impossible to serialize the parameter list. If you serialize
797 // the list and read it back in from the serialized
798 // representation, you might not even get a valid
799 // RCP<std::ostream>, let alone the same output stream that you
800 // serialized.
801 //
802 // However, existing Belos users depend on setting the output
803 // stream in the parameter list. We retain this behavior for
804 // backwards compatibility.
805 //
806 // In case the output stream can't be read back in, we default to
807 // stdout (std::cout), just to ensure reasonable behavior.
808 if (params_->isParameter ("Output Stream")) {
809 try {
810 outputStream_ = getParameter<RCP<std::ostream> > (*params_, "Output Stream");
811 }
812 catch (InvalidParameter&) {
813 outputStream_ = rcpFromRef (std::cout);
814 }
815 // We assume that a null output stream indicates that the user
816 // doesn't want to print anything, so we replace it with a
817 // "black hole" stream that prints nothing sent to it. (We
818 // can't use outputStream_ = Teuchos::null, since the output
819 // manager assumes that operator<< works on its output stream.)
820 if (outputStream_.is_null()) {
821 outputStream_ = rcp (new Teuchos::oblackholestream);
822 }
823 }
824
825 outputFreq_ = params_->get<int> ("Output Frequency");
826 if (verbosity_ & Belos::StatusTestDetails) {
827 // Update parameter in our output status test.
828 if (! outputTest_.is_null()) {
829 outputTest_->setOutputFrequency (outputFreq_);
830 }
831 }
832
833 // Create output manager if we need to, using the verbosity level
834 // and output stream that we fetched above. We do this here because
835 // instantiating an OrthoManager using OrthoManagerFactory requires
836 // a valid OutputManager.
837 if (printer_.is_null()) {
838 printer_ = rcp (new OutputManager<ScalarType> (verbosity_, outputStream_));
839 } else {
840 printer_->setVerbosity (verbosity_);
841 printer_->setOStream (outputStream_);
842 }
843
844 // Get the orthogonalization manager name ("Orthogonalization").
845 //
846 // Getting default values for the orthogonalization manager
847 // parameters ("Orthogonalization Parameters") requires knowing the
848 // orthogonalization manager name. Save it for later, and also
849 // record whether it's different than before.
850
851 //bool changedOrthoType = false; // silence warning about not being used
852 if (params_->isParameter ("Orthogonalization")) {
853 const std::string& tempOrthoType =
854 params_->get<std::string> ("Orthogonalization");
855 // Ensure that the specified orthogonalization type is valid.
856 if (! orthoFactory_.isValidName (tempOrthoType)) {
857 std::ostringstream os;
858 os << "Belos::BlockGCRODRSolMgr: Invalid orthogonalization name \""
859 << tempOrthoType << "\". The following are valid options "
860 << "for the \"Orthogonalization\" name parameter: ";
861 orthoFactory_.printValidNames (os);
862 TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument, os.str());
863 }
864 if (tempOrthoType != orthoType_) {
865 //changedOrthoType = true;
866 orthoType_ = tempOrthoType;
867 }
868 }
869
870 // Get any parameters for the orthogonalization
871 // ("Orthogonalization Parameters"). If not supplied, the
872 // orthogonalization manager factory will supply default values.
873 //
874 // NOTE (mfh 12 Jan 2011, summer 2011, 31 Dec 2011) For backwards
875 // compatibility with other Belos GMRES-type solvers, if params
876 // has an "Orthogonalization Constant" parameter and the DGKS
877 // orthogonalization manager is to be used, the value of this
878 // parameter will override DGKS's "depTol" parameter.
879 //
880 // Users must supply the orthogonalization manager parameters as
881 // a sublist (supplying it as an RCP<ParameterList> would make
882 // the resulting parameter list not serializable).
883 RCP<ParameterList> orthoParams = sublist (params, "Orthogonalization Parameters", true);
884 TEUCHOS_TEST_FOR_EXCEPTION(orthoParams.is_null(), std::logic_error,
885 "Failed to get orthogonalization parameters. "
886 "Please report this bug to the Belos developers.");
887
888 // Check if the desired orthogonalization method changed, or if
889 // the orthogonalization manager has not yet been instantiated.
890 // If either is the case, instantiate a new MatOrthoManager
891 // subclass instance corresponding to the desired
892 // orthogonalization method. We've already fetched the
893 // orthogonalization method name (orthoType_) and its parameters
894 // (orthoParams) above.
895 //
896 // As suggested (by mfh 12 Jan 2011 in Belos::GCRODRSolMgr) In
897 // order to ensure parameter changes get propagated to the
898 // orthomanager we simply reinstantiate the OrthoManager every
899 // time, whether or not the orthogonalization method name or
900 // parameters have changed. This is not efficient for some
901 // orthogonalization managers that keep a lot of internal state.
902 // A more general way to fix this inefficiency would be to
903 //
904 // 1. make each orthogonalization implement
905 // Teuchos::ParameterListAcceptor, and
906 //
907 // 2. make each orthogonalization implement a way to set the
908 // OutputManager instance and timer label.
909 //
910 // Create orthogonalization manager. This requires that the
911 // OutputManager (printer_) already be initialized.
912 ortho_ = orthoFactory_.makeMatOrthoManager (orthoType_, null, printer_,
913 label_, orthoParams);
914
915 // The DGKS orthogonalization accepts a "Orthogonalization
916 // Constant" parameter (also called kappa in the code, but not in
917 // the parameter list). If its value is provided in the given
918 // parameter list and its value is positive, use it. Ignore
919 // negative values.
920 //
921 // The "Orthogonalization Constant" parameter is retained only
922 // for backwards compatibility.
923
924 //bool gotValidOrthoKappa = false; // silence warning about not being used
925 if (params_->isParameter ("Orthogonalization Constant")) {
926 const MagnitudeType orthoKappa =
927 params_->get<MagnitudeType> ("Orthogonalization Constant");
928 if (orthoKappa > 0) {
929 orthoKappa_ = orthoKappa;
930 //gotValidOrthoKappa = true;
931 // Only DGKS currently accepts this parameter.
932 if (orthoType_ == "DGKS" && ! ortho_.is_null()) {
933 typedef DGKSOrthoManager<ScalarType, MV, OP> ortho_man_type;
934 // This cast should always succeed; it's a bug otherwise.
935 // (If the cast fails, then orthoType_ doesn't correspond
936 // to the OrthoManager subclass instance that we think we
937 // have, so we initialized the wrong subclass somehow.)
938 rcp_dynamic_cast<ortho_man_type>(ortho_)->setDepTol (orthoKappa_);
939 }
940 }
941 }
942
943 // Convergence
944 typedef Belos::StatusTestCombo<ScalarType,MV,OP> StatusTestCombo_t;
945 typedef Belos::StatusTestGenResNorm<ScalarType,MV,OP> StatusTestResNorm_t;
946
947 // Check for convergence tolerance
948 convTol_ = params_->get<MagnitudeType> ("Convergence Tolerance");
949 if (! impConvTest_.is_null()) {
950 impConvTest_->setTolerance (convTol_);
951 }
952 if (! expConvTest_.is_null()) {
953 expConvTest_->setTolerance (convTol_);
954 }
955
956 // Check for a change in scaling, if so we need to build new residual tests.
957 if (params_->isParameter ("Implicit Residual Scaling")) {
958 std::string tempImpResScale =
959 getParameter<std::string> (*params_, "Implicit Residual Scaling");
960
961 // Only update the scaling if it's different.
962 if (impResScale_ != tempImpResScale) {
963 ScaleType impResScaleType = convertStringToScaleType (tempImpResScale);
964 impResScale_ = tempImpResScale;
965
966 if (! impConvTest_.is_null()) {
967 try {
968 impConvTest_->defineScaleForm (impResScaleType, Belos::TwoNorm);
969 }
970 catch (StatusTestError&) {
971 // Delete the convergence test so it gets constructed again.
972 impConvTest_ = null;
973 convTest_ = null;
974 }
975 }
976 }
977 }
978
979 if (params_->isParameter("Explicit Residual Scaling")) {
980 std::string tempExpResScale =
981 getParameter<std::string> (*params_, "Explicit Residual Scaling");
982
983 // Only update the scaling if it's different.
984 if (expResScale_ != tempExpResScale) {
985 ScaleType expResScaleType = convertStringToScaleType (tempExpResScale);
986 expResScale_ = tempExpResScale;
987
988 if (! expConvTest_.is_null()) {
989 try {
990 expConvTest_->defineScaleForm (expResScaleType, Belos::TwoNorm);
991 }
992 catch (StatusTestError&) {
993 // Delete the convergence test so it gets constructed again.
994 expConvTest_ = null;
995 convTest_ = null;
996 }
997 }
998 }
999 }
1000
1001 //
1002 // Create iteration stopping criteria ("status tests") if we need
1003 // to, by combining three different stopping criteria.
1004 //
1005 // First, construct maximum-number-of-iterations stopping criterion.
1006 if (maxIterTest_.is_null()) {
1007 maxIterTest_ = rcp (new StatusTestMaxIters<ScalarType,MV,OP> (maxIters_));
1008 } else {
1009 maxIterTest_->setMaxIters (maxIters_);
1010 }
1011
1012 // Implicit residual test, using the native residual to determine if
1013 // convergence was achieved.
1014 if (impConvTest_.is_null()) {
1015 impConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1016 impConvTest_->defineScaleForm (convertStringToScaleType (impResScale_),
1018 }
1019
1020 // Explicit residual test once the native residual is below the tolerance
1021 if (expConvTest_.is_null()) {
1022 expConvTest_ = rcp (new StatusTestResNorm_t (convTol_));
1023 expConvTest_->defineResForm (StatusTestResNorm_t::Explicit, Belos::TwoNorm);
1024 expConvTest_->defineScaleForm (convertStringToScaleType (expResScale_),
1026 }
1027 // Convergence test first tests the implicit residual, then the
1028 // explicit residual if the implicit residual test passes.
1029 if (convTest_.is_null()) {
1030 convTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::SEQ,
1031 impConvTest_,
1032 expConvTest_));
1033 }
1034 // Construct the complete iteration stopping criterion:
1035 //
1036 // "Stop iterating if the maximum number of iterations has been
1037 // reached, or if the convergence test passes."
1038 sTest_ = rcp (new StatusTestCombo_t (StatusTestCombo_t::OR,
1039 maxIterTest_, convTest_));
1040 // Create the status test output class.
1041 // This class manages and formats the output from the status test.
1042 StatusTestOutputFactory<ScalarType,MV,OP> stoFactory (outputStyle_);
1043 outputTest_ = stoFactory.create (printer_, sTest_, outputFreq_,
1045
1046 // Set the solver string for the output test
1047 std::string solverDesc = "Block GCRODR ";
1048 outputTest_->setSolverDesc (solverDesc);
1049
1050 // Inform the solver manager that the current parameters were set.
1051 isSet_ = true;
1052 }
1053
1054 // initializeStateStorage.
1055 template<class ScalarType, class MV, class OP>
1056 void
1058 {
1059
1060 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1061
1062 // Check if there is any multivector to clone from.
1063 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1064
1065 //The Dimension of the Krylov Subspace
1066 int KrylSpaDim = (numBlocks_ - 1) * blockSize_;
1067
1068 //Number of columns in [U_ V_(:,1:KrylSpaDim)]
1069 int augSpaDim = KrylSpaDim + recycledBlocks_ + 1;// + 1 is for possible extra recycle vector
1070
1071 //Number of columns in [C_ V_]
1072 int augSpaImgDim = KrylSpaDim + blockSize_ + recycledBlocks_+1;
1073
1074 //TEMPORARY SKELETON DEFINITION OF THIS FUNCTION TO GET THINGS WORKING
1075 //NOT EVERYTHING IS INITIALIZE CORRECTLY YET.
1076
1077 //INITIALIZE RECYCLE SPACE VARIABLES HERE
1078
1079 //WE DO NOT ALLOCATE V HERE IN THIS SOLVERMANAGER. If no recycle space exists, then block_gmres_iter
1080 //will allocated V for us. If a recycle space already exists, then we will allocate V after updating the
1081 //recycle space for the current problem.
1082 // If the block Krylov subspace has not been initialized before, generate it using the RHS from lp_.
1083 /*if (V_ == Teuchos::null) {
1084 V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1085 }
1086 else{
1087 // Generate V_ by cloning itself ONLY if more space is needed.
1088 if (MVT::GetNumberVecs(*V_) < numBlocks_+1) {
1089 Teuchos::RCP<const MV> tmp = V_;
1090 V_ = MVT::Clone( *tmp, numBlocks_+1 );
1091 }
1092 }*/
1093
1094 //INTITIALIZE SPACE FOR THE NEWLY COMPUTED RECYCLE SPACE VARIABLES HERE
1095
1096 if (U_ == Teuchos::null) {
1097 U_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1098 }
1099 else {
1100 // Generate U_ by cloning itself ONLY if more space is needed.
1101 if (MVT::GetNumberVecs(*U_) < recycledBlocks_+1) {
1102 Teuchos::RCP<const MV> tmp = U_;
1103 U_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1104 }
1105 }
1106
1107 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1108 if (C_ == Teuchos::null) {
1109 C_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1110 }
1111 else {
1112 // Generate C_ by cloning itself ONLY if more space is needed.
1113 if (MVT::GetNumberVecs(*C_) < recycledBlocks_+1) {
1114 Teuchos::RCP<const MV> tmp = C_;
1115 C_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1116 }
1117 }
1118
1119 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1120 if (U1_ == Teuchos::null) {
1121 U1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1122 }
1123 else {
1124 // Generate U1_ by cloning itself ONLY if more space is needed.
1125 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1126 Teuchos::RCP<const MV> tmp = U1_;
1127 U1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1128 }
1129 }
1130
1131 // If the subspace has not been initialized before, generate it using the RHS from lp_.
1132 if (C1_ == Teuchos::null) {
1133 C1_ = MVT::Clone( *rhsMV, recycledBlocks_+1 );
1134 }
1135 else {
1136 // Generate C1_ by cloning itself ONLY if more space is needed.
1137 if (MVT::GetNumberVecs(*U1_) < recycledBlocks_+1) {
1138 Teuchos::RCP<const MV> tmp = C1_;
1139 C1_ = MVT::Clone( *tmp, recycledBlocks_+1 );
1140 }
1141 }
1142
1143 // Generate R_ only if it doesn't exist
1144 if (R_ == Teuchos::null){
1145 R_ = MVT::Clone( *rhsMV, blockSize_ );
1146 }
1147
1148 //INITIALIZE SOME WORK VARIABLES
1149
1150 // Generate G_ only if it doesn't exist, otherwise resize it.
1151 if (G_ == Teuchos::null){
1152 G_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
1153 }
1154 else{
1155 if ( (G_->numRows() != augSpaImgDim) || (G_->numCols() != augSpaDim) )
1156 {
1157 G_->reshape( augSpaImgDim, augSpaDim );
1158 }
1159 G_->putScalar(zero);
1160 }
1161
1162 // Generate H_ only if it doesn't exist by pointing it to a view of G_.
1163 if (H_ == Teuchos::null){
1164 H_ = Teuchos::rcp (new SDM ( Teuchos::View, *G_, KrylSpaDim + blockSize_, KrylSpaDim, recycledBlocks_+1 ,recycledBlocks_+1 ) );
1165 }
1166
1167 // Generate F_ only if it doesn't exist, otherwise resize it.
1168 if (F_ == Teuchos::null){
1169 F_ = Teuchos::rcp( new SDM( recycledBlocks_+1, recycledBlocks_+1 ) );
1170 }
1171 else {
1172 if ( (F_->numRows() != recycledBlocks_+1) || (F_->numCols() != recycledBlocks_+1) ){
1173 F_->reshape( recycledBlocks_+1, recycledBlocks_+1 );
1174 }
1175 }
1176 F_->putScalar(zero);
1177
1178 // Generate PP_ only if it doesn't exist, otherwise resize it.
1179 if (PP_ == Teuchos::null){
1180 PP_ = Teuchos::rcp( new SDM( augSpaImgDim, recycledBlocks_+1 ) );
1181 }
1182 else{
1183 if ( (PP_->numRows() != augSpaImgDim) || (PP_->numCols() != recycledBlocks_+1) ){
1184 PP_->reshape( augSpaImgDim, recycledBlocks_+1 );
1185 }
1186 }
1187
1188 // Generate HP_ only if it doesn't exist, otherwise resize it.
1189 if (HP_ == Teuchos::null)
1190 HP_ = Teuchos::rcp( new SDM( augSpaImgDim, augSpaDim ) );
1191 else{
1192 if ( (HP_->numRows() != augSpaImgDim) || (HP_->numCols() != augSpaDim) ){
1193 HP_->reshape( augSpaImgDim, augSpaDim );
1194 }
1195 }
1196
1197 // Size of tau_ will change during computation, so just be sure it starts with appropriate size
1198 tau_.resize(recycledBlocks_+1);
1199
1200 // Size of work_ will change during computation, so just be sure it starts with appropriate size
1201 work_.resize(recycledBlocks_+1);
1202
1203 // Size of ipiv_ will change during computation, so just be sure it starts with appropriate size
1204 ipiv_.resize(recycledBlocks_+1);
1205
1206 }
1207
1208template<class ScalarType, class MV, class OP>
1210
1211 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1212 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1213
1214 int p = block_gmres_iter->getState().curDim; //Dimension of the Krylov space generated
1215 std::vector<int> index(keff);//we use this to index certain columns of U, C, and V to
1216 //get views into pieces of these matrices.
1217
1218 //GET CORRECT PIECE OF MATRIX H APPROPRIATE TO SIZE OF KRYLOV SUBSPACE
1219 SDM HH(Teuchos::Copy, *H_, p+blockSize_, p);
1220 if(recycledBlocks_ >= p + blockSize_){//keep whole block Krylov subspace
1221 //IF THIS HAS HAPPENED, THIS MEANS WE CONVERGED DURING THIS CYCLE
1222 //THEREFORE, WE DO NOT CONSTRUCT C = A*U;
1223 index.resize(p);
1224 for (int ii=0; ii < p; ++ii) index[ii] = ii;
1225 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1226 MVT::SetBlock(*V_, index, *Utmp);
1227 keff = p;
1228 }
1229 else{ //use a subspace selection method to get recycle space
1230 int info = 0;
1231 Teuchos::RCP<SDM > PPtmp = Teuchos::rcp (new SDM ( Teuchos::View, *PP_, p, recycledBlocks_+1 ) );
1232 if(recycleMethod_ == "harmvecs"){
1233 keff = getHarmonicVecsKryl(p, HH, *PPtmp);
1234 printer_->stream(Debug) << "keff = " << keff << std::endl;
1235 }
1236// Hereafter, only keff columns of PP are needed
1237PPtmp = Teuchos::rcp (new SDM ( Teuchos::View, *PP_, p, keff ) );
1238// Now get views into C, U, V
1239index.resize(keff);
1240for (int ii=0; ii<keff; ++ii) index[ii] = ii;
1241Teuchos::RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1242Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1243Teuchos::RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1244index.resize(p);
1245for (int ii=0; ii < p; ++ii) index[ii] = ii;
1246Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1247
1248// Form U (the subspace to recycle)
1249// U = newstate.V(:,1:p) * PP;
1250MVT::MvTimesMatAddMv( one, *Vtmp, *PPtmp, zero, *U1tmp );
1251
1252// Step #1: Form HP = H*P
1253SDM HPtmp( Teuchos::View, *HP_, p+blockSize_, keff );
1254HPtmp.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *H_, *PPtmp, zero );
1255// Step #1.5: Perform workspace size query for QR factorization of HP
1256int lwork = -1;
1257tau_.resize(keff);
1258lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1259TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1260
1261// Step #2: Compute QR factorization of HP
1262 lwork = static_cast<int> (Teuchos::ScalarTraits<ScalarType>::magnitude (work_[0]));
1263work_.resize(lwork);
1264lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1265TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1266
1267// Step #3: Explicitly construct Q and R factors
1268// The upper triangular part of HP is copied into R and HP becomes Q.
1269SDM Rtmp( Teuchos::View, *F_, keff, keff );
1270for(int ii=0;ii<keff;ii++) { for(int jj=ii;jj<keff;jj++) Rtmp(ii,jj) = HPtmp(ii,jj); }
1271//lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1272lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1273TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1274 // Now we have [Q,R] = qr(H*P)
1275
1276 // Now compute C = V(:,1:p+blockSize_) * Q
1277 index.resize( p+blockSize_ );
1278 for (int ii=0; ii < (p+blockSize_); ++ii) { index[ii] = ii; }
1279 Vtmp = MVT::CloneView( *V_, index ); // need new view into V (p+blockSize_ vectors now; needed p above)
1280 MVT::MvTimesMatAddMv( one, *Vtmp, HPtmp, zero, *Ctmp );
1281
1282 // Finally, compute U = U*R^{-1}.
1283 // This unfortuntely requires me to form R^{-1} explicitly and execute U = U * R^{-1}, as
1284 // backsolve capabilities don't exist in the Belos::MultiVec class
1285
1286// Step #1: First, compute LU factorization of R
1287ipiv_.resize(Rtmp.numRows());
1288lapack.GETRF(Rtmp.numRows(),Rtmp.numCols(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&info);
1289TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1290// Step #2: Form inv(R)
1291lwork = Rtmp.numRows();
1292work_.resize(lwork);
1293lapack.GETRI(Rtmp.numRows(),Rtmp.values(),Rtmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1294//TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::GCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1295// Step #3: Let U = U * R^{-1}
1296MVT::MvTimesMatAddMv( one, *U1tmp, Rtmp, zero, *Utmp );
1297
1298} //end else from if(recycledBlocks_ >= p + 1)
1299return;
1300} // end buildRecycleSpaceKryl defnition
1301
1302template<class ScalarType, class MV, class OP>
1304 const MagnitudeType one = Teuchos::ScalarTraits<ScalarType>::one();
1305 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1306
1307 std::vector<MagnitudeType> d(keff);
1308 std::vector<ScalarType> dscalar(keff);
1309 std::vector<int> index(numBlocks_+1);
1310
1311 // Get the state
1312 BlockGCRODRIterState<ScalarType,MV> oldState = block_gcrodr_iter->getState();
1313 int p = oldState.curDim;
1314
1315 // insufficient new information to update recycle space
1316 if (p<1) return;
1317
1318 if(recycledBlocks_ >= keff + p){ //we add new Krylov vectors to existing recycle space
1319 // If this has happened, we have converged during this cycle. Therefore, do not construct C = A*C;
1320 index.resize(p);
1321 for (int ii=0; ii < p; ++ii) { index[ii] = keff+ii; } //get a view after current reycle vectors
1322 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1323 for (int ii=0; ii < p; ++ii) { index[ii] =ii; }
1324 MVT::SetBlock(*V_, index, *Utmp);
1325 keff += p;
1326 }
1327
1328 // Take the norm of the recycled vectors
1329 {
1330 index.resize(keff);
1331 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1332 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1333 d.resize(keff);
1334 dscalar.resize(keff);
1335 MVT::MvNorm( *Utmp, d );
1336 for (int i=0; i<keff; ++i) {
1337 d[i] = one / d[i];
1338 dscalar[i] = (ScalarType)d[i];
1339 }
1340 MVT::MvScale( *Utmp, dscalar );
1341 }
1342
1343 // Get view into current "full" upper Hessnburg matrix
1344 // note that p describes the dimension of the iter+1 block Krylov space so we have to adjust how we use it to index Gtmp
1345 Teuchos::RCP<SDM> Gtmp = Teuchos::rcp( new SDM( Teuchos::View, *G_, p+keff, p+keff-blockSize_ ) );
1346
1347 // Insert D into the leading keff x keff block of G
1348 for (int i=0; i<keff; ++i)
1349 (*Gtmp)(i,i) = d[i];
1350
1351 // Compute the harmoic Ritz pairs for the generalized eigenproblem
1352 // getHarmonicVecsKryl assumes PP has recycledBlocks_+1 columns available
1353 // See previous block of comments for why we subtract p-blockSize_
1354 int keff_new;
1355 {
1356 SDM PPtmp( Teuchos::View, *PP_, p+keff-blockSize_, recycledBlocks_+1 );
1357 keff_new = getHarmonicVecsAugKryl( keff, p-blockSize_, *Gtmp, oldState.V, PPtmp );
1358 }
1359
1360 // Code to form new U, C
1361 // U = [U V(:,1:p)] * P; (in two steps)
1362
1363 // U(:,1:keff) = matmul(U(:,1:keff_old),PP(1:keff_old,1:keff)) (step 1)
1364 Teuchos::RCP<MV> U1tmp;
1365 {
1366 index.resize( keff );
1367 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1368 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1369 index.resize( keff_new );
1370 for (int ii=0; ii<keff_new; ++ii) { index[ii] = ii; }
1371 U1tmp = MVT::CloneViewNonConst( *U1_, index );
1372 SDM PPtmp( Teuchos::View, *PP_, keff, keff_new );
1373 MVT::MvTimesMatAddMv( one, *Utmp, PPtmp, zero, *U1tmp );
1374 }
1375
1376 // U(:,1:keff) = U(:,1:keff) + matmul(V(:,1:m-k),PP(keff_old+1:m-k+keff_old,1:keff)) (step 2)
1377 {
1378 index.resize(p-blockSize_);
1379 for (int ii=0; ii < p-blockSize_; ii++) { index[ii] = ii; }
1380 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1381 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_, keff_new, keff );
1382 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *U1tmp );
1383 }
1384
1385 // Form GP = G*P
1386 SDM HPtmp( Teuchos::View, *HP_, p+keff, keff_new );
1387 {
1388 SDM PPtmp( Teuchos::View, *PP_, p-blockSize_+keff, keff_new );
1389 HPtmp.multiply(Teuchos::NO_TRANS,Teuchos::NO_TRANS,one,*Gtmp,PPtmp,zero);
1390 }
1391
1392 // Workspace size query for QR factorization HP= QF (the worksize will be placed in work_[0])
1393 int info = 0, lwork = -1;
1394 tau_.resize(keff_new);
1395 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1396 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a workspace size.");
1397
1398 lwork = static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work_[0]));
1399 work_.resize(lwork);
1400 lapack.GEQRF(HPtmp.numRows(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1401 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GEQRF failed to compute a QR factorization.");
1402
1403 // Explicitly construct Q and F factors
1404 // NOTE: The upper triangular part of HP is copied into F and HP becomes Q.
1405 SDM Ftmp( Teuchos::View, *F_, keff_new, keff_new );
1406 for(int i=0;i<keff_new;i++) { for(int j=i;j<keff_new;j++) Ftmp(i,j) = HPtmp(i,j); }
1407 //lapack.ORGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1408 lapack.UNGQR(HPtmp.numRows(),HPtmp.numCols(),HPtmp.numCols(),HPtmp.values(),HPtmp.stride(),&tau_[0],&work_[0],lwork,&info);
1409 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _UNGQR failed to construct the Q factor.");
1410
1411 // Form orthonormalized C and adjust U accordingly so that C = A*U
1412 // C = [C V] * Q;
1413
1414 // C(:,1:keff) = matmul(C(:,1:keff_old),QQ(1:keff_old,1:keff))
1415 {
1416 Teuchos::RCP<MV> C1tmp;
1417 {
1418 index.resize(keff);
1419 for (int i=0; i < keff; i++) { index[i] = i; }
1420 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1421 index.resize(keff_new);
1422 for (int i=0; i < keff_new; i++) { index[i] = i; }
1423 C1tmp = MVT::CloneViewNonConst( *C1_, index );
1424 SDM PPtmp( Teuchos::View, *HP_, keff, keff_new );
1425 MVT::MvTimesMatAddMv( one, *Ctmp, PPtmp, zero, *C1tmp );
1426 }
1427 // Now compute C += V(:,1:p+1) * Q
1428 {
1429 index.resize( p );
1430 for (int i=0; i < p; ++i) { index[i] = i; }
1431 Teuchos::RCP<const MV> Vtmp = MVT::CloneView( *V_, index );
1432 SDM PPtmp( Teuchos::View, *HP_, p, keff_new, keff, 0 );
1433 MVT::MvTimesMatAddMv( one, *Vtmp, PPtmp, one, *C1tmp );
1434 }
1435 }
1436
1437 // C_ = C1_; (via a swap)
1438 std::swap(C_, C1_);
1439
1440 // Finally, compute U_ = U_*R^{-1}
1441 // First, compute LU factorization of R
1442 ipiv_.resize(Ftmp.numRows());
1443 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1444 TEUCHOS_TEST_FOR_EXCEPTION(info != 0,BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1445
1446 // Now, form inv(R)
1447 lwork = Ftmp.numRows();
1448 work_.resize(lwork);
1449 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1450 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to compute an LU factorization.");
1451
1452 {
1453 index.resize(keff_new);
1454 for (int i=0; i < keff_new; i++) { index[i] = i; }
1455 Teuchos::RCP<MV> Utmp = MVT::CloneViewNonConst( *U_, index );
1456 MVT::MvTimesMatAddMv( one, *U1tmp, Ftmp, zero, *Utmp );
1457 }
1458
1459 // Set the current number of recycled blocks and subspace
1460 // dimension with the Block GCRO-DR iteration.
1461 if (keff != keff_new) {
1462 keff = keff_new;
1463 block_gcrodr_iter->setSize( keff, numBlocks_ );
1464 // Important to zero this out before next cyle
1465 SDM b1( Teuchos::View, *G_, recycledBlocks_+2, 1, 0, recycledBlocks_ );
1466 b1.putScalar(zero);
1467 }
1468
1469} //end buildRecycleSpaceAugKryl definition
1470
1471template<class ScalarType, class MV, class OP>
1472int BlockGCRODRSolMgr<ScalarType,MV,OP>::getHarmonicVecsAugKryl(int keff, int m, const SDM& GG, const Teuchos::RCP<const MV>& VV, SDM& PP){
1473 int i, j;
1474 int m2 = GG.numCols();
1475 bool xtraVec = false;
1476 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1477 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1478 std::vector<int> index;
1479
1480 // Real and imaginary eigenvalue components
1481 std::vector<MagnitudeType> wr(m2), wi(m2);
1482
1483 // Magnitude of harmonic Ritz values
1484 std::vector<MagnitudeType> w(m2);
1485
1486 // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
1487 SDM vr(m2,m2,false);
1488
1489 // Sorted order of harmonic Ritz values
1490 std::vector<int> iperm(m2);
1491
1492 // Set flag indicating recycle space has been generated this solve
1493 builtRecycleSpace_ = true;
1494
1495 // Form matrices for generalized eigenproblem
1496
1497 // B = G' * G; Don't zero out matrix when constructing
1498 SDM B(m2,m2,false);
1499 B.multiply(Teuchos::TRANS,Teuchos::NO_TRANS,one,GG,GG,zero);
1500
1501 // A_tmp = | C'*U 0 |
1502 // | V_{m+1}'*U I |
1503 SDM A_tmp( keff+m+blockSize_, keff+m );
1504
1505
1506 // A_tmp(1:keff,1:keff) = C' * U;
1507 index.resize(keff);
1508 for (int i=0; i<keff; ++i) { index[i] = i; }
1509 Teuchos::RCP<const MV> Ctmp = MVT::CloneView( *C_, index );
1510 Teuchos::RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1511 SDM A11( Teuchos::View, A_tmp, keff, keff );
1512 MVT::MvTransMv( one, *Ctmp, *Utmp, A11 );
1513
1514 // A_tmp(keff+1:m-k+keff+1,1:keff) = V' * U;
1515 SDM A21( Teuchos::View, A_tmp, m+blockSize_, keff, keff );
1516 index.resize(m+blockSize_);
1517 for (i=0; i < m+blockSize_; i++) { index[i] = i; }
1518 Teuchos::RCP<const MV> Vp = MVT::CloneView( *VV, index );
1519 MVT::MvTransMv( one, *Vp, *Utmp, A21 );
1520
1521 // A_tmp(keff+1:m-k+keff,keff+1:m-k+keff) = eye(m-k);
1522 for( i=keff; i<keff+m; i++)
1523 A_tmp(i,i) = one;
1524
1525 // A = G' * A_tmp;
1526 SDM A( m2, A_tmp.numCols() );
1527 A.multiply( Teuchos::TRANS, Teuchos::NO_TRANS, one, GG, A_tmp, zero );
1528
1529 // Compute k smallest harmonic Ritz pairs
1530 // SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
1531 // ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
1532 // IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
1533 // RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
1534 // MLP: 'SCALING' in DGGEVX generates incorrect eigenvalues. Therefore, only permuting
1535 char balanc='P', jobvl='N', jobvr='V', sense='N';
1536 int ld = A.numRows();
1537 int lwork = 6*ld;
1538 int ldvl = ld, ldvr = ld;
1539 int info = 0,ilo = 0,ihi = 0;
1540 MagnitudeType abnrm = 0.0, bbnrm = 0.0;
1541 ScalarType *vl = 0; // This is never referenced by dggevx if jobvl == 'N'
1542 std::vector<ScalarType> beta(ld);
1543 std::vector<ScalarType> work(lwork);
1544 std::vector<MagnitudeType> rwork(lwork);
1545 std::vector<MagnitudeType> lscale(ld), rscale(ld);
1546 std::vector<MagnitudeType> rconde(ld), rcondv(ld);
1547 std::vector<int> iwork(ld+6);
1548 int *bwork = 0; // If sense == 'N', bwork is never referenced
1549 //lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1550 // &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1551 // &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &iwork[0], bwork, &info);
1552 lapack.GGEVX(balanc, jobvl, jobvr, sense, ld, A.values(), ld, B.values(), ld, &wr[0], &wi[0],
1553 &beta[0], vl, ldvl, vr.values(), ldvr, &ilo, &ihi, &lscale[0], &rscale[0],
1554 &abnrm, &bbnrm, &rconde[0], &rcondv[0], &work[0], lwork, &rwork[0],
1555 &iwork[0], bwork, &info);
1556 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GGEVX failed to compute eigensolutions.");
1557
1558 // Construct magnitude of each harmonic Ritz value
1559 // NOTE : Forming alpha/beta *should* be okay here, given assumptions on construction of matrix pencil above
1560 for( i=0; i<ld; i++ ) // Construct magnitude of each harmonic Ritz value
1561 w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] ) / std::abs( beta[i] );
1562
1563 this->sort(w,ld,iperm);
1564
1565 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1566
1567 // Select recycledBlocks_ smallest eigenvectors
1568 for( i=0; i<recycledBlocks_; i++ )
1569 for( j=0; j<ld; j++ )
1570 PP(j,i) = vr(j,iperm[ld-recycledBlocks_+i]);
1571
1572 if(scalarTypeIsComplex==false) {
1573
1574 // Determine exact size for PP (i.e., determine if we need to store an additional vector)
1575 if (wi[iperm[ld-recycledBlocks_]] != 0.0) {
1576 int countImag = 0;
1577 for ( i=ld-recycledBlocks_; i<ld; i++ )
1578 if (wi[iperm[i]] != 0.0) countImag++;
1579 // Check to see if this count is even or odd:
1580 if (countImag % 2) xtraVec = true;
1581 }
1582
1583 if (xtraVec) { // we need to store one more vector
1584 if (wi[iperm[ld-recycledBlocks_]] > 0.0) { // I picked the "real" component
1585 for( j=0; j<ld; j++ ) // so get the "imag" component
1586 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]+1);
1587 }
1588 else { // I picked the "imag" component
1589 for( j=0; j<ld; j++ ) // so get the "real" component
1590 PP(j,recycledBlocks_) = vr(j,iperm[ld-recycledBlocks_]-1);
1591 }
1592 }
1593
1594 }
1595
1596 // Return whether we needed to store an additional vector
1597 if (xtraVec)
1598 return recycledBlocks_+1;
1599 else
1600 return recycledBlocks_;
1601
1602} //end getHarmonicVecsAugKryl definition
1603
1604template<class ScalarType, class MV, class OP>
1606 bool xtraVec = false;
1607 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1608 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1609
1610 // Real and imaginary eigenvalue components
1611 std::vector<MagnitudeType> wr(m), wi(m);
1612
1613 // Real and imaginary (right) eigenvectors; Don't zero out matrix when constructing
1614 SDM vr(m,m,false);
1615
1616 // Magnitude of harmonic Ritz values
1617 std::vector<MagnitudeType> w(m);
1618
1619 // Sorted order of harmonic Ritz values, also used for DGEEV
1620 std::vector<int> iperm(m);
1621
1622 // Output info
1623 int info = 0;
1624
1625 // Set flag indicating recycle space has been generated this solve
1626 builtRecycleSpace_ = true;
1627
1628 // Solve linear system: H_m^{-H}*E_m where E_m is the last blockSize_ columns of the identity matrix
1629 SDM HHt( HH, Teuchos::TRANS );
1630 Teuchos::RCP<SDM> harmRitzMatrix = Teuchos::rcp( new SDM( m, blockSize_));
1631
1632 //Initialize harmRitzMatrix as E_m
1633 for(int i=0; i<=blockSize_-1; i++) (*harmRitzMatrix)[blockSize_-1-i][harmRitzMatrix->numRows()-1-i] = 1;
1634
1635 //compute harmRitzMatrix <- H_m^{-H}*E_m
1636 lapack.GESV(m, blockSize_, HHt.values(), HHt.stride(), &iperm[0], harmRitzMatrix->values(), harmRitzMatrix->stride(), &info);
1637
1638 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure, "Belos::BlockGCRODRSolMgr::solve(): LAPACK GESV failed to compute a solution.");
1639 // Compute H_m + H_m^{-H}*E_m*H_lbl^{H}*H_lbl
1640 // H_lbl is bottom-right block of H_, which is a blockSize_ x blockSize_ matrix
1641
1642 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl(Teuchos::View, HH, blockSize_, blockSize_, (HH).numRows()-blockSize_, (HH).numCols()-blockSize_ );
1643 Teuchos::SerialDenseMatrix<int, ScalarType> H_lbl_t( H_lbl, Teuchos::TRANS );
1644
1645 { //So that HTemp will fall out of scope
1646
1647 // HH_lbl_t <- H_lbl_t*H_lbl
1648 Teuchos::RCP<SDM> Htemp = Teuchos::null;
1649 Htemp = Teuchos::rcp(new SDM(H_lbl_t.numRows(), H_lbl_t.numCols()));
1650 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, H_lbl_t, H_lbl, zero);
1651 H_lbl_t.assign(*Htemp);
1652 // harmRitzMatrix <- harmRitzMatrix*HH_lbl_t
1653 Htemp = Teuchos::rcp(new SDM(harmRitzMatrix -> numRows(), harmRitzMatrix -> numCols()));
1654 Htemp -> multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, *harmRitzMatrix, H_lbl_t, zero);
1655 harmRitzMatrix -> assign(*Htemp);
1656
1657 // We need to add harmRitzMatrix to the last blockSize_ columns of HH and store the result harmRitzMatrix
1658 int harmColIndex, HHColIndex;
1659 Htemp = Teuchos::rcp(new SDM(Teuchos::Copy,HH,HH.numRows()-blockSize_,HH.numCols()));
1660 for(int i = 0; i<blockSize_; i++){
1661 harmColIndex = harmRitzMatrix -> numCols() - i -1;
1662 HHColIndex = m-i-1;
1663 for(int j=0; j<m; j++) (*Htemp)[HHColIndex][j] += (*harmRitzMatrix)[harmColIndex][j];
1664 }
1665 harmRitzMatrix = Htemp;
1666 }
1667
1668 // Revise to do query for optimal workspace first
1669 // Create simple storage for the left eigenvectors, which we don't care about.
1670
1671 // Size of workspace and workspace for DGEEV
1672 int lwork = -1;
1673 std::vector<ScalarType> work(1);
1674 std::vector<MagnitudeType> rwork(2*m);
1675
1676 const int ldvl = 1;
1677 ScalarType* vl = 0;
1678
1679 // First query GEEV for the optimal workspace size
1680 lapack.GEEV('N', 'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1681 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1682
1683 lwork = std::abs (static_cast<int> (Teuchos::ScalarTraits<ScalarType>::real (work[0])));
1684 work.resize( lwork );
1685
1686 lapack.GEEV('N', 'V', m, harmRitzMatrix->values(), harmRitzMatrix->stride(), &wr[0], &wi[0],
1687 vl, ldvl, vr.values(), vr.stride(), &work[0], lwork, &rwork[0], &info);
1688
1689 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK GEEV failed to compute eigensolutions.");
1690
1691 // Construct magnitude of each harmonic Ritz value
1692 for( int i=0; i<m; ++i ) w[i] = Teuchos::ScalarTraits<MagnitudeType>::squareroot( wr[i]*wr[i] + wi[i]*wi[i] );
1693
1694 this->sort(w, m, iperm);
1695
1696 bool scalarTypeIsComplex = Teuchos::ScalarTraits<ScalarType>::isComplex;
1697
1698 // Select recycledBlocks_ smallest eigenvectors
1699 for( int i=0; i<recycledBlocks_; ++i )
1700 for(int j=0; j<m; j++ )
1701 PP(j,i) = vr(j,iperm[i]);
1702
1703 if(scalarTypeIsComplex==false) {
1704
1705 // Determine exact size for PP (i.e., determine if we need to store an additional vector)
1706 if (wi[iperm[recycledBlocks_-1]] != 0.0) {
1707 int countImag = 0;
1708 for (int i=0; i<recycledBlocks_; ++i )
1709 if (wi[iperm[i]] != 0.0) countImag++;
1710 // Check to see if this count is even or odd:
1711 if (countImag % 2) xtraVec = true;
1712 }
1713
1714 if (xtraVec) { // we need to store one more vector
1715 if (wi[iperm[recycledBlocks_-1]] > 0.0) { // I picked the "real" component
1716 for(int j=0; j<m; ++j ) // so get the "imag" component
1717 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]+1);
1718 }
1719 else{ // I picked the "imag" component
1720 for(int j=0; j<m; ++j ) // so get the "real" component
1721 PP(j,recycledBlocks_) = vr(j,iperm[recycledBlocks_-1]-1);
1722 }
1723 }
1724
1725 }
1726
1727 // Return whether we needed to store an additional vector
1728 if (xtraVec) {
1729 printer_->stream(Debug) << "Recycled " << recycledBlocks_+1 << " vectors" << std::endl;
1730 return recycledBlocks_+1;
1731 }
1732 else {
1733 printer_->stream(Debug) << "Recycled " << recycledBlocks_ << " vectors" << std::endl;
1734 return recycledBlocks_;
1735 }
1736} //end getHarmonicVecsKryl
1737
1738// This method sorts list of n floating-point numbers and return permutation vector
1739template<class ScalarType, class MV, class OP>
1740void BlockGCRODRSolMgr<ScalarType,MV,OP>::sort(std::vector<MagnitudeType>& dlist, int n, std::vector<int>& iperm) {
1741 int l, r, j, i, flag;
1742 int RR2;
1743 MagnitudeType dRR, dK;
1744
1745 // Initialize the permutation vector.
1746 for(j=0;j<n;j++)
1747 iperm[j] = j;
1748
1749 if (n <= 1) return;
1750
1751 l = n / 2 + 1;
1752 r = n - 1;
1753 l = l - 1;
1754 dRR = dlist[l - 1];
1755 dK = dlist[l - 1];
1756
1757 RR2 = iperm[l - 1];
1758 while (r != 0) {
1759 j = l;
1760 flag = 1;
1761 while (flag == 1) {
1762 i = j;
1763 j = j + j;
1764 if (j > r + 1)
1765 flag = 0;
1766 else {
1767 if (j < r + 1)
1768 if (dlist[j] > dlist[j - 1]) j = j + 1;
1769 if (dlist[j - 1] > dK) {
1770 dlist[i - 1] = dlist[j - 1];
1771 iperm[i - 1] = iperm[j - 1];
1772 }
1773 else {
1774 flag = 0;
1775 }
1776 }
1777 }
1778 dlist[i - 1] = dRR;
1779 iperm[i - 1] = RR2;
1780 if (l == 1) {
1781 dRR = dlist [r];
1782 RR2 = iperm[r];
1783 dK = dlist[r];
1784 dlist[r] = dlist[0];
1785 iperm[r] = iperm[0];
1786 r = r - 1;
1787 }
1788 else {
1789 l = l - 1;
1790 dRR = dlist[l - 1];
1791 RR2 = iperm[l - 1];
1792 dK = dlist[l - 1];
1793 }
1794 }
1795 dlist[0] = dRR;
1796 iperm[0] = RR2;
1797} //end sort() definition
1798
1799template<class ScalarType, class MV, class OP>
1801 using Teuchos::RCP;
1802 using Teuchos::rcp;
1803 using Teuchos::rcp_const_cast;
1804
1805 // MLP: NEED TO ADD CHECK IF PARAMETERS ARE SET LATER
1806
1807 ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
1808 ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
1809 std::vector<int> index(numBlocks_+1);
1810
1811 // MLP: EXCEPTION TESTING SHOULD GO HERE
1812
1813 //THE FOLLOWING BLOCK OF CODE IS TO INFORM THE PROBLEM HOW MANY RIGHT HAND
1814 //SIDES WILL BE SOLVED. IF THE CHOSEN BLOCK SIZE IS THE SAME AS THE NUMBER
1815 //OF RIGHT HAND SIDES IN THE PROBLEM THEN WE PASS THE PROBLEM
1816 //INDICES FOR ALL RIGHT HAND SIDES. IF BLOCK SIZE IS GREATER THAN
1817 //THE NUMBER OF RIGHT HAND SIDES AND THE USER HAS ADAPTIVE BLOCK SIZING
1818 //TURNED OFF, THEN THE PROBLEM WILL GENERATE RANDOM RIGHT HAND SIDES SO WE HAVE AS MANY
1819 //RIGHT HAND SIDES AS BLOCK SIZE INDICATES. THIS MAY NEED TO BE CHANGED
1820 //LATER TO ACCOMADATE SMARTER CHOOSING OF FICTITIOUS RIGHT HAND SIDES.
1821 //THIS CODE IS DIRECTLY LIFTED FROM BelosBlockGmresSolMgr.hpp.
1822
1823 // Create indices for the linear systems to be solved.
1824 int startPtr = 0;
1825 int numRHS2Solve = MVT::GetNumberVecs( *(problem_->getRHS()) );
1826 int numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
1827
1828 //currIdx holds indices to all right-hand sides to be solved
1829 //or -1's to indicate that random right-hand sides should be generated
1830 std::vector<int> currIdx;
1831
1832 if ( adaptiveBlockSize_ ) {
1833 blockSize_ = numCurrRHS;
1834 currIdx.resize( numCurrRHS );
1835 for (int i=0; i<numCurrRHS; ++i)
1836 currIdx[i] = startPtr+i;
1837 }
1838 else {
1839 currIdx.resize( blockSize_ );
1840 for (int i=0; i<numCurrRHS; ++i)
1841 currIdx[i] = startPtr+i;
1842 for (int i=numCurrRHS; i<blockSize_; ++i)
1843 currIdx[i] = -1;
1844 }
1845
1846 // Inform the linear problem of the linear systems to solve/generate.
1847 problem_->setLSIndex( currIdx );
1848
1849 //ADD ERROR CHECKING TO MAKE SURE SIZE OF BLOCK KRYLOV SUBSPACE NOT LARGER THAN dim
1850 //ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1851
1852 // reset loss of orthogonality flag
1853 loaDetected_ = false;
1854
1855 // Assume convergence is achieved, then let any failed convergence set this to false.
1856 bool isConverged = true;
1857
1858 // Initialize storage for all state variables
1859 initializeStateStorage();
1860
1861 // Parameter list MLP: WE WILL NEED TO ASSIGN SOME PARAMETER VALUES LATER
1862 Teuchos::ParameterList plist;
1863
1864 while (numRHS2Solve > 0){ //This loops through each block of RHS's being solved
1865 // **************************************************************************************************************************
1866 // Begin initial solver preparations. Either update U,C for new operator or generate an initial space using a cycle of GMRES
1867 // **************************************************************************************************************************
1868 //int prime_iterations; // silence warning about not being used
1869 if(keff > 0){//If there is already a subspace to recycle, then use it
1870 // Update U, C for the new operator
1871
1872// TEUCHOS_TEST_FOR_EXCEPTION(keff < recycledBlocks_,BlockGCRODRSolMgrRecyclingFailure,
1873// "Belos::BlockGCRODRSolMgr::solve(): Requested size of recycled subspace is not consistent with the current recycle subspace.");
1874 printer_->stream(Debug) << " Now solving RHS index " << currIdx[0] << " using recycled subspace of dimension " << keff << std::endl << std::endl;
1875
1876 // Compute image of U_ under the new operator
1877 index.resize(keff);
1878 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1879 RCP<const MV> Utmp = MVT::CloneView( *U_, index );
1880 RCP<MV> Ctmp = MVT::CloneViewNonConst( *C_, index );
1881 problem_->apply( *Utmp, *Ctmp );
1882
1883 RCP<MV> U1tmp = MVT::CloneViewNonConst( *U1_, index );
1884
1885 // Orthogonalize this block
1886 // Get a matrix to hold the orthonormalization coefficients.
1887 SDM Ftmp( Teuchos::View, *F_, keff, keff );
1888 int rank = ortho_->normalize(*Ctmp, rcp(&Ftmp,false));
1889 // Throw an error if we could not orthogonalize this block
1890 TEUCHOS_TEST_FOR_EXCEPTION(rank != keff,BlockGCRODRSolMgrOrthoFailure,"Belos::BlockGCRODRSolMgr::solve(): Failed to compute orthonormal basis for initial recycled subspace.");
1891
1892 // U_ = U_*F^{-1}
1893 // First, compute LU factorization of R
1894 int info = 0;
1895 ipiv_.resize(Ftmp.numRows());
1896 lapack.GETRF(Ftmp.numRows(),Ftmp.numCols(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&info);
1897 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRF failed to compute an LU factorization.");
1898
1899 // Now, form inv(F)
1900 int lwork = Ftmp.numRows();
1901 work_.resize(lwork);
1902 lapack.GETRI(Ftmp.numRows(),Ftmp.values(),Ftmp.stride(),&ipiv_[0],&work_[0],lwork,&info);
1903 TEUCHOS_TEST_FOR_EXCEPTION(info != 0, BlockGCRODRSolMgrLAPACKFailure,"Belos::BlockGCRODRSolMgr::solve(): LAPACK _GETRI failed to invert triangular matrix.");
1904
1905 // U_ = U1_; (via a swap)
1906 MVT::MvTimesMatAddMv( one, *Utmp, Ftmp, zero, *U1tmp );
1907 std::swap(U_, U1_);
1908
1909 // Must reinitialize after swap
1910 index.resize(keff);
1911 for (int ii=0; ii<keff; ++ii) { index[ii] = ii; }
1912 Ctmp = MVT::CloneViewNonConst( *C_, index );
1913 Utmp = MVT::CloneView( *U_, index );
1914
1915 // Compute C_'*R_
1916 SDM Ctr(keff,blockSize_);
1917 problem_->computeCurrPrecResVec( &*R_ );
1918 MVT::MvTransMv( one, *Ctmp, *R_, Ctr );
1919
1920 // Update solution ( x += U_*C_'*R_ )
1921 RCP<MV> update = MVT::Clone( *problem_->getCurrLHSVec(), blockSize_ );
1922 MVT::MvInit( *update, 0.0 );
1923 MVT::MvTimesMatAddMv( one, *Utmp, Ctr, one, *update );
1924 problem_->updateSolution( update, true );
1925
1926 // Update residual norm ( r -= C_*C_'*R_ )
1927 MVT::MvTimesMatAddMv( -one, *Ctmp, Ctr, one, *R_ );
1928
1929 // We recycled space from previous call
1930 //prime_iterations = 0;
1931
1932 // Since we have a previous recycle space, we do not need block_gmres_iter
1933 // and we must allocate V ourselves. See comments in initialize() routine for
1934 // further explanation.
1935 if (V_ == Teuchos::null) {
1936 // Check if there is a multivector to clone from.
1937 Teuchos::RCP<const MV> rhsMV = problem_->getRHS();
1938 V_ = MVT::Clone( *rhsMV, (numBlocks_+1)*blockSize_ );
1939 }
1940 else{
1941 // Generate V_ by cloning itself ONLY if more space is needed.
1942 if (MVT::GetNumberVecs(*V_) < (numBlocks_+1)*blockSize_ ) {
1943 Teuchos::RCP<const MV> tmp = V_;
1944 V_ = MVT::Clone( *tmp, (numBlocks_+1)*blockSize_ );
1945 }
1946 }
1947 } //end if(keff > 0)
1948 else{ //if there was no subspace to recycle, then build one
1949 printer_->stream(Debug) << " No recycled subspace available for RHS index " << std::endl << std::endl;
1950
1951 Teuchos::ParameterList primeList;
1952 //we set this as numBlocks-1 because that is the Krylov dimension we want
1953 //so that the size of the Hessenberg created matches that of what we were expecting.
1954 primeList.set("Num Blocks",numBlocks_-1);
1955 primeList.set("Block Size",blockSize_);
1956 primeList.set("Recycled Blocks",0);
1957 primeList.set("Keep Hessenberg",true);
1958 primeList.set("Initialize Hessenberg",true);
1959
1960 ptrdiff_t dim = MVT::GetGlobalLength( *(problem_->getRHS()) );
1961 if (blockSize_*static_cast<ptrdiff_t>(numBlocks_) > dim) {//if user has selected a total subspace dimension larger than system dimension
1962 ptrdiff_t tmpNumBlocks = 0;
1963 if (blockSize_ == 1)
1964 tmpNumBlocks = dim / blockSize_; // Allow for a good breakdown.
1965 else{
1966 tmpNumBlocks = ( dim - blockSize_) / blockSize_; // Allow for restarting.
1967 printer_->stream(Warnings) << "Belos::BlockGmresSolMgr::solve(): Warning! Requested Krylov subspace dimension is larger than operator dimension!"
1968 << std::endl << "The maximum number of blocks allowed for the Krylov subspace will be adjusted to " << tmpNumBlocks << std::endl;
1969 primeList.set("Num Blocks",Teuchos::as<int>(tmpNumBlocks));
1970 }
1971 }
1972 else{
1973 primeList.set("Num Blocks",numBlocks_-1);
1974 }
1975 //Create Block GMRES iteration object to perform one cycle of GMRES
1976 Teuchos::RCP<BlockGmresIter<ScalarType,MV,OP> > block_gmres_iter;
1977 block_gmres_iter = Teuchos::rcp( new BlockGmresIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,primeList) );
1978
1979 // MLP: ADD LOGIC TO DEAL WITH USER ASKING TO GENERATE A LARGER SPACE THAN dim AS IN HEIDI'S BlockGmresSolMgr CODE (DIDN'T WE ALREADY DO THIS SOMEWHERE?)
1980 block_gmres_iter->setSize( blockSize_, numBlocks_-1 );
1981
1982 //Create the first block in the current BLOCK Krylov basis (using residual)
1983 Teuchos::RCP<MV> V_0;
1984 if (currIdx[blockSize_-1] == -1) {//if augmented with random RHS's
1985 V_0 = MVT::Clone( *(problem_->getInitPrecResVec()), blockSize_ );
1986 problem_->computeCurrPrecResVec( &*V_0 );
1987 }
1988 else {
1989 V_0 = MVT::CloneCopy( *(problem_->getInitPrecResVec()), currIdx );
1990 }
1991
1992 // Get a matrix to hold the orthonormalization coefficients.
1993 Teuchos::RCP<SDM > z_0 = Teuchos::rcp( new SDM( blockSize_, blockSize_ ) );
1994
1995 // Orthonormalize the new V_0
1996 int rank = ortho_->normalize( *V_0, z_0 );
1997 (void) rank; // silence compiler warning for unused variable
1998 // MLP: ADD EXCEPTION IF INITIAL BLOCK IS RANK DEFFICIENT
1999
2000 // Set the new state and initialize the iteration.
2002 newstate.V = V_0;
2003 newstate.z = z_0;
2004 newstate.curDim = 0;
2005 block_gmres_iter->initializeGmres(newstate);
2006
2007 bool primeConverged = false;
2008
2009 try {
2010 printer_->stream(Debug) << " Preparing to Iterate!!!!" << std::endl << std::endl;
2011 block_gmres_iter->iterate();
2012
2013 // *********************
2014 // check for convergence
2015 // *********************
2016 if ( convTest_->getStatus() == Passed ) {
2017 printer_->stream(Debug) << "We converged during the prime the pump stage" << std::endl << std::endl;
2018 primeConverged = !(expConvTest_->getLOADetected());
2019 if ( expConvTest_->getLOADetected() ) {
2020 // we don't have convergence
2021 loaDetected_ = true;
2022 printer_->stream(Warnings) << "Belos::BlockGmresSolMgr::solve(): Warning! Solver has experienced a loss of accuracy!" << std::endl;
2023 }
2024 }
2025 // *******************************************
2026 // check for maximum iterations of block GMRES
2027 // *******************************************
2028 else if( maxIterTest_->getStatus() == Passed ) {
2029 // we don't have convergence
2030 primeConverged = false;
2031 }
2032 // ****************************************************************
2033 // We need to recycle and continue, print a message indicating this
2034 // ****************************************************************
2035 else{ //debug statements
2036 printer_->stream(Debug) << " We did not converge on priming cycle of Block GMRES. Therefore we recycle and restart. " << std::endl << std::endl;
2037 }
2038 } //end try
2039 catch (const GmresIterationOrthoFailure &e) {
2040 // If the block size is not one, it's not considered a lucky breakdown.
2041 if (blockSize_ != 1) {
2042 printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2043 << block_gmres_iter->getNumIters() << std::endl
2044 << e.what() << std::endl;
2045 if (convTest_->getStatus() != Passed)
2046 primeConverged = false;
2047 }
2048 else {
2049 // If the block size is one, try to recover the most recent least-squares solution
2050 block_gmres_iter->updateLSQR( block_gmres_iter->getCurSubspaceDim() );
2051 // Check to see if the most recent least-squares solution yielded convergence.
2052 sTest_->checkStatus( &*block_gmres_iter );
2053 if (convTest_->getStatus() != Passed)
2054 isConverged = false;
2055 }
2056 } // end catch (const GmresIterationOrthoFailure &e)
2057 catch (const std::exception &e) {
2058 printer_->stream(Errors) << "Error! Caught std::exception in BlockGmresIter::iterate() at iteration "
2059 << block_gmres_iter->getNumIters() << std::endl
2060 << e.what() << std::endl;
2061 throw;
2062 }
2063
2064 // Record number of iterations in generating initial recycle spacec
2065 //prime_iterations = block_gmres_iter->getNumIters();//instantiated here because it is not needed outside of else{} scope; we'll see if this is true or not
2066
2067 // Update the linear problem.
2068 RCP<MV> update = block_gmres_iter->getCurrentUpdate();
2069 problem_->updateSolution( update, true );
2070
2071 // Update the block residual
2072
2073 problem_->computeCurrPrecResVec( &*R_ );
2074
2075 // Get the state.
2076 newstate = block_gmres_iter->getState();
2077 int p = newstate.curDim;
2078 (void) p; // silence compiler warning for unused variable
2079 H_->assign(*(newstate.H));//IS THIS A GOOD IDEA? I DID IT SO MEMBERS WOULD HAVE ACCESS TO H.
2080 // Get new Krylov vectors, store in V_
2081
2082 // Since the block_gmres_iter returns the state
2083 // with const RCP's we need to cast newstate.V as
2084 // a non const RCP. This is okay because block_gmres_iter
2085 // is about to fall out of scope, so we are simply
2086 // rescuing *newstate.V from being destroyed so we can
2087 // use it for future block recycled GMRES cycles
2088 V_ = rcp_const_cast<MV>(newstate.V);
2089 newstate.V.release();
2090 //COMPUTE NEW RECYCLE SPACE SOMEHOW
2091 buildRecycleSpaceKryl(keff, block_gmres_iter);
2092 printer_->stream(Debug) << "Generated recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2093
2094 // Return to outer loop if the priming solve
2095 // converged, set the next linear system.
2096 if (primeConverged) {
2097 /* MLP: POSSIBLY INCORRECT CODE WE ARE REPLACING *********************************
2098 // Inform the linear problem that we are
2099 // finished with this block linear system.
2100 problem_->setCurrLS();
2101
2102 // Update indices for the linear systems to be solved.
2103 // KMS: Fix the numRHS2Solve; Copy from Heidi's BlockGmres code
2104 numRHS2Solve -= 1;
2105 printer_->stream(Debug) << numRHS2Solve << " Right hand sides left to solve" << std::endl;
2106 if ( numRHS2Solve > 0 ) {
2107 currIdx[0]++;
2108 // Set the next indices
2109 problem_->setLSIndex( currIdx );
2110 }
2111 else {
2112 currIdx.resize( numRHS2Solve );
2113 }
2114 *****************************************************************************/
2115
2116 // Inform the linear problem that we are finished with this block linear system.
2117 problem_->setCurrLS();
2118
2119 // Update indices for the linear systems to be solved.
2120 startPtr += numCurrRHS;
2121 numRHS2Solve -= numCurrRHS;
2122 if ( numRHS2Solve > 0 ) {
2123 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2124 if ( adaptiveBlockSize_ ) {
2125 blockSize_ = numCurrRHS;
2126 currIdx.resize( numCurrRHS );
2127 for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2128 }
2129 else {
2130 currIdx.resize( blockSize_ );
2131 for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2132 for (int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2133 }
2134 // Set the next indices.
2135 problem_->setLSIndex( currIdx );
2136 }
2137 else {
2138 currIdx.resize( numRHS2Solve );
2139 }
2140 continue;//Begin solving the next block of RHS's
2141 } //end if (primeConverged)
2142 } //end else [if(keff > 0)]
2143
2144 // *****************************************
2145 // Initial subspace update/construction done
2146 // Start cycles of recycled GMRES
2147 // *****************************************
2148 Teuchos::ParameterList blockgcrodrList;
2149 blockgcrodrList.set("Num Blocks",numBlocks_);
2150 blockgcrodrList.set("Block Size",blockSize_);
2151 blockgcrodrList.set("Recycled Blocks",keff);
2152
2153 Teuchos::RCP<BlockGCRODRIter<ScalarType,MV,OP> > block_gcrodr_iter;
2154 block_gcrodr_iter = Teuchos::rcp( new BlockGCRODRIter<ScalarType,MV,OP>(problem_,printer_,outputTest_,ortho_,blockgcrodrList) );
2156
2157 index.resize( blockSize_ );
2158 for(int ii = 0; ii < blockSize_; ii++) index[ii] = ii;
2159 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2160
2161 // MLP: MVT::SetBlock(*R_,index,*V0);
2162 MVT::Assign(*R_,*V0);
2163
2164 index.resize(keff);//resize to get appropriate recycle space vectors
2165 for(int i=0; i < keff; i++){ index[i] = i;};
2166 B_ = rcp(new SDM(Teuchos::View, *G_, keff, numBlocks_*blockSize_, 0, keff));
2167 H_ = rcp(new SDM(Teuchos::View, *G_, (numBlocks_-1)*blockSize_ + blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2168
2169 newstate.V = V_;
2170 newstate.B= B_;
2171 newstate.U = MVT::CloneViewNonConst(*U_, index);
2172 newstate.C = MVT::CloneViewNonConst(*C_, index);
2173 newstate.H = H_;
2174 newstate.curDim = blockSize_;
2175 block_gcrodr_iter -> initialize(newstate);
2176
2177 int numRestarts = 0;
2178
2179 while(1){//Each execution of this loop is a cycle of block GCRODR
2180 try{
2181 block_gcrodr_iter -> iterate();
2182
2183 // ***********************
2184 // Check Convergence First
2185 // ***********************
2186 if( convTest_->getStatus() == Passed ) {
2187 // we have convergence
2188 break;//from while(1)
2189 } //end if converged
2190
2191 // ***********************************
2192 // Check if maximum iterations reached
2193 // ***********************************
2194 else if(maxIterTest_->getStatus() == Passed ){
2195 // no convergence; hit maxit
2196 isConverged = false;
2197 break; // from while(1)
2198 } //end elseif reached maxit
2199
2200 // **********************************************
2201 // Check if subspace full; do we need to restart?
2202 // **********************************************
2203 else if (block_gcrodr_iter->getCurSubspaceDim() == block_gcrodr_iter->getMaxSubspaceDim()){
2204
2205 // Update recycled space even if we have reached max number of restarts
2206
2207 // Update linear problem
2208 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2209 problem_->updateSolution(update, true);
2210 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2211
2212 printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2213 // NOTE: If we have hit the maximum number of restarts, then we will quit
2214 if(numRestarts >= maxRestarts_) {
2215 isConverged = false;
2216 break; //from while(1)
2217 } //end if max restarts
2218
2219 numRestarts++;
2220
2221 printer_ -> stream(Debug) << " Performing restart number " << numRestarts << " of " << maxRestarts_ << std::endl << std::endl;
2222
2223 // Create the restart vector (first block in the current Krylov basis)
2224 problem_->computeCurrPrecResVec( &*R_ );
2225 index.resize( blockSize_ );
2226 for (int ii=0; ii<blockSize_; ++ii) index[ii] = ii;
2227 Teuchos::RCP<MV> V0 = MVT::CloneViewNonConst( *V_, index );
2228 MVT::SetBlock(*R_,index,*V0);
2229
2230 // Set the new state and initialize the solver.
2232 index.resize( numBlocks_*blockSize_ );
2233 for (int ii=0; ii<(numBlocks_*blockSize_); ++ii) index[ii] = ii;
2234 restartState.V = MVT::CloneViewNonConst( *V_, index );
2235 index.resize( keff );
2236 for (int ii=0; ii<keff; ++ii) index[ii] = ii;
2237 restartState.U = MVT::CloneViewNonConst( *U_, index );
2238 restartState.C = MVT::CloneViewNonConst( *C_, index );
2239 B_ = rcp(new SDM(Teuchos::View, *G_, keff, (numBlocks_-1)*blockSize_, 0, keff));
2240 H_ = rcp(new SDM(Teuchos::View, *G_, numBlocks_*blockSize_, (numBlocks_-1)*blockSize_, keff ,keff ));
2241 restartState.B = B_;
2242 restartState.H = H_;
2243 restartState.curDim = blockSize_;
2244 block_gcrodr_iter->initialize(restartState);
2245
2246 } //end else if need to restart
2247
2248 // ****************************************************************
2249 // We returned from iterate(), but none of our status tests passed.
2250 // Something is wrong, and it is probably our fault.
2251 // ****************************************************************
2252 else {
2253 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"Belos::BlockGCRODRSolMgr::solve(): Invalid return from BlockGCRODRIter::iterate().");
2254 } //end else (no status test passed)
2255
2256 }//end try
2257 catch(const BlockGCRODRIterOrthoFailure &e){ //there was an exception
2258 // Try to recover the most recent least-squares solution
2259 block_gcrodr_iter->updateLSQR( block_gcrodr_iter->getCurSubspaceDim() );
2260 // Check to see if the most recent least-squares solution yielded convergence.
2261 sTest_->checkStatus( &*block_gcrodr_iter );
2262 if (convTest_->getStatus() != Passed) isConverged = false;
2263 break;
2264 } // end catch orthogonalization failure
2265 catch(const std::exception &e){
2266 printer_->stream(Errors) << "Error! Caught exception in BlockGCRODRIter::iterate() at iteration "
2267 << block_gcrodr_iter->getNumIters() << std::endl
2268 << e.what() << std::endl;
2269 throw;
2270 } //end catch standard exception
2271 } //end while(1)
2272
2273 // Compute the current solution.
2274 // Update the linear problem.
2275 Teuchos::RCP<MV> update = block_gcrodr_iter->getCurrentUpdate();
2276 problem_->updateSolution( update, true );
2277
2278 /* MLP: POSSIBLY INCORRECT CODE WE ARE REPLACING *********************************
2279 // Inform the linear problem that we are finished with this block linear system.
2280 problem_->setCurrLS();
2281
2282 // Update indices for the linear systems to be solved.
2283 numRHS2Solve -= 1;
2284 if ( numRHS2Solve > 0 ) {
2285 currIdx[0]++;
2286 // Set the next indices.
2287 problem_->setLSIndex( currIdx );
2288 }
2289 else {
2290 currIdx.resize( numRHS2Solve );
2291 }
2292 ******************************************************************************/
2293
2294 // Inform the linear problem that we are finished with this block linear system.
2295 problem_->setCurrLS();
2296
2297 // Update indices for the linear systems to be solved.
2298 startPtr += numCurrRHS;
2299 numRHS2Solve -= numCurrRHS;
2300 if ( numRHS2Solve > 0 ) {
2301 numCurrRHS = ( numRHS2Solve < blockSize_) ? numRHS2Solve : blockSize_;
2302 if ( adaptiveBlockSize_ ) {
2303 blockSize_ = numCurrRHS;
2304 currIdx.resize( numCurrRHS );
2305 for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2306 }
2307 else {
2308 currIdx.resize( blockSize_ );
2309 for (int i=0; i<numCurrRHS; ++i) currIdx[i] = startPtr+i;
2310 for (int i=numCurrRHS; i<blockSize_; ++i) currIdx[i] = -1;
2311 }
2312 // Set the next indices.
2313 problem_->setLSIndex( currIdx );
2314 }
2315 else {
2316 currIdx.resize( numRHS2Solve );
2317 }
2318
2319 // If we didn't build a recycle space this solve but ran at least k iterations, force build of new recycle space
2320 if (!builtRecycleSpace_) {
2321 buildRecycleSpaceAugKryl(block_gcrodr_iter);
2322 printer_->stream(Debug) << " Generated new recycled subspace using RHS index " << currIdx[0] << " of dimension " << keff << std::endl << std::endl;
2323 }
2324 }//end while (numRHS2Solve > 0)
2325
2326 // print final summary
2327 sTest_->print( printer_->stream(FinalSummary) );
2328
2329 // print timing information
2330 #ifdef BELOS_TEUCHOS_TIME_MONITOR
2331 // Calling summarize() can be expensive, so don't call unless the user wants to print out timing details.
2332 // summarize() will do all the work even if it's passed a "black hole" output stream.
2333 if (verbosity_ & TimingDetails) Teuchos::TimeMonitor::summarize( printer_->stream(TimingDetails) );
2334 #endif
2335 // get iteration information for this solve
2336 numIters_ = maxIterTest_->getNumIters();
2337
2338 // get residual information for this solve
2339 const std::vector<MagnitudeType>* pTestValues = NULL;
2340 pTestValues = impConvTest_->getTestValue();
2341 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues == NULL, std::logic_error,
2342 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2343 "getTestValue() method returned NULL. Please report this bug to the "
2344 "Belos developers.");
2345 TEUCHOS_TEST_FOR_EXCEPTION(pTestValues->size() < 1, std::logic_error,
2346 "Belos::BlockGCRODRSolMgr::solve(): The implicit convergence test's "
2347 "getTestValue() method returned a vector of length zero. Please report "
2348 "this bug to the Belos developers.");
2349 // FIXME (mfh 12 Dec 2011) Does pTestValues really contain the
2350 // achieved tolerances for all vectors in the current solve(), or
2351 // just for the vectors from the last deflation?
2352 achievedTol_ = *std::max_element (pTestValues->begin(), pTestValues->end());
2353
2354 if (!isConverged) return Unconverged; // return from BlockGCRODRSolMgr::solve()
2355 return Converged; // return from BlockGCRODRSolMgr::solve()
2356} //end solve()
2357
2358} //End Belos Namespace
2359
2360#endif /* BELOS_BLOCK_GCRODR_SOLMGR_HPP */
Belos concrete class for performing the block, flexible GMRES iteration.
Belos concrete class for performing the block GCRO-DR (block GMRES with recycling) iteration.
Belos concrete class for performing the block GMRES iteration.
Belos header file which uses auto-configuration information to include necessary C++ headers.
Pure virtual base class which augments the basic interface for a Gmres linear solver 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.
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
Thrown when an LAPACK call fails.
BlockGCRODRSolMgrLAPACKFailure(const std::string &what_arg)
Thrown when the linear problem was not set up correctly.
BlockGCRODRSolMgrLinearProblemFailure(const std::string &what_arg)
Thrown when the solution manager was unable to orthogonalize the basis vectors.
BlockGCRODRSolMgrOrthoFailure(const std::string &what_arg)
Thrown if any problem occurs in using or creating the recycle subspace.
BlockGCRODRSolMgrRecyclingFailure(const std::string &what_arg)
A solver manager for the Block GCRO-DR (Block Recycling GMRES) linear solver.
Teuchos::LAPACK< int, ScalarType > lapack
void setParameters(const Teuchos::RCP< Teuchos::ParameterList > &params)
Set the parameters the solver should use to solve the linear problem.
BlockGCRODRSolMgr()
Default constructor.
Teuchos::RCP< const Teuchos::ParameterList > defaultParams_
Default parameter list.
Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > problem_
The current linear problem to solve.
Teuchos::RCP< const Teuchos::ParameterList > getCurrentParameters() const
Get a parameter list containing the current parameters for this object.
void buildRecycleSpaceKryl(int &keff, Teuchos::RCP< BlockGmresIter< ScalarType, MV, OP > > block_gmres_iter)
Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > ortho_
Orthogonalization manager.
std::string description() const
A description of the Block GCRODR solver manager.
Teuchos::RCP< StatusTestMaxIters< ScalarType, MV, OP > > maxIterTest_
MagnitudeType achievedTol() const
Get the residual for the most recent call to solve().
OrthoManagerFactory< ScalarType, MV, OP > ortho_factory_type
bool isSet_
Whether setParameters() successfully finished setting parameters.
Teuchos::ScalarTraits< MagnitudeType > SMT
OperatorTraits< ScalarType, MV, OP > OPT
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > expConvTest_
Teuchos::ScalarTraits< MagnitudeType > MT
Teuchos::SerialDenseVector< int, ScalarType > SDV
Teuchos::RCP< StatusTestOutput< ScalarType, MV, OP > > outputTest_
Teuchos::RCP< StatusTestGenResNorm< ScalarType, MV, OP > > impConvTest_
bool loaDetected_
Whether a loss of accuracy was detected during the solve.
Teuchos::RCP< const Teuchos::ParameterList > getValidParameters() const
Get a parameter list containing the valid parameters for this object.
static const bool adaptiveBlockSize_default_
MultiVecTraits< ScalarType, MV > MVT
void reset(const ResetType type)
Performs a reset of the solver manager specified by the ResetType.
static const std::string recycleMethod_default_
Teuchos::RCP< Teuchos::ParameterList > params_
This solver's current parameter list.
Teuchos::SerialDenseMatrix< int, ScalarType > SDM
virtual ~BlockGCRODRSolMgr()
Destructor.
Teuchos::RCP< std::ostream > outputStream_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > sTest_
std::vector< ScalarType > tau_
void buildRecycleSpaceAugKryl(Teuchos::RCP< BlockGCRODRIter< ScalarType, MV, OP > > gcrodr_iter)
int getHarmonicVecsAugKryl(int keff, int m, const SDM &HH, const Teuchos::RCP< const MV > &VV, SDM &PP)
bool isLOADetected() const
Whether a loss of accuracy was detected during the most recent solve.
ReturnType solve()
Solve the current linear problem.
Teuchos::ScalarTraits< ScalarType >::magnitudeType MagnitudeType
Teuchos::RCP< OutputManager< ScalarType > > printer_
void setProblem(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem)
Set the linear problem to solve on the next call to solve().
Teuchos::ScalarTraits< ScalarType > SCT
int getNumIters() const
Get the iteration count for the most recent call to solve().
std::vector< ScalarType > work_
Teuchos::RCP< StatusTest< ScalarType, MV, OP > > convTest_
Teuchos::RCP< Teuchos::Time > timerSolve_
Timer for solve().
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get current linear problem being solved for in this object.
bool builtRecycleSpace_
Whether we have generated or regenerated a recycle space yet this solve.
int getHarmonicVecsKryl(int m, const SDM &HH, SDM &PP)
void sort(std::vector< MagnitudeType > &dlist, int n, std::vector< int > &iperm)
ortho_factory_type orthoFactory_
Factory for creating MatOrthoManager subclass instances.
This class implements the block GMRES iteration, where a block Krylov subspace is constructed....
An implementation of the Belos::MatOrthoManager that performs orthogonalization using (potentially) m...
void setDepTol(const MagnitudeType dep_tol)
Set parameter for re-orthogonalization threshhold.
GmresIterationOrthoFailure is thrown when the GmresIteration object is unable to compute independent ...
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.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
The Belos::SolverManager is a templated virtual base class that defines the basic interface that any ...
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.
@ StatusTestDetails
@ FinalSummary
@ TimingDetails
ReturnType
Whether the Belos solve converged for all linear systems.
@ Unconverged
OutputType
Style of output used to display status test information.
ScaleType
The type of scaling to use on the residual norm value.
ResetType
How to reset the solver.
@ RecycleSubspace
Structure to contain pointers to BlockGCRODRIter state variables.
Teuchos::RCP< MV > V
The current Krylov basis.
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.
int curDim
The current dimension of the reduction.
Structure to contain pointers to GmresIteration state variables.
Teuchos::RCP< const MV > V
The current Krylov basis.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > z
The current right-hand side of the least squares system RY = Z.
int curDim
The current dimension of the reduction.
Teuchos::RCP< const Teuchos::SerialDenseMatrix< int, ScalarType > > H
The current Hessenberg matrix.