lsqr_solver Derived Type

type, public :: lsqr_solver

main class to access the lsqr solver.

You can use this class directory by extending it and specifying aprod, or you can use the lsqr_solver_ez class that has an easier interface.


Inherited by

type~~lsqr_solver~~InheritedByGraph type~lsqr_solver lsqr_solver type~lsqr_solver_ez lsqr_solver_ez type~lsqr_solver_ez->type~lsqr_solver

Type-Bound Procedures

procedure(aprod_func), public, deferred :: aprod

User function to access the sparse matrix A.

  • subroutine aprod_func(me, mode, m, n, x, y) Prototype

    User function to access the sparse matrix A.

    Arguments

    Type IntentOptional Attributes Name
    class(lsqr_solver), intent(inout) :: me
    integer, intent(in) :: mode
    • If mode = 1, compute y = y + A*x. y should be altered without changing x.
    • If mode = 2, compute x = x + A(transpose)*y. x should be altered without changing y.
    integer, intent(in) :: m

    number of rows in A matrix
    integer, intent(in) :: n

    number of columns in A matrix
    real(kind=wp), intent(inout), dimension(:) :: x
    real(kind=wp), intent(inout), dimension(:) :: y

procedure, public :: lsqr => LSQR

main solver routine

  • private subroutine LSQR(me, m, n, damp, wantse, u, v, w, x, se, atol, btol, conlim, itnlim, nout, istop, itn, anorm, acond, rnorm, arnorm, xnorm)

    LSQR finds a solution to the following problems:

    Read more…

    Arguments

    Type IntentOptional Attributes Name
    class(lsqr_solver), intent(inout) :: me
    integer, intent(in) :: m

    the number of rows in A.
    integer, intent(in) :: n

    the number of columns in A.
    real(kind=wp), intent(in) :: damp

    The damping parameter for problem 3 above. (damp should be 0.0 for problems 1 and 2.) If the system A*x = b is incompatible, values of damp in the range 0 to sqrt(relpr)*norm(A) will probably have a negligible effect. Larger values of damp will tend to decrease the norm of x and reduce the number of iterations required by LSQR.

    Read more…
    logical, intent(in) :: wantse

    A logical variable to say if the array se(*) of standard error estimates should be computed. If m > n or damp > 0, the system is overdetermined and the standard errors may be useful. (See the first LSQR reference.) Otherwise (m <= n and damp = 0) they do not mean much. Some time and storage can be saved by setting wantse = .false. and using any convenient array for se(*), which won't be touched.
    real(kind=wp), intent(inout) :: u(m)

    The rhs vector b. Beware that u is over-written by LSQR.
    real(kind=wp), intent(inout) :: v(n)

    workspace
    real(kind=wp), intent(inout) :: w(n)

    workspace
    real(kind=wp), intent(out) :: x(n)

    Returns the computed solution x.
    real(kind=wp), intent(out), dimension(*) :: se

    If wantse is true, the dimension of se must be n or more. se(*) then returns standard error estimates for the components of x. For each i, se(i) is set to the value rnorm * sqrt( sigma(i,i) / t ), where sigma(i,i) is an estimate of the i-th diagonal of the inverse of Abar(transpose)*Abar and: * t = 1 if m <= n * t = m - n if m > n and damp = 0 * t = m if damp /= 0

    Read more…
    real(kind=wp), intent(in) :: atol

    An estimate of the relative error in the data defining the matrix A. For example, if A is accurate to about 6 digits, set atol = 1.0e-6.
    real(kind=wp), intent(in) :: btol

    An estimate of the relative error in the data defining the rhs vector b. For example, if b is accurate to about 6 digits, set btol = 1.0e-6.
    real(kind=wp), intent(in) :: conlim

    An upper limit on cond(Abar), the apparent condition number of the matrix Abar. Iterations will be terminated if a computed estimate of cond(Abar) exceeds conlim. This is intended to prevent certain small or zero singular values of A or Abar from coming into effect and causing unwanted growth in the computed solution.

    Read more…
    integer, intent(in) :: itnlim

    An upper limit on the number of iterations. Suggested value: * itnlim = n/2 for well-conditioned systems with clustered singular values, * itnlim = 4*n otherwise.
    integer, intent(in) :: nout

    File number for printed output. If nonzero, a summary will be printed on file nout.
    integer, intent(out) :: istop

    An integer giving the reason for termination:

    Read more…
    integer, intent(out) :: itn

    The number of iterations performed.
    real(kind=wp), intent(out) :: anorm

    An estimate of the Frobenius norm of Abar. This is the square-root of the sum of squares of the elements of Abar. If damp is small and if the columns of A have all been scaled to have length 1.0, anorm should increase to roughly sqrt(n). A radically different value for anorm may indicate an error in subroutine aprod (there may be an inconsistency between modes 1 and 2).
    real(kind=wp), intent(out) :: acond

    An estimate of cond(Abar), the condition number of Abar. A very high value of acond may again indicate an error in aprod.
    real(kind=wp), intent(out) :: rnorm

    An estimate of the final value of norm(rbar), the function being minimized (see notation above). This will be small if A*x = b has a solution.
    real(kind=wp), intent(out) :: arnorm

    An estimate of the final value of norm( Abar(transpose)*rbar ), the norm of the residual for the usual normal equations. This should be small in all cases. (arnorm will often be smaller than the true value computed from the output vector x.)
    real(kind=wp), intent(out) :: xnorm

    An estimate of the norm of the final solution vector x.

