Stokhos Package Browser (Single Doxygen Collection) Version of the Day
Loading...
Searching...
No Matches
Belos_PseudoBlockGmresIter_MP_Vector.hpp
Go to the documentation of this file.
1//@HEADER
2// ************************************************************************
3//
4// Belos: Block Linear Solvers Package
5// Copyright 2004 Sandia Corporation
6//
7// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8// the U.S. Government retains certain rights in this software.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38//
39// ************************************************************************
40//@HEADER
41
42#ifndef BELOS_PSEUDO_BLOCK_GMRES_ITER_MP_VECTOR_HPP
43#define BELOS_PSEUDO_BLOCK_GMRES_ITER_MP_VECTOR_HPP
44
54#include "BelosPseudoBlockGmresIter.hpp"
55
69namespace Belos {
70
71 template<class Storage, class MV, class OP>
72 class PseudoBlockGmresIter<Sacado::MP::Vector<Storage>, MV, OP> :
73 virtual public Iteration<Sacado::MP::Vector<Storage>, MV, OP> {
74
75 public:
76
77 //
78 // Convenience typedefs
79 //
81 typedef MultiVecTraits<ScalarType,MV> MVT;
82 typedef OperatorTraits<ScalarType,MV,OP> OPT;
83 typedef Teuchos::ScalarTraits<ScalarType> SCT;
84 typedef typename SCT::magnitudeType MagnitudeType;
85 typedef Teuchos::ScalarTraits<typename Storage::value_type> SVT;
86
88
89
98 PseudoBlockGmresIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
99 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
100 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
101 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
102 Teuchos::ParameterList &params );
103
105 virtual ~PseudoBlockGmresIter() = default;
107
108
110
111
133 void iterate();
134
156 void initialize(const PseudoBlockGmresIterState<ScalarType,MV> & newstate);
157
162 {
163 PseudoBlockGmresIterState<ScalarType,MV> empty;
164 initialize(empty);
165 }
166
174 PseudoBlockGmresIterState<ScalarType,MV> getState() const {
175 PseudoBlockGmresIterState<ScalarType,MV> state;
176 state.curDim = curDim_;
177 state.V.resize(numRHS_);
178 state.H.resize(numRHS_);
179 state.Z.resize(numRHS_);
180 state.sn.resize(numRHS_);
181 state.cs.resize(numRHS_);
182 for (int i=0; i<numRHS_; ++i) {
183 state.V[i] = V_[i];
184 state.H[i] = H_[i];
185 state.Z[i] = Z_[i];
186 state.sn[i] = sn_[i];
187 state.cs[i] = cs_[i];
188 }
189 return state;
190 }
191
193
194
196
197
199 int getNumIters() const { return iter_; }
200
202 void resetNumIters( int iter = 0 ) { iter_ = iter; }
203
221 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
222
224
229 Teuchos::RCP<MV> getCurrentUpdate() const;
230
232
235 void updateLSQR( int dim = -1 );
236
238 int getCurSubspaceDim() const {
239 if (!initialized_) return 0;
240 return curDim_;
241 }
242
244 int getMaxSubspaceDim() const { return numBlocks_; }
245
247
248
250
251
253 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
254
256 int getBlockSize() const { return 1; }
257
259 void setBlockSize(int blockSize) {
260 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
261 "Belos::PseudoBlockGmresIter::setBlockSize(): Cannot use a block size that is not one.");
262 }
263
265 int getNumBlocks() const { return numBlocks_; }
266
268 void setNumBlocks(int numBlocks);
269
271 bool isInitialized() { return initialized_; }
272
274
275 private:
276
277 //
278 // Classes inputed through constructor that define the linear problem to be solved.
279 //
280 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
281 const Teuchos::RCP<OutputManager<ScalarType> > om_;
282 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
283 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
284
285 //
286 // Algorithmic parameters
287 //
288 // numRHS_ is the current number of linear systems being solved.
290 // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
292
293 // Mask used to store whether the ensemble GMRES has faced lucky breakdown or not.
295
296 // Storage for QR factorization of the least squares system.
297 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > sn_;
298 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,MagnitudeType> > > cs_;
299
300 // Pointers to a work vector used to improve aggregate performance.
301 Teuchos::RCP<MV> U_vec_, AU_vec_;
302
303 // Pointers to the current right-hand side and solution multivecs being solved for.
304 Teuchos::RCP<MV> cur_block_rhs_, cur_block_sol_;
305
306 //
307 // Current solver state
308 //
309 // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
310 // is capable of running; _initialize is controlled by the initialize() member method
311 // For the implications of the state of initialized_, please see documentation for initialize()
313
314 // Current subspace dimension, and number of iterations performed.
315 int curDim_, iter_;
316
317 //
318 // State Storage
319 //
320 std::vector<Teuchos::RCP<MV> > V_;
321 //
322 // Projected matrices
323 // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
324 //
325 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > H_;
326 //
327 // QR decomposition of Projected matrices for solving the least squares system HY = B.
328 // R_: Upper triangular reduction of H
329 // Z_: Q applied to right-hand side of the least squares system
330 std::vector<Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > R_;
331 std::vector<Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > > Z_;
332
333 // Tolerance for ensemble breakdown
334 typename SVT::magnitudeType breakDownTol_;
335
336#ifdef BELOS_TEUCHOS_TIME_MONITOR
337 Teuchos::RCP<Teuchos::Time> timerUpdateLSQR_, timerSolveLSQR_;
338#endif
339 };
340
342 // Constructor.
343 template<class StorageType, class MV, class OP>
345 PseudoBlockGmresIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
346 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
347 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
348 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
349 Teuchos::ParameterList &params ):
350 lp_(problem),
351 om_(printer),
352 stest_(tester),
353 ortho_(ortho),
354 numRHS_(0),
355 numBlocks_(0),
356 initialized_(false),
357 curDim_(0),
358 iter_(0),
359 breakDownTol_(params.get("Ensemble Breakdown Tolerance", 0.0))
360 {
361 // Get the maximum number of blocks allowed for each Krylov subspace
362 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument,
363 "Belos::PseudoBlockGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
364 int nb = Teuchos::getParameter<int>(params, "Num Blocks");
365
366 setNumBlocks( nb );
367 }
368
370 // Set the block size and make necessary adjustments.
371 template <class StorageType, class MV, class OP>
372 void PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::setNumBlocks (int numBlocks)
373 {
374 // This routine only allocates space; it doesn't not perform any computation
375 // any change in size will invalidate the state of the solver.
376
377 TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0, std::invalid_argument, "Belos::PseudoBlockGmresIter::setNumBlocks was passed a non-positive argument.");
378
379 numBlocks_ = numBlocks;
380 curDim_ = 0;
381
382 initialized_ = false;
383 }
384
386 // Get the current update from this subspace.
387 template <class StorageType, class MV, class OP>
388 Teuchos::RCP<MV> PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::getCurrentUpdate() const
389 {
390#ifdef BELOS_TEUCHOS_TIME_MONITOR
391 Teuchos::TimeMonitor updateTimer( *timerSolveLSQR_ );
392#endif
393 //
394 // If this is the first iteration of the Arnoldi factorization,
395 // there is no update, so return Teuchos::null.
396 //
397 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
398 if (curDim_==0) {
399 return currentUpdate;
400 } else {
401 currentUpdate = MVT::Clone(*(V_[0]), numRHS_);
402 std::vector<int> index(1), index2(curDim_);
403 for (int i=0; i<curDim_; ++i) {
404 index2[i] = i;
405 }
406 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
407 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
408 Teuchos::BLAS<int,ScalarType> blas;
409
410 for (int i=0; i<numRHS_; ++i) {
411 index[0] = i;
412 Teuchos::RCP<MV> cur_block_copy_vec = MVT::CloneViewNonConst( *currentUpdate, index );
413 //
414 // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
415 //
416 Teuchos::SerialDenseVector<int,ScalarType> y( Teuchos::Copy, Z_[i]->values(), curDim_ );
417 //
418 // Solve the least squares problem and compute current solutions.
419 //
420 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
421 Teuchos::NON_UNIT_DIAG, curDim_, 1, one,
422 H_[i]->values(), H_[i]->stride(), y.values(), y.stride() );
423
424
425 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_[i], index2 );
426 MVT::MvTimesMatAddMv( one, *Vjp1, y, zero, *cur_block_copy_vec );
427 }
428 }
429 return currentUpdate;
430 }
431
432
434 // Get the native residuals stored in this iteration.
435 // Note: No residual vector will be returned by Gmres.
436 template <class StorageType, class MV, class OP>
437 Teuchos::RCP<const MV>
439 getNativeResiduals (std::vector<MagnitudeType> *norms) const
440 {
441 typedef typename Teuchos::ScalarTraits<ScalarType> STS;
442
443 if (norms)
444 { // Resize the incoming std::vector if necessary. The type
445 // cast avoids the compiler warning resulting from a signed /
446 // unsigned integer comparison.
447 if (static_cast<int> (norms->size()) < numRHS_)
448 norms->resize (numRHS_);
449
450 Teuchos::BLAS<int, ScalarType> blas;
451 for (int j = 0; j < numRHS_; ++j)
452 {
453 const ScalarType curNativeResid = (*Z_[j])(curDim_);
454 (*norms)[j] = STS::magnitude (curNativeResid);
455 }
456 }
457 return Teuchos::null;
458 }
459
460
461 template <class StorageType, class MV, class OP>
462 void
464 initialize (const PseudoBlockGmresIterState<ScalarType,MV> & newstate)
465 {
466 using Teuchos::RCP;
467
468 // (Re)set the number of right-hand sides, by interrogating the
469 // current LinearProblem to solve.
470 this->numRHS_ = MVT::GetNumberVecs (*(lp_->getCurrLHSVec()));
471
472#ifdef BELOS_TEUCHOS_TIME_MONITOR
473 std::stringstream ss;
474 ss << "Belos";
475
476 std::string updateLabel = ss.str() + ": Update LSQR";
477 timerUpdateLSQR_ = Teuchos::TimeMonitor::getNewCounter(updateLabel);
478
479 std::string solveLabel = ss.str() + ": Solve LSQR";
480 timerSolveLSQR_ = Teuchos::TimeMonitor::getNewCounter(solveLabel);
481#endif
482
483 // NOTE: In PseudoBlockGmresIter, V and Z are required!!!
484 // Inconsistent multivectors widths and lengths will not be tolerated, and
485 // will be treated with exceptions.
486 //
487 std::string errstr ("Belos::PseudoBlockGmresIter::initialize(): "
488 "Specified multivectors must have a consistent "
489 "length and width.");
490
491 // Check that newstate has V and Z arrays with nonzero length.
492 TEUCHOS_TEST_FOR_EXCEPTION((int)newstate.V.size()==0 || (int)newstate.Z.size()==0,
493 std::invalid_argument,
494 "Belos::PseudoBlockGmresIter::initialize(): "
495 "V and/or Z was not specified in the input state; "
496 "the V and/or Z arrays have length zero.");
497
498 // In order to create basis multivectors, we have to clone them
499 // from some already existing multivector. We require that at
500 // least one of the right-hand side B and left-hand side X in the
501 // LinearProblem be non-null. Thus, we can clone from either B or
502 // X. We prefer to close from B, since B is in the range of the
503 // operator A and the basis vectors should also be in the range of
504 // A (the first basis vector is a scaled residual vector).
505 // However, if B is not given, we will try our best by cloning
506 // from X.
507 RCP<const MV> lhsMV = lp_->getLHS();
508 RCP<const MV> rhsMV = lp_->getRHS();
509
510 // If the right-hand side is null, we make do with the left-hand
511 // side, otherwise we use the right-hand side.
512 RCP<const MV> vectorInBasisSpace = rhsMV.is_null() ? lhsMV : rhsMV;
513 //RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
514
515 TEUCHOS_TEST_FOR_EXCEPTION(vectorInBasisSpace.is_null(),
516 std::invalid_argument,
517 "Belos::PseudoBlockGmresIter::initialize(): "
518 "The linear problem to solve does not specify multi"
519 "vectors from which we can clone basis vectors. The "
520 "right-hand side(s), left-hand side(s), or both should "
521 "be nonnull.");
522
523 // Check the new dimension is not more that the maximum number of
524 // allowable blocks.
525 TEUCHOS_TEST_FOR_EXCEPTION(newstate.curDim > numBlocks_+1,
526 std::invalid_argument,
527 errstr);
528 curDim_ = newstate.curDim;
529
530 // Initialize the state storage. If the subspace has not be
531 // initialized before, generate it using the right-hand side or
532 // left-hand side from the LinearProblem lp_ to solve.
533 V_.resize(numRHS_);
534 for (int i=0; i<numRHS_; ++i) {
535 // Create a new vector if we need to. We "need to" if the
536 // current vector V_[i] is null, or if it doesn't have enough
537 // columns.
538 if (V_[i].is_null() || MVT::GetNumberVecs(*V_[i]) < numBlocks_ + 1) {
539 V_[i] = MVT::Clone (*vectorInBasisSpace, numBlocks_ + 1);
540 }
541 // Check that the newstate vector newstate.V[i] has dimensions
542 // consistent with those of V_[i].
543 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetGlobalLength(*newstate.V[i]) != MVT::GetGlobalLength(*V_[i]),
544 std::invalid_argument, errstr );
545 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V[i]) < newstate.curDim,
546 std::invalid_argument, errstr );
547 //
548 // If newstate.V[i] and V_[i] are not identically the same
549 // vector, then copy newstate.V[i] into V_[i].
550 //
551 int lclDim = MVT::GetNumberVecs(*newstate.V[i]);
552 if (newstate.V[i] != V_[i]) {
553 // Only copy over the first block and print a warning.
554 if (curDim_ == 0 && lclDim > 1) {
555 om_->stream(Warnings)
556 << "Belos::PseudoBlockGmresIter::initialize(): the solver was "
557 << "initialized with a kernel of " << lclDim
558 << std::endl
559 << "The block size however is only " << 1
560 << std::endl
561 << "The last " << lclDim - 1 << " vectors will be discarded."
562 << std::endl;
563 }
564 std::vector<int> nevind (curDim_ + 1);
565 for (int j = 0; j < curDim_ + 1; ++j)
566 nevind[j] = j;
567
568 RCP<const MV> newV = MVT::CloneView (*newstate.V[i], nevind);
569 RCP<MV> lclV = MVT::CloneViewNonConst( *V_[i], nevind );
570 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
571 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
572 MVT::MvAddMv (one, *newV, zero, *newV, *lclV);
573
574 // Done with local pointers
575 lclV = Teuchos::null;
576 }
577 }
578
579
580 // Check size of Z
581 Z_.resize(numRHS_);
582 for (int i=0; i<numRHS_; ++i) {
583 // Create a vector if we need to.
584 if (Z_[i] == Teuchos::null) {
585 Z_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>() );
586 }
587 if (Z_[i]->length() < numBlocks_+1) {
588 Z_[i]->shapeUninitialized(numBlocks_+1, 1);
589 }
590
591 // Check that the newstate vector is consistent.
592 TEUCHOS_TEST_FOR_EXCEPTION(newstate.Z[i]->numRows() < curDim_, std::invalid_argument, errstr);
593
594 // Put data into Z_, make sure old information is not still hanging around.
595 if (newstate.Z[i] != Z_[i]) {
596 if (curDim_==0)
597 Z_[i]->putScalar();
598
599 Teuchos::SerialDenseVector<int,ScalarType> newZ(Teuchos::View,newstate.Z[i]->values(),curDim_+1);
600 Teuchos::RCP<Teuchos::SerialDenseVector<int,ScalarType> > lclZ;
601 lclZ = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(Teuchos::View,Z_[i]->values(),curDim_+1) );
602 lclZ->assign(newZ);
603
604 // Done with local pointers
605 lclZ = Teuchos::null;
606 }
607 }
608
609
610 // Check size of H
611 H_.resize(numRHS_);
612 for (int i=0; i<numRHS_; ++i) {
613 // Create a matrix if we need to.
614 if (H_[i] == Teuchos::null) {
615 H_[i] = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
616 }
617 if (H_[i]->numRows() < numBlocks_+1 || H_[i]->numCols() < numBlocks_) {
618 H_[i]->shapeUninitialized(numBlocks_+1, numBlocks_);
619 }
620
621 // Put data into H_ if it exists, make sure old information is not still hanging around.
622 if ((int)newstate.H.size() == numRHS_) {
623
624 // Check that the newstate matrix is consistent.
625 TEUCHOS_TEST_FOR_EXCEPTION((newstate.H[i]->numRows() < curDim_ || newstate.H[i]->numCols() < curDim_), std::invalid_argument,
626 "Belos::PseudoBlockGmresIter::initialize(): Specified Hessenberg matrices must have a consistent size to the current subspace dimension");
627
628 if (newstate.H[i] != H_[i]) {
629 // H_[i]->putScalar();
630
631 Teuchos::SerialDenseMatrix<int,ScalarType> newH(Teuchos::View,*newstate.H[i],curDim_+1, curDim_);
632 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > lclH;
633 lclH = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>(Teuchos::View,*H_[i],curDim_+1, curDim_) );
634 lclH->assign(newH);
635
636 // Done with local pointers
637 lclH = Teuchos::null;
638 }
639 }
640 }
641
643 // Reinitialize storage for least squares solve
644 //
645 cs_.resize(numRHS_);
646 sn_.resize(numRHS_);
647
648 // Copy over rotation angles if they exist
649 if ((int)newstate.cs.size() == numRHS_ && (int)newstate.sn.size() == numRHS_) {
650 for (int i=0; i<numRHS_; ++i) {
651 if (cs_[i] != newstate.cs[i])
652 cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(*newstate.cs[i]) );
653 if (sn_[i] != newstate.sn[i])
654 sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(*newstate.sn[i]) );
655 }
656 }
657
658 // Resize or create the vectors as necessary
659 for (int i=0; i<numRHS_; ++i) {
660 if (cs_[i] == Teuchos::null)
661 cs_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,MagnitudeType>(numBlocks_+1) );
662 else
663 cs_[i]->resize(numBlocks_+1);
664 if (sn_[i] == Teuchos::null)
665 sn_[i] = Teuchos::rcp( new Teuchos::SerialDenseVector<int,ScalarType>(numBlocks_+1) );
666 else
667 sn_[i]->resize(numBlocks_+1);
668 }
669
670 // the solver is initialized
671 initialized_ = true;
672
673 /*
674 if (om_->isVerbosity( Debug ) ) {
675 // Check almost everything here
676 CheckList chk;
677 chk.checkV = true;
678 chk.checkArn = true;
679 chk.checkAux = true;
680 om_->print( Debug, accuracyCheck(chk, ": after initialize()") );
681 }
682 */
683
684 }
685
686
688 // Iterate until the status test informs us we should stop.
689 template <class StorageType, class MV, class OP>
691 {
692 //
693 // Allocate/initialize data structures
694 //
695 if (initialized_ == false) {
696 initialize();
697 }
698
699 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
700 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
701
702 // Compute the current search dimension.
703 int searchDim = numBlocks_;
704 //
705 // Associate each initial block of V_[i] with U_vec[i]
706 // Reset the index vector (this might have been changed if there was a restart)
707 //
708 std::vector<int> index(1);
709 std::vector<int> index2(1);
710 index[0] = curDim_;
711 Teuchos::RCP<MV> U_vec = MVT::Clone( *V_[0], numRHS_ );
712
713 // Create AU_vec to hold A*U_vec.
714 Teuchos::RCP<MV> AU_vec = MVT::Clone( *V_[0], numRHS_ );
715
716 for (int i=0; i<numRHS_; ++i) {
717 index2[0] = i;
718 Teuchos::RCP<const MV> tmp_vec = MVT::CloneView( *V_[i], index );
719 Teuchos::RCP<MV> U_vec_view = MVT::CloneViewNonConst( *U_vec, index2 );
720 MVT::MvAddMv( one, *tmp_vec, zero, *tmp_vec, *U_vec_view );
721 }
722
724 // iterate until the status test tells us to stop.
725 //
726 // also break if our basis is full
727 //
728 while (stest_->checkStatus(this) != Passed && curDim_ < searchDim) {
729 iter_++;
730 //
731 // Apply the operator to _work_vector
732 //
733
734 lp_->apply( *U_vec, *AU_vec );
735 //
736 //
737 // Resize index.
738 //
739 int num_prev = curDim_+1;
740 index.resize( num_prev );
741 for (int i=0; i<num_prev; ++i) {
742 index[i] = i;
743 }
744 //
745 // Orthogonalize next Krylov vector for each right-hand side.
746 //
747 for (int i=0; i<numRHS_; ++i) {
748 //
749 // Get previous Krylov vectors.
750 //
751 Teuchos::RCP<const MV> V_prev = MVT::CloneView( *V_[i], index );
752 Teuchos::Array< Teuchos::RCP<const MV> > V_array( 1, V_prev );
753
754 //
755 // Get a view of the new candidate std::vector.
756 //
757 index2[0] = i;
758 Teuchos::RCP<MV> V_new = MVT::CloneViewNonConst( *AU_vec, index2 );
759
760 //
761 // Get a view of the current part of the upper-hessenberg matrix.
762 //
763 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > h_new
764 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
765 ( Teuchos::View, *H_[i], num_prev, 1, 0, curDim_ ) );
766 Teuchos::Array< Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > > h_array( 1, h_new );
767
768 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > r_new
769 = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>
770 ( Teuchos::View, *H_[i], 1, 1, num_prev, curDim_ ) );
771
772 //
773 // Orthonormalize the new block of the Krylov expansion
774 // NOTE: Rank deficiencies are not checked because this is a single-std::vector Krylov method.
775 //
776 ortho_->projectAndNormalize( *V_new, h_array, r_new, V_array );
777
778 //
779 // NOTE: V_new is a copy of the iter+1 vector in V_[i], so the normalized vector has to be
780 // be copied back in when V_new is changed.
781 //
782 index2[0] = curDim_+1;
783 Teuchos::RCP<MV> tmp_vec = MVT::CloneViewNonConst( *V_[i], index2 );
784 MVT::MvAddMv( one, *V_new, zero, *V_new, *tmp_vec );
785 }
786 //
787 // Now _AU_vec is the new _U_vec, so swap these two vectors.
788 // NOTE: This alleviates the need for allocating a vector for AU_vec each iteration.
789 //
790 Teuchos::RCP<MV> tmp_AU_vec = U_vec;
791 U_vec = AU_vec;
792 AU_vec = tmp_AU_vec;
793 //
794 // V has been extended, and H has been extended.
795 //
796 // Update the QR factorization of the upper Hessenberg matrix
797 //
798 {
799#ifdef BELOS_TEUCHOS_TIME_MONITOR
800 Teuchos::TimeMonitor updateTimer( *timerUpdateLSQR_ );
801#endif
802 updateLSQR();
803 }
804 //
805 // Update basis dim and release all pointers.
806 //
807 curDim_ += 1;
808 //
809 /*
810 // When required, monitor some orthogonalities
811 if (om_->isVerbosity( Debug ) ) {
812 // Check almost everything here
813 CheckList chk;
814 chk.checkV = true;
815 chk.checkArn = true;
816 om_->print( Debug, accuracyCheck(chk, ": after local update") );
817 }
818 else if (om_->isVerbosity( OrthoDetails ) ) {
819 CheckList chk;
820 chk.checkV = true;
821 om_->print( OrthoDetails, accuracyCheck(chk, ": after local update") );
822 }
823 */
824
825 } // end while (statusTest == false)
826
827 }
828
830 // Update the least squares solution for each right-hand side.
831 template<class StorageType, class MV, class OP>
832 void PseudoBlockGmresIter<Sacado::MP::Vector<StorageType>,MV,OP>::updateLSQR( int dim )
833 {
834 // Get correct dimension based on input "dim"
835 // Remember that ortho failures result in an exit before updateLSQR() is called.
836 // Therefore, it is possible that dim == curDim_.
837 int curDim = curDim_;
838 if (dim >= curDim_ && dim < getMaxSubspaceDim()) {
839 curDim = dim;
840 }
841
842 int i, j;
843 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
844 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
845
846 Teuchos::BLAS<int, ScalarType> blas;
847
848 for (i=0; i<numRHS_; ++i) {
849 //
850 // Update the least-squares QR for each linear system.
851 //
852 // QR factorization of Least-Squares system with Givens rotations
853 //
854 for (j=0; j<curDim; j++) {
855 //
856 // Apply previous Givens rotations to new column of Hessenberg matrix
857 //
858 blas.ROT( 1, &(*H_[i])(j,curDim), 1, &(*H_[i])(j+1, curDim), 1, &(*cs_[i])[j], &(*sn_[i])[j] );
859 }
860 //
861 // Calculate new Givens rotation
862 //
863 if ( curDim == 0)
864 {
865 // Initialize the lucky_breakdown_ mask to take account of perfect initial guess
866 lucky_breakdown_ = ((*Z_[i])(curDim) == zero);
867 }
868 auto lucky_breakdown_tmp = ((*H_[i])(curDim+1,curDim) == zero);
869 mask_assign(lucky_breakdown_,(*H_[i])(curDim,curDim)) = one;
870 lucky_breakdown_ = lucky_breakdown_ || lucky_breakdown_tmp;
871 //Teuchos::my_GivensRotator<ScalarType> GR;
872 blas.ROTG( &(*H_[i])(curDim,curDim), &(*H_[i])(curDim+1,curDim), &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
873 (*H_[i])(curDim+1,curDim) = zero;
874 //
875 // Update RHS w/ new transformation
876 //
877 blas.ROT( 1, &(*Z_[i])(curDim), 1, &(*Z_[i])(curDim+1), 1, &(*cs_[i])[curDim], &(*sn_[i])[curDim] );
878 }
879
880 } // end updateLSQR()
881
882} // end Belos namespace
883
884#endif /* BELOS_PSEUDO_BLOCK_GMRES_ITER_MP_VECTOR_HPP */
KOKKOS_INLINE_FUNCTION MaskedAssign< scalar > mask_assign(bool b, scalar *s)
const LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
int getCurSubspaceDim() const
Get the dimension of the search subspace used to generate the current solution to the linear problem.
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
PseudoBlockGmresIter(const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< OutputManager< ScalarType > > &printer, const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > &tester, const Teuchos::RCP< MatOrthoManager< ScalarType, MV, OP > > &ortho, Teuchos::ParameterList &params)
PseudoBlockGmresIter constructor with linear problem, solver utilities, and parameter list of solver ...
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the "native" residual vectors.
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > > Z_
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, ScalarType > > > sn_
void iterate()
This method performs block Gmres iterations until the status test indicates the need to stop or an er...
std::vector< Teuchos::RCP< Teuchos::SerialDenseVector< int, MagnitudeType > > > cs_
std::vector< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > H_
std::vector< Teuchos::RCP< Teuchos::SerialDenseMatrix< int, ScalarType > > > R_
PseudoBlockGmresIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
void initialize(const PseudoBlockGmresIterState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, providing a complete state.
void updateLSQR(int dim=-1)
Method for updating QR factorization of upper Hessenberg matrix.
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
bool isInitialized()
States whether the solver has been initialized or not.
void initialize()
Initialize the solver with the initial vectors from the linear problem or random data.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
This class implements the pseudo-block GMRES iteration, where a block Krylov subspace is constructed ...
const IndexType const IndexType const IndexType const IndexType const ValueType const ValueType ValueType * y
Definition: csr_vector.h:267