lsmrModule Module

LSMR solves Ax = b or min ||Ax - b|| with or without damping, using the iterative algorithm of David Fong and Michael Saunders.

Authors

  • David Fong & Michael Saunders, Systems Optimization Laboratory (SOL)

See also

  • http://www.stanford.edu/group/SOL/software/lsmr.html

History

  • 17 Jul 2010: F90 LSMR derived from F90 LSQR and lsqr.m.
  • 07 Sep 2010: Local reorthogonalization now works (localSize > 0).
  • 28 Jan 2014: In lsmrDataModule.f90: ip added for integer(ip) declarations. dnrm2 and dscal coded directly (no longer use lsmrblasInterface.f90 or lsmrblas.f90).

Uses

  • module~~lsmrmodule~~UsesGraph module~lsmrmodule lsmrModule module~lsmrdatamodule lsmrDataModule module~lsmrmodule->module~lsmrdatamodule iso_fortran_env iso_fortran_env module~lsmrdatamodule->iso_fortran_env

Abstract Interfaces

abstract interface

  • private subroutine Aprod1_f(m, n, x, y)

    y := y + A*x

    Arguments

    Type IntentOptional Attributes Name
    integer(kind=ip), intent(in) :: m
    integer(kind=ip), intent(in) :: n
    real(kind=dp), intent(in) :: x(n)
    real(kind=dp), intent(inout) :: y(m)

abstract interface

  • private subroutine Aprod2_f(m, n, x, y)

    x := x + A'*y

    Arguments

    Type IntentOptional Attributes Name
    integer(kind=ip), intent(in) :: m
    integer(kind=ip), intent(in) :: n
    real(kind=dp), intent(inout) :: x(n)
    real(kind=dp), intent(in) :: y(m)

Subroutines

public subroutine lsmr(m, n, Aprod1, Aprod2, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx)

LSMR finds a solution x to the following problems:

Read more…

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: m

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

the number of columns in A.
procedure(Aprod1_f) :: Aprod1

See above.
procedure(Aprod2_f) :: Aprod2

See above.
real(kind=dp), intent(in) :: b(m)

The rhs vector b.
real(kind=dp), 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(eps)*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 LSMR.

Read more…
real(kind=dp), 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=dp), intent(in) :: btol

An estimate of the relative error in the data defining the rhs b. For example, if b is accurate to about 6 digits, set btol = 1.0e-6.
real(kind=dp), 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(kind=ip), intent(in) :: itnlim

An upper limit on the number of iterations. Suggested value:

Read more…
integer(kind=ip), intent(in) :: localSize

No. of vectors for local reorthogonalization:

Read more…
integer(kind=ip), intent(in) :: nout

File number for printed output. If positive, a summary will be printed on file nout.
real(kind=dp), intent(out) :: x(n)

Returns the computed solution x.
integer(kind=ip), intent(out) :: istop

An integer giving the reason for termination:

Read more…
integer(kind=ip), intent(out) :: itn

The number of iterations performed.
real(kind=dp), intent(out) :: normA

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 the columns of A have all been scaled to have length 1.0, normA should increase to roughly sqrt(n). A radically different value for normA may indicate an error in Aprod1 or Aprod2.
real(kind=dp), intent(out) :: condA

An estimate of cond(Abar), the condition number of Abar. A very high value of condA may again indicate an error in Aprod1 or Aprod2.
real(kind=dp), intent(out) :: normr

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=dp), intent(out) :: normAr

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

An estimate of norm(x) for the final solution x.

public subroutine lsmr_ez(m, n, irow, icol, a, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx)

Easy interface to lsmr.

Read more…

Arguments

Type IntentOptional Attributes Name
integer(kind=ip), intent(in) :: m
integer(kind=ip), intent(in) :: n
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=dp), intent(in), dimension(:) :: a

nonzero elements of A
real(kind=dp), intent(in) :: b(m)
real(kind=dp), intent(in) :: damp
real(kind=dp), intent(in) :: atol
real(kind=dp), intent(in) :: btol
real(kind=dp), intent(in) :: conlim
integer(kind=ip), intent(in) :: itnlim
integer(kind=ip), intent(in) :: localSize
integer(kind=ip), intent(in) :: nout
real(kind=dp), intent(out) :: x(n)
integer(kind=ip), intent(out) :: istop
integer(kind=ip), intent(out) :: itn
real(kind=dp), intent(out) :: normA
real(kind=dp), intent(out) :: condA
real(kind=dp), intent(out) :: normr
real(kind=dp), intent(out) :: normAr
real(kind=dp), intent(out) :: normx