procedure, public :: acheck

  • private subroutine acheck(me, m, n, nout, eps, v, w, x, y, inform)

    Checks the two modes of aprod for LSQR.

    Read more…

    Arguments

    Type IntentOptional Attributes Name
    class(lsqr_solver), intent(inout) :: me
    integer, intent(in) :: m

    No. of rows of A.
    integer, intent(in) :: n

    No. of columns of A.
    integer, intent(in) :: nout

    A file number for printed output.
    real(kind=wp), intent(in) :: eps

    The machine precision.
    real(kind=wp) :: v(n)
    real(kind=wp) :: w(m)
    real(kind=wp) :: x(n)
    real(kind=wp) :: y(m)
    integer, intent(out) :: inform

    Error indicator. inform = 0 if aprod seems to be consistent. inform = 1 otherwise.

procedure, public :: xcheck

  • private subroutine xcheck(me, m, n, nout, anorm, damp, eps, b, u, v, w, x, inform, test1, test2, test3)

    Tests if x solves a certain least-squares problem.

    Read more…

    Arguments

    Type IntentOptional Attributes Name
    class(lsqr_solver), intent(inout) :: me
    integer, intent(in) :: m

    The number of rows in A.
    integer, intent(in) :: n

    The number of columns in A.
    integer, intent(in) :: nout

    A file number for printed output. If nout = 0, nothing is printed.
    real(kind=wp), intent(in) :: anorm

    An estimate of norm(A) or norm( A, delta*I ) if delta > 0. Normally this will be available from LSQR or CRAIG.
    real(kind=wp), intent(in) :: damp

    Possibly defines a damped problem.
    real(kind=wp), intent(in) :: eps

    Machine precision.
    real(kind=wp), intent(in) :: b(m)

    The right-hand side of Ax = b etc.
    real(kind=wp), intent(out) :: u(m)

    On exit, u = r (where r = b - Ax).
    real(kind=wp), intent(out) :: v(n)

    On exit, v = A'r.
    real(kind=wp), intent(out) :: w(n)

    On exit, w = A'r - damp^2 x.
    real(kind=wp), intent(in) :: x(n)

    The given estimate of a solution.
    integer, intent(out) :: inform

    inform = 0 if b = 0 and x = 0. inform = 1, 2 or 3 if x seems to solve systems 1 2 or 3 above.
    real(kind=wp), intent(out) :: test1

    These are dimensionless quantities that should be "small" if x does seem to solve one of the systems. "small" means less than tol = eps**power, where power is defined as a parameter below.
    real(kind=wp), intent(out) :: test2

    These are dimensionless quantities that should be "small" if x does seem to solve one of the systems. "small" means less than tol = eps**power, where power is defined as a parameter below.
    real(kind=wp), intent(out) :: test3

    These are dimensionless quantities that should be "small" if x does seem to solve one of the systems. "small" means less than tol = eps**power, where power is defined as a parameter below.

Source Code

   type,abstract,public :: lsqr_solver
      !! main class to access the [[lsqr]] solver.
      !!
      !! You can use this class directory by extending it
      !! and specifying `aprod`, or you can use the
      !! [[lsqr_solver_ez]] class that has an easier
      !! interface.
      private
   contains
      private
      procedure(aprod_func),deferred,public :: aprod !! User function to access the sparse matrix `A`.
      procedure,public :: lsqr  !! main solver routine
      procedure,public :: acheck
      procedure,public :: xcheck
   end type lsqr_solver