Belos Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
BelosBlockGCRODRIter.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_BLOCK_GCRODR_ITER_HPP
43#define BELOS_BLOCK_GCRODR_ITER_HPP
44
45
50#include "BelosConfigDefs.hpp"
51#include "BelosTypes.hpp"
52
56#include "BelosStatusTest.hpp"
59
60#include "Teuchos_BLAS.hpp"
61#include "Teuchos_SerialDenseMatrix.hpp"
62#include "Teuchos_SerialDenseVector.hpp"
63#include "Teuchos_ScalarTraits.hpp"
64#include "Teuchos_ParameterList.hpp"
65#include "Teuchos_TimeMonitor.hpp"
66
67// MLP
68#include <unistd.h>
69
84namespace Belos{
85
87
88
94 template <class ScalarType, class MV>
102
104 Teuchos::RCP<MV> V;
105
107 Teuchos::RCP<MV> U, C;
108
114 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > H;
115
118 Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > B;
119
120 BlockGCRODRIterState() : curDim(0), V(Teuchos::null),
121 U(Teuchos::null), C(Teuchos::null),
122 H(Teuchos::null), B(Teuchos::null)
123 {}
124
125 };
126
128
129
130
132
133
147 public:
148 BlockGCRODRIterInitFailure(const std::string& what_arg) : BelosError(what_arg) {}
149 };
150
159 public:
160 BlockGCRODRIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
161 };
162
164
165
166 template<class ScalarType, class MV, class OP>
167 class BlockGCRODRIter : virtual public Iteration<ScalarType,MV,OP> {
168 public:
169
170 //
171 //Convenience typedefs
172 //
175 typedef Teuchos::ScalarTraits<ScalarType> SCT;
176 typedef typename SCT::magnitudeType MagnitudeType;
177 typedef Teuchos::SerialDenseMatrix<int,ScalarType> SDM;
178 typedef Teuchos::SerialDenseVector<int,ScalarType> SDV;
179
181
182
191 BlockGCRODRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem,
192 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
193 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
194 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
195 Teuchos::ParameterList &params );
196
198 virtual ~BlockGCRODRIter() {};
200
202
203
225 void iterate();
226
249 void initialize() {
251 initialize(empty);
252 }
253
258
268 state.curDim = curDim_;
269 state.V = V_;
270 state.U = U_;
271 state.C = C_;
272 state.H = H_;
273 state.B = B_;
274 return state;
275 }
277
279
280
282 bool isInitialized(){ return initialized_;};
283
285 int getNumIters() const { return iter_; };
286
288 void resetNumIters( int iter = 0 ) { iter_ = iter; };
289
292 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
293
295
300 Teuchos::RCP<MV> getCurrentUpdate() const;
301
302
304
305
307
308
309
311 const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; };
312
314 int getNumBlocks() const { return numBlocks_; }
315
317 int getBlockSize() const { return blockSize_; };
318
320 int getCurSubspaceDim() const {
321 if (!initialized_) return 0;
322 return curDim_;
323 };
324
326 int getMaxSubspaceDim() const { return numBlocks_*blockSize_; };
327
329 int getRecycledBlocks() const { return recycledBlocks_; };
330
332
333
335
336
337 void updateLSQR( int dim = -1);
338
340 void setBlockSize(int blockSize){ blockSize_ = blockSize; }
341
343 void setRecycledBlocks(int recycledBlocks) { setSize( recycledBlocks, numBlocks_ ); };
344
346 void setNumBlocks(int numBlocks) { setSize( recycledBlocks_, numBlocks ); };
347
349 void setSize( int recycledBlocks, int numBlocks ) {
350 // only call resize if size changed
351 if ( (recycledBlocks_ != recycledBlocks) || (numBlocks_ != numBlocks) ) {
352 recycledBlocks_ = recycledBlocks;
353 numBlocks_ = numBlocks;
354 cs_.sizeUninitialized( numBlocks_ );
355 sn_.sizeUninitialized( numBlocks_ );
356 Z_.shapeUninitialized( numBlocks_*blockSize_, blockSize_ );
357 }
358 }
359
361
362 private:
363
364 //
365 // Internal methods
366 //
367
368
369
370 //Classes inputed through constructor that define the linear problem to be solved
371 //
372 const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > lp_;
373 const Teuchos::RCP<OutputManager<ScalarType> > om_;
374 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > stest_;
375 const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
376
377 //
378 //Algorithmic Parameters
379 //
380
381 //numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
382 //blockSize_ is the number of columns in each block Krylov vector.
384
385 //boolean vector indicating which right-hand sides we care about
386 //when we are testing residual norms. THIS IS NOT IMPLEMENTED. RIGHT NOW JUST
387 //SELECTS ALL RIGHT HANDS SIDES FOR NORM COMPUTATION.
388 std::vector<bool> trueRHSIndices_;
389
390 // recycledBlocks_ is the size of the allocated space for the recycled subspace, in vectors.
392
393 //Storage for QR factorization of the least squares system if using plane rotations.
395 Teuchos::SerialDenseVector<int,MagnitudeType> cs_;
396
397 //Storage for QR factorization of the least squares system if using Householder reflections
398 //Per block Krylov vector, we actually construct a 2*blockSize_ by 2*blockSize_ matrix which
399 //is the product of all Householder transformations for that block. This has been shown to yield
400 //speed ups without losing accuracy because we can apply previous Householder transformations
401 //with BLAS3 operations.
402 std::vector< SDM >House_;
404
405 //
406 //Current Solver State
407 //
408 //initialized_ specifies that the basis vectors have been initialized and the iterate() routine
409 //is capable of running; _initialize is controlled by the initialize() member method
410 //For the implications of the state of initialized_, please see documentation for initialize()
412
413 // Current subspace dimension, number of iterations performed, and number of iterations performed in this cycle.
415
416 //
417 // Recycled Krylov Method Storage
418 //
419
420
422 Teuchos::RCP<MV> V_;
423
425 Teuchos::RCP<MV> U_, C_;
426
427
429
430
434 Teuchos::RCP<SDM > H_;
435
439 Teuchos::RCP<SDM > B_;
440
447 Teuchos::RCP<SDM> R_;
448
451
453
454 // File stream variables to use Mike Parks' Matlab output codes.
455 std::ofstream ofs;
456 char filename[30];
457
458 };//End BlockGCRODRIter Class Definition
459
461 //Constructor.
462 template<class ScalarType, class MV, class OP>
464 const Teuchos::RCP<OutputManager<ScalarType> > &printer,
465 const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
466 const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
467 Teuchos::ParameterList &params ):lp_(problem),
468 om_(printer), stest_(tester), ortho_(ortho) {
469 numBlocks_ = 0;
470 blockSize_ = 0;
471 recycledBlocks_ = 0;
472 initialized_ = false;
473 curDim_ = 0;
474 iter_ = 0;
475 lclIter_ = 0;
476 V_ = Teuchos::null;
477 U_ = Teuchos::null;
478 C_ = Teuchos::null;
479 H_ = Teuchos::null;
480 B_ = Teuchos::null;
481 R_ = Teuchos::null;
482 // Get the maximum number of blocks allowed for this Krylov subspace
483 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument, "Belos::BlockGCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
484 int nb = Teuchos::getParameter<int>(params, "Num Blocks");
485
486 TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,"Belos::BlockGCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
487 int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
488
489 TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument, "Belos::BlockGCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
490 TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument, "Belos::BlockGCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
491
492
493 int bs = Teuchos::getParameter<int>(params, "Block Size");
494
495 TEUCHOS_TEST_FOR_EXCEPTION(bs <= 0, std::invalid_argument, "Belos::BlockGCRODRIter() the block size was passed a non-postitive argument.");
496
497
498 numBlocks_ = nb;
499 recycledBlocks_ = rb;
500 blockSize_ = bs;
501
502 //INITIALIZE ENTRIES OF trueRHSIndices_ TO CHECK EVERY NORM FOR NOW. LATER, THE USER
503 //SHOULD SPECIFY WHICH RIGHT HAND SIDES ARE IMPORTANT FOR CONVERGENCE TESTING
505 int i;
506 for(i=0; i<blockSize_; i++){
507 trueRHSIndices_[i] = true;
508 }
509
510 //THIS MAKES SPACE FOR GIVENS ROTATIONS BUT IN REALITY WE NEED TO DO TESTING ON BLOCK SIZE
511 //AND CHOOSE BETWEEN GIVENS ROTATIONS AND HOUSEHOLDER TRANSFORMATIONS.
512 cs_.sizeUninitialized( numBlocks_+1 );
513 sn_.sizeUninitialized( numBlocks_+1 );
514 Z_.shapeUninitialized( (numBlocks_+1)*blockSize_,blockSize_ );
515
516 House_.resize(numBlocks_);
517
518 for(i=0; i<numBlocks_;i++){
519 House_[i].shapeUninitialized(2*blockSize_, 2*blockSize_);
520 }
521 }//End Constructor Definition
522
524 // Iterate until the status test informs us we should stop.
525 template <class ScalarType, class MV, class OP>
527 TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, BlockGCRODRIterInitFailure,"Belos::BlockGCRODRIter::iterate(): GCRODRIter class not initialized." );
528
529// MLP
530//sleep(1);
531//std::cout << "Calling setSize" << std::endl;
532 // Force call to setsize to ensure internal storage is correct dimension
533 setSize( recycledBlocks_, numBlocks_ );
534
535 Teuchos::RCP<MV> Vnext;
536 Teuchos::RCP<const MV> Vprev;
537 std::vector<int> curind(blockSize_);
538
539 // z_ must be zeroed out in order to compute Givens rotations correctly
540 Z_.putScalar(0.0);
541
542 // Orthonormalize the new V_0
543 for(int i = 0; i<blockSize_; i++){curind[i] = i;};
544
545// MLP
546//sleep(1);
547//std::cout << "Calling normalize" << std::endl;
548 Vnext = MVT::CloneViewNonConst(*V_,curind);
549 //Orthonormalize Initial Columns
550 //Store orthogonalization coefficients in Z0
551 Teuchos::RCP<SDM > Z0 =
552 Teuchos::rcp( new SDM(blockSize_,blockSize_) );
553 int rank = ortho_->normalize(*Vnext,Z0);
554
555// MLP
556//sleep(1);
557//std::cout << "Assigning Z" << std::endl;
558 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockGCRODRIterOrthoFailure, "Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank at the initial step.");
559 // Copy Z0 into the leading blockSize_ by blockSize_ block of Z_
560 Teuchos::RCP<SDM > Z_block = Teuchos::rcp( new SDM(Teuchos::View, Z_, blockSize_,blockSize_) );
561 Z_block->assign(*Z0);
562
563 std::vector<int> prevind(blockSize_*(numBlocks_ + 1));
564
566 // iterate until the status test tells us to stop.
567 //
568 // also break if the basis is full
569 //
570 while( (stest_->checkStatus(this) != Passed) && (curDim_+blockSize_-1) < (numBlocks_*blockSize_)) {
571 lclIter_++;
572 iter_++;
573//KMS
574//std::cout << "Iter=" << iter_ << std::endl << "lclIter=" << lclIter_ << std::endl;
575
576 int HFirstCol = curDim_-blockSize_;//First column of H we need view of
577 int HLastCol = HFirstCol + blockSize_-1 ;//last column of H we need a view of
578 int HLastOrthRow = HLastCol;//Last row of H we will put orthog coefficients in
579 int HFirstNormRow = HLastOrthRow + 1;//First row of H where normalization matrix goes
580//KMS
581//std::cout << "curDim_ = " << curDim_ << ", HFirstCol = " << HFirstCol << ", HLastCol = " << HLastCol <<", HLastOrthRow = " << HLastOrthRow << ", HFirstNormRow = " << HFirstNormRow << std::endl;
582 // Get next basis indices
583 for(int i = 0; i< blockSize_; i++){
584 curind[i] = curDim_ + i;
585 }
586 Vnext = MVT::CloneViewNonConst(*V_,curind);
587
588 //Get a view of the previous block Krylov vector.
589 //This is used for orthogonalization and for computing V^H K H
590 // Get next basis indices
591 for(int i = 0; i< blockSize_; i++){
592 curind[blockSize_ - 1 - i] = curDim_ - i - 1;
593 }
594 Vprev = MVT::CloneView(*V_,curind);
595 // Compute the next vector in the Krylov basis: Vnext = Op*Vprev
596 lp_->apply(*Vprev,*Vnext);
597 Vprev = Teuchos::null;
598
599 //First, remove the recycled subspace (C) from Vnext and put coefficients in B.
600
601 //Get a view of the matrix B and put the pointer into an array
602 //Put a pointer to C in another array
603 Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
604
605 Teuchos::RCP<SDM >
606 subB = Teuchos::rcp( new SDM ( Teuchos::View,*B_,recycledBlocks_,blockSize_,0, HFirstCol ) );
607
608 Teuchos::Array<Teuchos::RCP<SDM > > AsubB;
609 AsubB.append( subB );
610 // Project out the recycled subspace.
611 ortho_->project( *Vnext, AsubB, C );
612 //Now, remove block Krylov Subspace from Vnext and store coefficients in H_ and R_
613
614 // Get a view of all the previous vectors
615 prevind.resize(curDim_);
616 for (int i=0; i<curDim_; i++) { prevind[i] = i; }
617 Vprev = MVT::CloneView(*V_,prevind);
618 Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
619
620 // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
621 Teuchos::RCP<SDM> subH = Teuchos::rcp( new SDM ( Teuchos::View,*H_,curDim_,blockSize_,0,HFirstCol ) );
622 Teuchos::Array<Teuchos::RCP<SDM > > AsubH;
623 AsubH.append( subH );
624 // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
625 Teuchos::RCP<SDM > subR = Teuchos::rcp( new SDM ( Teuchos::View,*H_,blockSize_,blockSize_,HFirstNormRow,HFirstCol ) );
626 // Project out the previous Krylov vectors and normalize the next vector.
627 int rank = ortho_->projectAndNormalize(*Vnext,AsubH,subR,AVprev);
628
629 // Copy over the coefficients to R just in case we run into an error.
630 SDM subR2( Teuchos::View,*R_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
631 SDM subH2( Teuchos::View,*H_,(lclIter_+1)*blockSize_,blockSize_,0,HFirstCol);
632 subR2.assign(subH2);
633
634 TEUCHOS_TEST_FOR_EXCEPTION(rank != blockSize_,BlockGCRODRIterOrthoFailure, "Belos::BlockGCRODRIter::iterate(): couldn't generate basis of full rank.");
635
636 // Update the QR factorization of the upper Hessenberg matrix
637 updateLSQR();
638
639 // Update basis dim
640 curDim_ = curDim_ + blockSize_;
641
642
643
644 }//end while(stest_->checkStatus(this) ~= Passed && curDim_+1 <= numBlocks_*blockSize_)
645
646 }//end iterate() defintition
647
649 //Initialize this iteration object.
650 template <class ScalarType, class MV, class OP>
652 if (newstate.V != Teuchos::null && newstate.H != Teuchos::null) {
653 curDim_ = newstate.curDim;
654 V_ = newstate.V;
655 U_ = newstate.U;
656 C_ = newstate.C;
657 H_ = newstate.H;
658 B_ = newstate.B;
659 lclIter_ = 0;//resets the count of local iterations for the new cycle
660 R_ = Teuchos::rcp(new SDM(H_->numRows(), H_->numCols() )); //R_ should look like H but point to separate memory
661
662 //All Householder product matrices start out as identity matrices.
663 //We construct an identity from which to copy.
664 SDM Identity(2*blockSize_, 2*blockSize_);
665 for(int i=0;i<2*blockSize_; i++){
666 Identity[i][i] = 1;
667 }
668 for(int i=0; i<numBlocks_;i++){
669 House_[i].assign(Identity);
670 }
671 }
672 else {
673 TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have V initialized.");
674 TEUCHOS_TEST_FOR_EXCEPTION(newstate.H == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): BlockGCRODRIterState does not have H initialized.");
675 }
676 // the solver is initialized
677 initialized_ = true;
678 }//end initialize() defintition
679
681 //Get the native residuals stored in this iteration.
682 //This function will only compute the native residuals for
683 //right-hand sides we are interested in, as dictated by
684 //std::vector<int> trueRHSIndices_ (THIS IS NOT YET IMPLEMENTED. JUST GETS ALL RESIDUALS)
685 //A norm of -1 is entered for all residuals about which we do not care.
686 template <class ScalarType, class MV, class OP>
687 Teuchos::RCP<const MV>
688 BlockGCRODRIter<ScalarType,MV,OP>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const
689 {
690 //
691 // NOTE: Make sure the incoming std::vector is the correct size!
692 //
693 if (norms != NULL) {
694 if (static_cast<int> (norms->size()) < blockSize_) {
695 norms->resize( blockSize_ );
696 }
697 Teuchos::BLAS<int,ScalarType> blas;
698 for (int j=0; j<blockSize_; j++) {
699 if(trueRHSIndices_[j]){
700 (*norms)[j] = blas.NRM2( blockSize_, &Z_(curDim_-blockSize_+j, j), 1);
701 }
702 else{
703 (*norms)[j] = -1;
704 }
705 }
706 return Teuchos::null;
707 } else { // norms is NULL
708 // FIXME If norms is NULL, return residual vectors.
709 return Teuchos::null;
710 }
711 }//end getNativeResiduals() definition
712
714 //Get the current update from this subspace.
715 template <class ScalarType, class MV, class OP>
717 //
718 // If this is the first iteration of the Arnoldi factorization,
719 // there is no update, so return Teuchos::null.
720 //
721 Teuchos::RCP<MV> currentUpdate = Teuchos::null;
722//KMS if(curDim_==0) {
723 if(curDim_<=blockSize_) {
724 return currentUpdate;
725 }
726 else{
727 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
728 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
729 Teuchos::BLAS<int,ScalarType> blas;
730 currentUpdate = MVT::Clone( *V_, blockSize_ );
731 //
732 // Make a view and then copy the RHS of the least squares problem. DON'T OVERWRITE IT!
733 //
734 SDM Y( Teuchos::Copy, Z_, curDim_-blockSize_, blockSize_ );
735 Teuchos::RCP<SDM> Rtmp = Teuchos::rcp(new SDM(Teuchos::View, *R_, curDim_, curDim_-blockSize_));
736//KMS
737//sleep(1);
738//std::cout << "Before TRSM" << std::endl;
739//sleep(1);
740//std::cout << "The size of Rtmp is " << Rtmp -> numRows() << " by " << Rtmp -> numCols() << std::endl;
741//std::cout << "The size of Y is " << Y.numRows() << " by " << Y.numCols() << std::endl;
742//std::cout << "blockSize_ = " << blockSize_ << std::endl;
743//std::cout << "curDim_ = " << curDim_ << std::endl;
744//std::cout << "curDim_ - blockSize_ = " << curDim_ - blockSize_ << std::endl;
745 //
746 // Solve the least squares problem.
747 // Observe that in calling TRSM, we use the value
748 // curDim_ -blockSize_. This is because curDim_ has
749 // already been incremented but the proper size of R is still
750 // based on the previous value.
751 //
752 blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
753 Teuchos::NON_UNIT_DIAG, curDim_-blockSize_, blockSize_, one,
754 Rtmp->values(), Rtmp->stride(), Y.values(), Y.stride() );
755//KMS
756//sleep(1);
757//std::cout << "After TRSM" << std::endl;
758 //
759 // Compute the current update from the Krylov basis; V(:,1:curDim_)*y.
760 //
761 std::vector<int> index(curDim_-blockSize_);
762 for ( int i=0; i<curDim_-blockSize_; i++ ) index[i] = i;
763 Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
764 MVT::MvTimesMatAddMv( one, *Vjp1, Y, zero, *currentUpdate );
765
766
767
768 //
769 // Add in portion of update from recycled subspace U; U(:,1:recycledBlocks_)*B*y.
770 //
771 if (U_ != Teuchos::null) {
772 SDM z(recycledBlocks_,blockSize_);
773 SDM subB( Teuchos::View, *B_, recycledBlocks_, curDim_-blockSize_ );
774 z.multiply( Teuchos::NO_TRANS, Teuchos::NO_TRANS, one, subB, Y, zero );
775
776 //std::cout << (*U_).MyLength() << " " << (*U_).NumVectors() << " " << subB.numRows() << " " << subB.numCols() << " " << Y.numRows() << " " << Y.numCols()<< " " << curDim_ << " " << recycledBlocks_;
777 MVT::MvTimesMatAddMv( -one, *U_, z, one, *currentUpdate );
778 }
779 }
780
781
782
783 return currentUpdate;
784 }//end getCurrentUpdate() definition
785
786 template<class ScalarType, class MV, class OP>
788
789 int i;
790 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
791 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
792
793 using Teuchos::rcp;
794
795 // Get correct dimension based on input "dim"
796 // Remember that ortho failures result in an exit before updateLSQR() is called.
797 // Therefore, it is possible that dim == curDim_.
798 int curDim = curDim_;
799 if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) ){
800 curDim = dim;
801 }
802
803 Teuchos::BLAS<int, ScalarType> blas;
804
805 if(blockSize_ == 1){//if only one right-hand side then use Givens rotations
806 //
807 // Apply previous rotations and compute new rotations to reduce upper-Hessenberg
808 // system to upper-triangular form.
809 //
810 // QR factorization of Least-Squares system with Givens rotations
811 //
812 for (i=0; i<curDim-1; i++) {
813 //
814 // Apply previous Givens rotations to new column of Hessenberg matrix
815 //
816 blas.ROT( 1, &(*R_)(i,curDim-1), 1, &(*R_)(i+1, curDim-1), 1, &cs_[i], &sn_[i] );
817 }
818 //
819 // Calculate new Givens rotation
820 //
821 blas.ROTG( &(*R_)(curDim-1,curDim-1), &(*R_)(curDim,curDim-1), &cs_[curDim-1], &sn_[curDim-1] );
822 (*R_)(curDim,curDim-1) = zero;
823 //
824 // Update RHS w/ new transformation
825 //
826 blas.ROT( 1, &Z_(curDim-1,0), 1, &Z_(curDim,0), 1, &cs_[curDim-1], &sn_[curDim-1] );
827 }
828 else{//if multiple right-hand sides then use Householder transormations
829 //
830 //apply previous reflections and compute new reflections to reduce upper-Hessenberg
831 //system to upper-triagular form.
832
833 //In Matlab, applying the reflection to a matrix M would look like
834 // M_refl = M - 2*v_refl*(v_refl'*M)/(norm(v_refl)^2)
835
836 //In order to take advantage of BLAS while applying reflections to a matrix, we
837 //perform it in a two step process
838 //1. workvec = M'*v_refl {using BLAS.GEMV()}
839 //2. M_refl = M_refl - 2*v_refl*workvec'/(norm(v_refl)^2) {using BLAS.GER()}
840
841 Teuchos::RCP< SDM > workmatrix = Teuchos::null;//matrix of column vectors onto which we apply the reflection
842 Teuchos::RCP< SDV > workvec = Teuchos::null;//where we store the result of the first step of the 2-step reflection process
843 Teuchos::RCP<SDV> v_refl = Teuchos::null;//the reflection vector
844 int R_colStart = curDim_-blockSize_;
845 Teuchos::RCP< SDM >Rblock = Teuchos::null;
846
847 //
848 //Apply previous reflections
849 //
850 for(i=0; i<lclIter_-1; i++){
851 int R_rowStart = i*blockSize_;
852 //get a view of the part of R_ effected by these reflections.
853 Teuchos::RCP< SDM > RblockCopy = rcp(new SDM (Teuchos::Copy, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
854 Teuchos::RCP< SDM > RblockView = rcp(new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, R_rowStart, R_colStart));
855 blas.GEMM(Teuchos::NO_TRANS,Teuchos::NO_TRANS, 2*blockSize_,blockSize_,2*blockSize_,one,House_[i].values(),House_[i].stride(), RblockCopy->values(),RblockCopy -> stride(), zero, RblockView->values(),RblockView -> stride());
856
857 }
858
859
860 //Get a view of last 2*blockSize entries of entire block to
861 //generate new reflections.
862 Rblock = rcp(new SDM (Teuchos::View, *R_, 2*blockSize_,blockSize_, curDim_-blockSize_, curDim_-blockSize_));
863
864 //Calculate and apply the new reflections
865 for(i=0; i<blockSize_; i++){
866 //
867 //Calculating Reflection
868 //
869 int curcol = (lclIter_ - 1)*blockSize_ + i;//current column of R_
870 int lclCurcol = i;//current column of Rblock
871 ScalarType signDiag = (*R_)(curcol,curcol) / Teuchos::ScalarTraits<ScalarType>::magnitude((*R_)(curcol,curcol));
872
873 // Norm of the vector to be reflected.
874 // BLAS returns a ScalarType, but it really should be a magnitude.
875 ScalarType nvs = blas.NRM2(blockSize_+1,&((*R_)[curcol][curcol]),1);
876 ScalarType alpha = -signDiag*nvs;
877
878 //norm of reflection vector which is just the vector being reflected
879 //i.e. v = R_(curcol:curcol+blockSize_,curcol))
880 //v_refl = v - alpha*e1
881 //norm(v_refl) = norm(v) + alpha^2 - 2*v*alpha
882 //store in v_refl
883
884 // Beware, nvs should have a zero imaginary part (since
885 // it is a norm of a vector), but it may not due to rounding
886 // error.
887 //nvs = nvs + alpha*alpha - 2*(*R_)(curcol,curcol)*alpha;
888 //(*R_)(curcol,curcol) -= alpha;
889
890 //Copy relevant values of the current column of R_ into the reflection vector
891 //Modify first entry
892 //Take norm of reflection vector
893 //Square the norm
894 v_refl = rcp(new SDV(Teuchos::Copy, &((*R_)(curcol,curcol)), blockSize_ + 1 ));
895 (*v_refl)[0] -= alpha;
896 nvs = blas.NRM2(blockSize_+1,v_refl -> values(),1);
897 nvs *= nvs;
898
899 //
900 //Apply new reflector to:
901 //1. To subsequent columns of R_ in the current block
902 //2. To House[iter_] to store product of reflections for this column
903 //3. To the least-squares right-hand side.
904 //4. The current column
905 //
906 //
907
908
909 //
910 //1.
911 //
912 if(i < blockSize_ - 1){//only do this when there are subsquent columns in the block to apply to
913 workvec = Teuchos::rcp(new SDV(blockSize_ - i -1));
914 //workvec = Teuchos::rcp(new SDV(2*blockSize_));
915 workmatrix = Teuchos::rcp(new SDM (Teuchos::View, *Rblock, blockSize_+1, blockSize_ - i -1, lclCurcol, lclCurcol +1 ) );
916 blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix->values(), workmatrix->stride(), v_refl->values(), 1, zero, workvec->values(), 1);
917 blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl->values(),1,workvec->values(),1,workmatrix->values(),workmatrix->stride());
918 }
919
920
921 //
922 //2.
923 //
924 workvec = Teuchos::rcp(new SDV(2*blockSize_));
925 workmatrix = Teuchos::rcp(new SDM (Teuchos::View, House_[lclIter_ -1], blockSize_+1, 2*blockSize_, i, 0 ) );
926 blas.GEMV(Teuchos::TRANS,workmatrix->numRows(),workmatrix->numCols(),one,workmatrix->values(),workmatrix->stride(), v_refl->values(), 1,zero,workvec->values(),1);
927 blas.GER(workmatrix->numRows(),workmatrix->numCols(), -2.0*one/nvs, v_refl -> values(),1,workvec->values(),1,workmatrix->values(),(*workmatrix).stride());
928
929 //
930 //3.
931 //
932 workvec = Teuchos::rcp(new SDV(blockSize_));
933 workmatrix = Teuchos::rcp(new SDM (Teuchos::View, Z_, blockSize_+1, blockSize_, curcol, 0 ) );
934 blas.GEMV(Teuchos::TRANS, workmatrix->numRows(), workmatrix->numCols(), one, workmatrix-> values(), workmatrix->stride(), v_refl -> values(), 1, zero, workvec->values(), 1);
935 blas.GER((*workmatrix).numRows(),(*workmatrix).numCols(), -2.0*one/nvs,v_refl -> values(), 1,&((*workvec)[0]),1,(*workmatrix)[0],(*workmatrix).stride());
936
937 //
938 //4.
939 //
940 (*R_)[curcol][curcol] = alpha;
941 for(int ii=1; ii<= blockSize_; ii++){
942 (*R_)[curcol][curcol+ii] = 0;
943 }
944 }
945
946 }
947
948 } // end updateLSQR()
949
950
951}//End Belos Namespace
952
953#endif /* BELOS_BLOCK_GCRODR_ITER_HPP */
Belos header file which uses auto-configuration information to include necessary C++ headers.
Class which describes the linear problem to be solved by the iterative solver.
Templated virtual class for providing orthogonalization/orthonormalization methods with matrix-based ...
Declaration of basic traits for the multivector type.
Class which defines basic traits for the operator type.
Class which manages the output and verbosity of the Belos solvers.
Pure virtual base class for defining the status testing capabilities of Belos.
Collection of types and exceptions used within the Belos solvers.
Parent class to all Belos exceptions.
Definition: BelosTypes.hpp:60
BlockGCRODRIterInitFailure is thrown when the BlockGCRODRIter object is unable to generate an initial...
BlockGCRODRIterInitFailure(const std::string &what_arg)
BlockGCRODRIterOrthoFailure is thrown when the BlockGCRODRIter object is unable to compute independen...
BlockGCRODRIterOrthoFailure(const std::string &what_arg)
Implementation of the Block GCRO-DR (Block Recycling GMRES) iteration.
BlockGCRODRIter(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)
BlockGCRODRIter constructor with linear problem, solver utilities, and parameter list of solver optio...
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
const Teuchos::RCP< OrthoManager< ScalarType, MV > > ortho_
Teuchos::RCP< MV > V_
The Krylov basis vectors.
Teuchos::RCP< SDM > B_
Projected matrix from the recycled subspace.
const Teuchos::RCP< OutputManager< ScalarType > > om_
void setSize(int recycledBlocks, int numBlocks)
Set the maximum number of blocks used by the iterative solver and the number of recycled vectors.
bool isInitialized()
States whether the solver has been initialized or not.
Teuchos::SerialDenseVector< int, MagnitudeType > cs_
int getNumBlocks() const
Get the maximum number of blocks used by the iterative solver in solving this linear problem.
SDM Z_
Q applied to right-hand side of the least squares system.
void initialize()
Initialize the solver to an iterate, providing a complete state.
Teuchos::SerialDenseVector< int, ScalarType > SDV
void setRecycledBlocks(int recycledBlocks)
Set the maximum number of recycled blocks used by the iterative solver.
int getRecycledBlocks() const
Set the maximum number of recycled blocks used by the iterative solver.
void setBlockSize(int blockSize)
Set the blocksize.
Teuchos::SerialDenseMatrix< int, ScalarType > SDM
const Teuchos::RCP< StatusTest< ScalarType, MV, OP > > stest_
MultiVecTraits< ScalarType, MV > MVT
BlockGCRODRIterState< ScalarType, MV > getState() const
Get the current state of the linear solver.
const Teuchos::RCP< LinearProblem< ScalarType, MV, OP > > lp_
void setNumBlocks(int numBlocks)
Set the maximum number of blocks used by the iterative solver.
SCT::magnitudeType MagnitudeType
Teuchos::RCP< SDM > R_
Upper triangular reduction of H_ (see above).
Teuchos::ScalarTraits< ScalarType > SCT
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *norms) const
Get the norms of the residuals native to the solver.
virtual ~BlockGCRODRIter()
Destructor.
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.
void iterate()
This method performs block GCRODR iterations until the status test indicates the need to stop or an e...
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::RCP< SDM > H_
Projected matrix from the Krylov factorization.
std::vector< bool > trueRHSIndices_
OperatorTraits< ScalarType, MV, OP > OPT
int getMaxSubspaceDim() const
Get the maximum dimension allocated for the search subspace.
int getNumIters() const
Get the current iteration count.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Teuchos::RCP< MV > U_
Recycled subspace vectors.
A linear system to solve, and its associated information.
Belos's templated virtual class for providing routines for orthogonalization and orthonormzalition of...
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A pure virtual class for defining the status tests for the Belos iterative solvers.
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.