Belos Version of the Day
Loading...
Searching...
No Matches
BelosLSQRIter.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_LSQR_ITER_HPP
43#define BELOS_LSQR_ITER_HPP
44
49#include "BelosConfigDefs.hpp"
50#include "BelosTypes.hpp"
52
55#include "BelosStatusTest.hpp"
58
59#include "Teuchos_SerialDenseMatrix.hpp"
60#include "Teuchos_SerialDenseVector.hpp"
61#include "Teuchos_ScalarTraits.hpp"
62#include "Teuchos_ParameterList.hpp"
63#include "Teuchos_TimeMonitor.hpp"
64
72namespace Belos {
73
74template<class ScalarType, class MV, class OP>
75class LSQRIter : virtual public Belos::Iteration<ScalarType,MV,OP> {
76
77 public:
78
79 //
80 // Convenience typedefs
81 //
84 typedef Teuchos::ScalarTraits<ScalarType> SCT;
85 typedef typename SCT::magnitudeType MagnitudeType;
86
88
89
96 LSQRIter( const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> > &problem,
97 const Teuchos::RCP<Belos::OutputManager<ScalarType> > &printer,
98 const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> > &tester,
99 Teuchos::ParameterList &params );
100
101// If either blocks or reorthogonalization exist, then
102// const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
103
104
106 virtual ~LSQRIter() {};
108
109
111
112
124 void iterate();
125
136
140 {
142 initializeLSQR(empty);
143 }
144
154 state.U = U_; // right Lanczos vector
155 state.V = V_; // left Lanczos vector
156 state.W = W_; // OP * V
157 state.lambda = lambda_;
158 state.resid_norm = resid_norm_;
159 state.frob_mat_norm = frob_mat_norm_;
160 state.mat_cond_num = mat_cond_num_;
161 state.mat_resid_norm = mat_resid_norm_;
162 state.sol_norm = sol_norm_;
163 state.bnorm = bnorm_;
164 return state;
165 }
166
168
169
171
172
174 int getNumIters() const { return iter_; }
175
177 void resetNumIters( int iter = 0 ) { iter_ = iter; }
178
181 Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> * /* norms */ ) const { return Teuchos::null; }
182
184
186 Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
187
189
191
192
195
197 int getBlockSize() const { return 1; }
198
199
201 //This is unique to single vector methods.
202 void setBlockSize(int blockSize) {
203 TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
204 "LSQRIter::setBlockSize(): Cannot use a block size that is not one.");
205 }
206
208 bool isInitialized() { return initialized_; }
209
211
212 private:
213
214 //
215 // Internal methods
216 //
218 void setStateSize();
219
220 //
221 // Classes inputed through constructor that define the linear problem to be solved.
222 //
223 const Teuchos::RCP<Belos::LinearProblem<ScalarType,MV,OP> > lp_;
224 const Teuchos::RCP<Belos::OutputManager<ScalarType> > om_;
225 const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> > stest_;
226// const Teuchos::RCP<OrthoManager<ScalarType,MV> > ortho_;
227
228 //
229 // Current solver state
230 //
231 // initialized_ specifies that the basis vectors have been initialized and the iterate()
232 // routine is capable of running; _initialize is controlled by the initialize() member
233 // method. For the implications of the state of initialized_, please see documentation
234 // for initialize()
235 bool initialized_;
236
237 // stateStorageInitialized_ specifies that the state storage has been initialized.
238 // This initialization may be postponed if the linear problem was generated without
239 // the right-hand side or solution vectors.
240 bool stateStorageInitialized_;
241
242 // Current number of iterations performed.
243 int iter_;
244
245 //
246 // State Storage
247 //
248 //
249 // Bidiagonalization vector
250 Teuchos::RCP<MV> U_;
251 //
252 // Bidiagonalization vector
253 Teuchos::RCP<MV> V_;
254 //
255 // Direction vector
256 Teuchos::RCP<MV> W_;
257 //
258 // Damping value
259 MagnitudeType lambda_;
260 //
261 // Residual norm estimate
262 ScalarType resid_norm_;
263 //
264 // Frobenius norm estimate
265 ScalarType frob_mat_norm_;
266 //
267 // Condition number estimate
268 ScalarType mat_cond_num_;
269 //
270 // A^T*resid norm estimate
271 ScalarType mat_resid_norm_;
272 //
273 // Solution norm estimate
274 ScalarType sol_norm_;
275 //
276 // RHS norm
277 ScalarType bnorm_;
278
279};
280
282 // Constructor.
283
284// const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
285
286 template<class ScalarType, class MV, class OP>
288 const Teuchos::RCP<Belos::OutputManager<ScalarType> > &printer,
289 const Teuchos::RCP<Belos::StatusTest<ScalarType,MV,OP> > &tester,
290 Teuchos::ParameterList &params ):
291 lp_(problem),
292 om_(printer),
293 stest_(tester),
294 initialized_(false),
295 stateStorageInitialized_(false),
296 iter_(0),
297 lambda_(params.get<MagnitudeType> ("Lambda")),
298 resid_norm_(0.0),
299 frob_mat_norm_(0.0),
300 mat_cond_num_(0.0),
301 mat_resid_norm_(0.0),
302 sol_norm_(0.0),
303 bnorm_(0.0)
304 {
305 }
306
308 // Setup the state storage.
309 template <class ScalarType, class MV, class OP>
311 {
312 if (!stateStorageInitialized_) {
313 // Check if there is any multivector to clone from.
314 Teuchos::RCP<const MV> rhsMV = lp_->getInitPrecResVec();
315 Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
316 if (lhsMV == Teuchos::null || rhsMV == Teuchos::null) {
317 stateStorageInitialized_ = false;
318 return;
319 }
320 else {
321 // Initialize the state storage
322 // If the subspace has not been initialized before, generate it
323 // using the LHS and RHS from lp_.
324 if (U_ == Teuchos::null) {
325 // Get the multivectors.
326 TEUCHOS_TEST_FOR_EXCEPTION(rhsMV == Teuchos::null, std::invalid_argument, "LSQRIter::setStateSize(): linear problem does not specify right hand multivector to clone from.");
327 TEUCHOS_TEST_FOR_EXCEPTION(lhsMV == Teuchos::null, std::invalid_argument, "LSQRIter::setStateSize(): linear problem does not specify left hand multivector to clone from.");
328
329 U_ = MVT::Clone( *rhsMV, 1 ); // LeftPrecond * rhs
330 V_ = MVT::Clone( *lhsMV, 1 ); // zero, overwrittein in
331 W_ = MVT::Clone( *lhsMV, 1 ); // zero, initializeLSQR
332 }
333
334 // State storage has now been initialized.
335 stateStorageInitialized_ = true;
336 }
337 }
338 }
339
340
342 // Initialize this iteration object
343 template <class ScalarType, class MV, class OP>
345 {
346 using Teuchos::RCP;
347
348 // Initialize the state storage if it isn't already.
349 if (!stateStorageInitialized_)
350 setStateSize();
351
352 TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
353 "LSQRIter::initialize(): Cannot initialize state storage!");
354
355 std::string errstr("LSQRIter::initialize(): Specified multivectors must have a consistent length and width.");
356
357
358 // Compute initial bidiagonalization vectors and search direction
359 //
360 RCP<const MV> lhsMV = lp_->getLHS(); // contains initial guess,
361
362 const bool debugSerialLSQR = false;
363
364 if( debugSerialLSQR )
365 {
366 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> lhsNorm(1);
367 MVT::MvNorm( *lhsMV, lhsNorm );
368 std::cout << "initializeLSQR lhsNorm " << lhsNorm[0] << std::endl;
369 }
370
371 // LinearProlbem provides right-hand side vectors including RHS CurrRHSVec InitResVec.
372 RCP<const MV> rhsMV = lp_->getInitPrecResVec();
373 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
374 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
375 MVT::MvAddMv( one, *rhsMV, zero, *U_, *U_);
376
377 RCP<const OP> M_left = lp_->getLeftPrec();
378 RCP<const OP> A = lp_->getOperator();
379 RCP<const OP> M_right = lp_->getRightPrec();
380
381 if( debugSerialLSQR )
382 {
383 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> rhsNorm(1);
384 MVT::MvNorm( *U_, rhsNorm );
385 std::cout << "initializeLSQR | U_ | : " << rhsNorm[0] << std::endl;
386 //U_->Print(std::cout);
387 }
388
389 //MVT::MvScale( *V_, zero );
390
391 // Apply the (conjugate) transpose of the preconditioned operator:
392 //
393 // V := (M_L A M_R)^* U, which means
394 // V := M_R^* (A^* (M_L^* U)).
395 //
396 //OPT::Apply(*(lp_->getOperator()), *U_, *V_, CONJTRANS);
397 if ( M_left.is_null())
398 {
399 OPT::Apply (*A, *U_, *V_, CONJTRANS); // V_ = A' U_
400 //std::cout << "*************** V_ ****************" << std::endl;
401 //V_->Print(std::cout);
402 }
403 else
404 {
405 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
406 OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS);
407 OPT::Apply (*A, *tempInRangeOfA, *V_, CONJTRANS); // V_ = A' LeftPrec' U_
408 //std::cout << "mLeft V_ = " << *V_ << std::endl;
409 }
410 if (! M_right.is_null())
411 {
412 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
413
414 OPT::Apply (*M_right, *tempInDomainOfA, *V_, CONJTRANS); // V:= RtPrec' A' LeftPrec' U
415 //std::cout << "mRight V_ = " << *V_ << std::endl;
416 }
417
418 // W := V (copy the vector)
419 MVT::MvAddMv( one, *V_, zero, *V_, *W_);
420
421 frob_mat_norm_ = zero; // These are
422 mat_cond_num_ = one; // lower
423 sol_norm_ = zero; // bounds.
424
425 // The solver is initialized.
426 initialized_ = true;
427 }
428
429
431 // Iterate until the status test informs us we should stop.
432 template <class ScalarType, class MV, class OP>
434 {
435 //
436 // Allocate/initialize data structures
437 //
438 if (initialized_ == false) {
439 initialize();
440 }
441
442 // Create convenience variables for zero and one.
443 const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
444 const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
445 const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
446
447 // Allocate memory for scalars.
448 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> alpha(1);
449 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> beta(1);
450 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> xi(1);
451 // xi is a dumb scalar used for storing inner products.
452 // Eventually SDM will replace the "vectors".
453 std::vector<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> wnorm2(1);
454 ScalarType rhobar, phibar, cs1, phi, rho, cs, sn, theta, xxnorm = MTzero, common;
455 ScalarType zetabar, sn1, psi, res = zero, ddnorm = zero, gamma, tau;
456 ScalarType anorm2 = zero;
457 ScalarType cs2 = -one, sn2 = zero, gammabar, zeta = zero, delta;
458
459 // The pair of work vectors AV and AtU are
460 Teuchos::RCP<MV> AV; // used in applying A to V_ and
461 AV = MVT::Clone( *U_, 1);
462 Teuchos::RCP<MV> AtU; // used in applying A^TRANS to U_ respectively.
463 AtU = MVT::Clone( *V_, 1);
464 const bool debugSerialLSQR = false;
465
466 // Get the current solution vector.
467 Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
468
469
470 // Check that the current solution vector only has one column.
471 TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, LSQRIterateFailure,
472 "LSQRIter::iterate(): current linear system has more than one vector!" );
473
474 // In initializeLSQR among other things V = A' U.
475 // alpha and beta normalize these vectors.
476 MVT::MvNorm( *U_, beta );
477 if( SCT::real(beta[0]) > zero )
478 {
479 MVT::MvScale( *U_, one / beta[0] );
480
481 //std::cout << "*************** U/beta ****************" << std::endl;
482 //U_->Print(std::cout);
483
484 MVT::MvScale( *V_, one / beta[0] ); // scale V = A'U to normalize U
485
486 //std::cout << "*************** V/beta ****************" << std::endl;
487 //V_->Print(std::cout);
488 }
489 MVT::MvNorm( *V_, alpha );
490 if( debugSerialLSQR )
491 {
492 // used to compare with implementations
493 // initializing mat_resid_norm to alpha/beta
494 std::cout << iter_ << " First alpha " << alpha[0] << " beta " << beta[0] << " lambda " << lambda_ << std::endl;
495 }
496 if( SCT::real(alpha[0]) > zero )
497 {
498 MVT::MvScale( *V_, one / alpha[0] ); // V alpha = A' U to normalize V
499 }
500 if( beta[0] * alpha[0] > zero )
501 {
502 MVT::MvScale( *W_, one / (beta[0] * alpha[0]) ); // W = V
503 }
504 else
505 {
506 MVT::MvScale( *W_, zero );
507 }
508
509 using Teuchos::RCP;
510 RCP<const OP> M_left = lp_->getLeftPrec();
511 RCP<const OP> A = lp_->getOperator();
512 RCP<const OP> M_right = lp_->getRightPrec();
513
514 rhobar = alpha[0];
515 phibar = beta[0];
516 bnorm_ = beta[0];
517 resid_norm_ = beta[0];
518 mat_resid_norm_ = alpha[0] * beta[0];
519
520
522 // Iterate until the status test tells us to stop.
523 //
524 while (stest_->checkStatus(this) != Belos::Passed) {
525 // Increment the iteration
526 iter_++;
527
528 // Perform the next step of the bidiagonalization.
529 // The next U_ and V_ vectors and scalars alpha and beta satisfy
530 // U_ betaNew := AV - U_ alphaOld ...
531
532 if ( M_right.is_null() )
533 {
534 OPT::Apply(*A, *V_, *AV); // AV := A * V_
535 }
536 else
537 {
538 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*V_);
539 OPT::Apply (*M_right, *V_, *tempInDomainOfA);
540 OPT::Apply(*A, *tempInDomainOfA, *AV);
541 }
542
543 if (! M_left.is_null())
544 {
545 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*AV);
546 OPT::Apply (*M_left, *tempInRangeOfA, *AV); // AV may change
547 }
548
549
550 if ( !( M_left.is_null() && M_right.is_null() )
551 && debugSerialLSQR && iter_ == 1)
552 {
553 // In practice, LSQR may reveal bugs in transposed preconditioners.
554 // This is the test that catches this type of bug.
555 // 1. confirm that V alpha = A' U
556
557 if (! M_left.is_null())
558 {
559 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
560 OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS);
561 OPT::Apply (*A, *tempInRangeOfA, *AtU, CONJTRANS); // AtU = B'L'U
562 }
563 else
564 {
565 OPT::Apply (*A, *U_, *AtU, CONJTRANS); // AtU = B'U
566 }
567 if ( !( M_right.is_null() ) )
568 {
569 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
570 OPT::Apply (*M_right, *tempInDomainOfA, *AtU, CONJTRANS); // AtU := R' AtU
571 }
572
573 MVT::MvAddMv( one, *AtU, -alpha[0], *V_, *AtU );
574 MVT::MvNorm( *AtU, xi );
575 std::cout << "| V alpha - A' u |= " << xi[0] << std::endl;
576 // 2. confirm that U is a unit vector
577 Teuchos::SerialDenseMatrix<int,ScalarType> uotuo(1,1);
578 MVT::MvTransMv( one, *U_, *U_, uotuo );
579 std::cout << "<U, U> = " << printMat(uotuo) << std::endl;
580 // 3. print alpha = <V, A'U>
581 std::cout << "alpha = " << alpha[0] << std::endl;
582 // 4. compute < AV, U> which ought to be alpha
583 Teuchos::SerialDenseMatrix<int,ScalarType> utav(1,1);
584 MVT::MvTransMv( one, *AV, *U_, utav );
585 std::cout << "<AV, U> = alpha = " << printMat(utav) << std::endl;
586 }
587
588 MVT::MvAddMv( one, *AV, -alpha[0], *U_, *U_ ); // uNew := Av - uOld alphaOld
589 MVT::MvNorm( *U_, beta);
590
591 anorm2 += alpha[0]*alpha[0] + beta[0]*beta[0] + lambda_*lambda_;
592
593
594 if ( SCT::real(beta[0]) > zero )
595 {
596
597 MVT::MvScale( *U_, one / beta[0] );
598
599 if (M_left.is_null())
600 { // ... and V_ alphaNew := AtU - V_ betaNew
601 OPT::Apply(*A, *U_, *AtU, CONJTRANS);
602 }
603 else
604 {
605 RCP<MV> tempInRangeOfA = MVT::CloneCopy (*U_);
606 OPT::Apply (*M_left, *U_, *tempInRangeOfA, CONJTRANS);
607 OPT::Apply(*A, *tempInRangeOfA, *AtU, CONJTRANS);
608 }
609 if (! M_right.is_null())
610 {
611 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*AtU);
612 OPT::Apply (*M_right, *tempInDomainOfA, *AtU, CONJTRANS); // AtU may change
613 }
614
615 MVT::MvAddMv( one, *AtU, -beta[0], *V_, *V_ );
616 MVT::MvNorm( *V_, alpha );
617 }
618 else // beta = 0
619 {
620 alpha[0] = zero;
621 }
622
623 if ( SCT::real(alpha[0]) > zero )
624 {
625 MVT::MvScale( *V_, one / alpha[0] );
626 }
627
628 // Use a plane rotation to eliminate the damping parameter.
629 // This alters the diagonal (rhobar) of the lower-bidiagonal matrix.
630 common = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_);
631 cs1 = rhobar / common;
632 sn1 = lambda_ / common;
633 psi = sn1 * phibar;
634 phibar = cs1 * phibar;
635
636 // Use a plane rotation to eliminate the subdiagonal element (beta)
637 // of the lower-bidiagonal matrix, giving an upper-bidiagonal matrix.
638 rho = Teuchos::ScalarTraits< ScalarType >::squareroot(rhobar*rhobar + lambda_*lambda_ + beta[0]*beta[0]);
639 cs = common / rho;
640 sn = beta[0] / rho;
641 theta = sn * alpha[0];
642 rhobar = -cs * alpha[0];
643 phi = cs * phibar;
644 phibar = sn * phibar; // If beta vanishes, so do sn, theta, phibar and eventually resid_norm
645 tau = sn * phi;
646
647 delta = sn2 * rho;
648 gammabar = -cs2 * rho;
649 zetabar = (phi - delta*zeta) / gammabar;
650 sol_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(xxnorm + zetabar*zetabar);
651 gamma = Teuchos::ScalarTraits< ScalarType >::squareroot(gammabar*gammabar + theta*theta);
652 cs2 = gammabar / gamma;
653 sn2 = theta / gamma;
654 zeta = (phi - delta*zeta) / gamma;
655 xxnorm += zeta*zeta;
656
657 //The next task may be addressed by some form of lp_->updateSolution.
658 if ( M_right.is_null())
659 {
660 MVT::MvAddMv( phi / rho, *W_, one, *cur_soln_vec, *cur_soln_vec);
661 }
662 else
663 {
664 RCP<MV> tempInDomainOfA = MVT::CloneCopy (*W_);
665 OPT::Apply (*M_right, *W_, *tempInDomainOfA);
666 MVT::MvAddMv( phi / rho, *tempInDomainOfA, one, *cur_soln_vec, *cur_soln_vec);
667 }
668
669 MVT::MvNorm( *W_, wnorm2 );
670 ddnorm += (one / rho)*(one / rho) * wnorm2[0]*wnorm2[0];
671 MVT::MvAddMv( one, *V_, -theta / rho, *W_, *W_ );
672
673 frob_mat_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(anorm2);
674 mat_cond_num_ = frob_mat_norm_ * Teuchos::ScalarTraits< ScalarType >::squareroot(ddnorm);
675 res+= psi*psi;
676 resid_norm_ = Teuchos::ScalarTraits< ScalarType >::squareroot(phibar*phibar + res);
677 mat_resid_norm_ = alpha[0] * Teuchos::ScalarTraits< ScalarType >::magnitude(tau);
678
679 } // end while (sTest_->checkStatus(this) != Passed)
680 } // iterate()
681
682} // end Belos namespace
683
684#endif /* BELOS_LSQR_ITER_HPP */
Belos header file which uses auto-configuration information to include necessary C++ headers.
IterationState contains the data that defines the state of the LSQR solver at any given time.
Class which describes the linear problem to be solved by the iterative solver.
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.
Implementation of the LSQR iteration.
int getBlockSize() const
Get the blocksize to be used by the iterative solver in solving this linear problem.
virtual ~LSQRIter()
Destructor.
LSQRIter(const Teuchos::RCP< Belos::LinearProblem< ScalarType, MV, OP > > &problem, const Teuchos::RCP< Belos::OutputManager< ScalarType > > &printer, const Teuchos::RCP< Belos::StatusTest< ScalarType, MV, OP > > &tester, Teuchos::ParameterList &params)
LSQRIter constructor with linear problem, solver utilities, and parameter list of solver options.
bool isInitialized()
States whether the solver has been initialized or not.
void iterate()
This method performs LSQR iterations until the status test indicates the need to stop or an error occ...
Teuchos::RCP< const MV > getNativeResiduals(std::vector< MagnitudeType > *) const
Get the norms of the residuals native to the solver.
LSQRIterationState< ScalarType, MV > getState() const
Get the current state of the linear solver.
Teuchos::RCP< MV > getCurrentUpdate() const
Get the current update to the linear system.
Belos::OperatorTraits< ScalarType, MV, OP > OPT
int getNumIters() const
Get the current iteration count.
void resetNumIters(int iter=0)
Reset the iteration count.
Teuchos::ScalarTraits< ScalarType > SCT
Belos::MultiVecTraits< ScalarType, MV > MVT
void initializeLSQR(LSQRIterationState< ScalarType, MV > &newstate)
Initialize the solver to an iterate, completing the initial state.
SCT::magnitudeType MagnitudeType
void setBlockSize(int blockSize)
Set the blocksize to be used by the iterative solver to solve this linear problem.
void initialize()
The solver is initialized using initializeLSQR.
const Belos::LinearProblem< ScalarType, MV, OP > & getProblem() const
Get a constant reference to the linear problem.
LSQRIterateFailure is thrown when the LSQRIteration object is unable to compute the next iterate in t...
A linear system to solve, and its associated information.
Traits class which defines basic operations on multivectors.
Class which defines basic traits for the operator type.
Belos's basic output manager for sending information of select verbosity levels to the appropriate ou...
A pure virtual class for defining the status tests for the Belos iterative solvers.
@ CONJTRANS
Definition: BelosTypes.hpp:83
Structure to contain pointers to LSQRIteration state variables, ...
ScalarType sol_norm
An estimate of the norm of the solution.
ScalarType frob_mat_norm
An approximation to the Frobenius norm of A.
Teuchos::RCP< const MV > V
Bidiagonalization vector.
ScalarType mat_resid_norm
An estimate of the norm of A^T*resid.
Teuchos::RCP< const MV > U
Bidiagonalization vector.
Teuchos::RCP< const MV > W
The search direction vector.
Teuchos::ScalarTraits< ScalarType >::magnitudeType lambda
The damping value.
ScalarType bnorm
The norm of the RHS vector b.
ScalarType resid_norm
The current residual norm.
ScalarType mat_cond_num
An approximation to the condition number of A.

Generated for Belos by doxygen 1.9.6