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

Generated on Fri Aug 14 2020 10:48:33 for Belos by doxygen 1.8.5