lsqr_solver_ez

type, public, extends(lsqr_solver) :: lsqr_solver_ez

a simplier version of lsqr_solver where the aprod function is provided internally. To use, first call the initialize method to set the matrix and other inputs.

Inherits

Help

Components

Type	Visibility	Attributes		Name		Initial
integer,	private		::	m	=	0	number of rows in `A` matrix
integer,	private		::	n	=	0	number of columns in `A` matrix
integer,	private		::	num_nonzero_elements	=	0	number of nonzero elements in `A` matrix
integer,	private,	dimension(:), allocatable	::	irow			sparsity row indices
integer,	private,	dimension(:), allocatable	::	icol			sparsity column indices
real(kind=wp),	private,	dimension(:), allocatable	::	a			sparse `A` matrix
real(kind=wp),	private		::	atol	=	zero	relative error in definition of `A`
real(kind=wp),	private		::	btol	=	zero	relative error in definition of `b`
real(kind=wp),	private		::	conlim	=	zero	An upper limit on `cond(Abar)`, the apparent condition number of the matrix `Abar`.
integer,	private		::	itnlim	=	100	max iterations
integer,	private		::	nout	=	0	output unit for printing
real(kind=wp),	private,	dimension(:), allocatable	::	Ax			`A*x` (dimension m)
real(kind=wp),	private,	dimension(:), allocatable	::	Aty			`A(transpose)*y` (dimension n)
real(kind=wp),	private,	dimension(:), allocatable	::	v			workspace array (dimension n)
real(kind=wp),	private,	dimension(:), allocatable	::	w			workspace array (dimension n)

Type-Bound Procedures

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 $x$ to the following problems:

Arguments

Type	Intent	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 `Ax = 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.

Arguments

Type	Intent		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.

Arguments

Type	Intent		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.

procedure, public :: initialize => initialize_ez

Constructor. Must be call first.

private subroutine initialize_ez(me, m, n, a, irow, icol, atol, btol, conlim, itnlim, nout)

Constructor for lsqr_solver_ez.

Arguments

Type	Intent	Optional	Attributes		Name
class(lsqr_solver_ez),	intent(out)			::	me
integer,	intent(in)			::	m	number of rows in `A` matrix
integer,	intent(in)			::	n	number of columns in `A` matrix
real(kind=wp),	intent(in),		dimension(:)	::	a	nonzero elements of `A`
integer,	intent(in),		dimension(:)	::	irow	row indices of nonzero elements of `A`
integer,	intent(in),		dimension(:)	::	icol	column indices of nonzero elements of `A`
real(kind=wp),	intent(in),	optional		::	atol	relative error in definition of `A`
real(kind=wp),	intent(in),	optional		::	btol	relative error in definition of `b`
real(kind=wp),	intent(in),	optional		::	conlim	An upper limit on `cond(Abar)`, the apparent condition number of the matrix `Abar`.
integer,	intent(in),	optional		::	itnlim	max iterations
integer,	intent(in),	optional		::	nout	output unit for printing

procedure, public :: solve => solve_ez

private subroutine solve_ez(me, b, damp, x, istop, se, itn, anorm, acond, rnorm, arnorm, xnorm)

Wrapper for LSQR for the easy version of the class.

Arguments

Type	Intent	Optional	Attributes		Name
class(lsqr_solver_ez),	intent(inout)			::	me
real(kind=wp),	intent(in),		dimension(me%m)	::	b
real(kind=wp),	intent(in)			::	damp
real(kind=wp),	intent(out),		dimension(me%n)	::	x	the computed solution `x`.
integer,	intent(out)			::	istop	exit code (see LSQR).
real(kind=wp),	intent(out),	optional,	dimension(me%n)	::	se
integer,	intent(out),	optional		::	itn
real(kind=wp),	intent(out),	optional		::	anorm
real(kind=wp),	intent(out),	optional		::	acond
real(kind=wp),	intent(out),	optional		::	rnorm
real(kind=wp),	intent(out),	optional		::	arnorm
real(kind=wp),	intent(out),	optional		::	xnorm

procedure, public :: aprod => aprod_ez

internal routine

private subroutine aprod_ez(me, mode, m, n, x, y)

The internal aprod function for the lsqr_solver_ez class.

Arguments

Type	Intent	Attributes		Name
class(lsqr_solver_ez),	intent(inout)		::	me
integer,	intent(in)		::	mode	Read more…
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	[n]
real(kind=wp),	intent(inout),	dimension(:)	::	y	[m]

Source Code

   type,public,extends(lsqr_solver) :: lsqr_solver_ez
      !! a simplier version of [[lsqr_solver]] where
      !! the `aprod` function is provided internally.
      !! To use, first call the `initialize` method
      !! to set the matrix and other inputs.
      private

      integer :: m = 0 !! number of rows in `A` matrix
      integer :: n = 0 !! number of columns in `A` matrix
      integer :: num_nonzero_elements = 0 !! number of nonzero elements in `A` matrix
      integer,dimension(:),allocatable  :: irow !! sparsity row indices
      integer,dimension(:),allocatable  :: icol !! sparsity column indices
      real(wp),dimension(:),allocatable :: a    !! sparse `A` matrix

      real(wp) :: atol    = zero  !! relative error in definition of `A`
      real(wp) :: btol    = zero  !! relative error in definition of `b`
      real(wp) :: conlim  = zero  !! An upper limit on `cond(Abar)`, the apparent
                                  !! condition number of the matrix `Abar`.
      integer  :: itnlim  = 100   !! max iterations
      integer  :: nout    = 0     !! output unit for printing

      ! used in aprod_ez:
      real(wp),dimension(:),allocatable :: Ax   !! `A*x` (dimension m)
      real(wp),dimension(:),allocatable :: Aty  !! `A(transpose)*y` (dimension n)

      real(wp),dimension(:),allocatable :: v  !! workspace array (dimension n)
      real(wp),dimension(:),allocatable :: w  !! workspace array (dimension n)

   contains
      private
      procedure,public :: initialize => initialize_ez  !! Constructor. Must be call first.
      procedure,public :: solve => solve_ez
      procedure,public :: aprod => aprod_ez !! internal routine
   end type lsqr_solver_ez

lsqr_solver_ez Derived Type

Contents

Variables

Type-Bound Procedures

Source Code

type, public, extends(lsqr_solver) :: lsqr_solver_ez

Inherits

Components

Type-Bound Procedures

procedure, public :: lsqr => LSQR

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)

Arguments

procedure, public :: acheck

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

Arguments

procedure, public :: xcheck

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

Arguments

procedure, public :: initialize => initialize_ez

private subroutine initialize_ez(me, m, n, a, irow, icol, atol, btol, conlim, itnlim, nout)

Arguments

procedure, public :: solve => solve_ez

private subroutine solve_ez(me, b, damp, x, istop, se, itn, anorm, acond, rnorm, arnorm, xnorm)

Arguments

procedure, public :: aprod => aprod_ez

private subroutine aprod_ez(me, mode, m, n, x, y)

Arguments

Source Code