Skip to content

Commit

Permalink
Merge pull request #17 from jacobwilliams/16-lsmr
Browse files Browse the repository at this point in the history 
added the LSMR solver as an option
  • Loading branch information
@jacobwilliams
jacobwilliams authored Jan 27, 2024
2 parents 3367ca7 + 6d34f47 commit 0584391
Show file tree
Hide file tree
Showing 5 changed files with 76 additions and 11 deletions.
3 changes: 2 additions & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@ A work in progress.
1. [LAPACK](https://www.netlib.org/lapack/) routines (`dgesv` or `dgels`) for dense systems: If `n=m`, uses `dgesv` (LU decomposition). If `n/=m`, uses `dgels` (if `m>n` uses QR factorization, if `m<n` uses LQ factorization).
2. [lsqr](https://github.com/jacobwilliams/LSQR) -- a conjugate-gradient type method for solving sparse linear equations and sparse least-squares problems.
3. [lusol](https://github.com/jacobwilliams/lusol) -- A sparse LU factorization for square and rectangular matrices, with Bartels-Golub-Reid updates for column replacement and other rank-1 modifications.
4. [lsmr](https://github.com/jacobwilliams/LSMR) -- a conjugate-gradient type method for solving sparse linear equations and sparse least-squares problems
* Has a Broyden update option.
* Has various line search options.
* use a specified constant step size (0,1]
Expand Down Expand Up @@ -66,7 +67,7 @@ Or to use a specific version:
nlesolver-fortran = { git="https://github.com/jacobwilliams/nlesolver-fortran.git", tag="1.1.0" }
```

Note that LAPACK is required to build. The [fmin](https://github.com/jacobwilliams/fmin), [lsqr](https://github.com/jacobwilliams/LSQR), and [lusol](https://github.com/jacobwilliams/lusol) libraries are also dependencies (which will be automatically downloaded by fpm).
Note that LAPACK is required to build. The [fmin](https://github.com/jacobwilliams/fmin), [lsqr](https://github.com/jacobwilliams/LSQR), [lusol](https://github.com/jacobwilliams/lusol), and [lsmr](https://github.com/jacobwilliams/LSMR) libraries are also dependencies (which will be automatically downloaded by fpm).

### Documentation

Expand Down
2 changes: 1 addition & 1 deletion codecov.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ coverage:
status:
patch:
default:
target: 70%
target: 40%
project:
default:
target: 60%
1 change: 1 addition & 0 deletions fpm.toml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@ source-dir = "src"
fmin = { git="https://github.com/jacobwilliams/fmin.git", tag="1.1.1" }
LSQR = { git="https://github.com/jacobwilliams/LSQR", tag="1.1.0" }
lusol = { git="https://github.com/jacobwilliams/lusol", tag="1.0.0" }
LSMR = { git="https://github.com/jacobwilliams/LSMR" }

[install]
library = true
51 changes: 50 additions & 1 deletion src/nlesolver_module.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@ module nlesolver_module
use fmin_module, only: fmin
use lsqr_module, only: lsqr_solver_ez
use lusol_ez_module, only: solve, lusol_settings
use lsmrModule, only: lsmr_ez

implicit none

Expand Down Expand Up @@ -109,6 +110,7 @@ module nlesolver_module
!! * 1 - assume dense (use dense solver)
!! * 2 - assume sparse (use LSQR sparse solver).
!! * 3 - assume sparse (use LUSOL sparse solver).
!! * 4 - assume sparse (use LSMR sparse solver).
integer :: n_nonzeros = -1 !! number of nonzero Jacobian elements (used for `sparsity_mode > 1`)
integer,dimension(:),allocatable :: irow !! sparsity pattern nonzero elements row indices.
integer,dimension(:),allocatable :: icol !! sparsity pattern nonzero elements column indices
Expand All @@ -122,6 +124,15 @@ module nlesolver_module
integer :: nout = 0 !! output unit for printing
real(wp) :: damp = zero !! damp parameter for LSQR

integer :: localSize = 0 !! LSMR: Number of vectors for local reorthogonalization:
!!
!! * 0 No reorthogonalization is performed.
!! * >0 This many n-vectors "v" (the most recent ones)
!! are saved for reorthogonalizing the next v.
!!
!! localSize need not be more than `min(m,n)`.
!! At most `min(m,n)` vectors will be allocated.

! LUSOL parameters:
integer :: lusol_method = 0 !! * 0 => TPP: Threshold Partial Pivoting.
!! * 1 => TRP: Threshold Rook Pivoting.
Expand Down Expand Up @@ -321,7 +332,7 @@ subroutine initialize_nlesolver_variables(me,&
verbose,iunit,n_uphill_max,n_intervals,&
sparsity_mode,irow,icol,&
atol,btol,conlim,damp,itnlim,nout,&
lusol_method )
lusol_method,localSize )

implicit none

Expand Down Expand Up @@ -366,6 +377,7 @@ subroutine initialize_nlesolver_variables(me,&
!! Must specify `grad` for this mode.
!! * 2 - assume sparse (use LSQR sparse solver).
!! * 3 - assume sparse (use LUSOL sparse solver).
!! * 4 - assume sparse (use LSMR sparse solver).
!! Must also specify `grad_sparse` and the `irow` and `icol`
!! sparsity pattern for `mode>1`.
integer,dimension(:),intent(in),optional :: irow !! sparsity pattern nonzero elements row indices.
Expand All @@ -384,6 +396,14 @@ subroutine initialize_nlesolver_variables(me,&
!! * 0 => TPP: Threshold Partial Pivoting.
!! * 1 => TRP: Threshold Rook Pivoting.
!! * 2 => TCP: Threshold Complete Pivoting.
integer,intent(in),optional :: localSize !! `LSMR`: Number of vectors for local reorthogonalization:
!!
!! * 0 No reorthogonalization is performed.
!! * >0 This many n-vectors "v" (the most recent ones)
!! are saved for reorthogonalizing the next v.
!!
!! localSize need not be more than `min(m,n)`.
!! At most `min(m,n)` vectors will be allocated.

logical :: status_ok !! true if there were no errors

Expand Down Expand Up @@ -500,6 +520,9 @@ subroutine initialize_nlesolver_variables(me,&
! LUSOL method
if (present(nout)) me%lusol_method = lusol_method

! LSMR optional inputs:
if (present(localSize)) me%localSize = localSize

end if
end if

Expand Down Expand Up @@ -549,6 +572,8 @@ subroutine nlesolver_solver(me,x)
integer :: i !! counter
integer,dimension(:),allocatable :: idx, index_array !! for sparse indexing
character(len=10) :: i_str !! string version of `i` for row string
real(wp) :: normA, condA, normr, normAr, normx !! for LSMR
integer :: itn !! for LSMR

if (me%istat<0) return ! class was not initialized properly

Expand Down Expand Up @@ -746,6 +771,30 @@ subroutine nlesolver_solver(me,x)
lusol_options%method = me%lusol_method
call solve(me%n,me%m,me%n_nonzeros,me%irow,me%icol,fjac_sparse,rhs,p,info,&
settings=lusol_options)
case (4)

! use LSMR solver:
!me%conlim = 1.0_wp/(100*epsilon(1.0_wp))
call lsmr_ez ( me%m, me%n, me%irow, me%icol, fjac_sparse, rhs, me%damp, &
me%atol, me%btol, me%conlim, me%itnlim, me%localSize, me%nout, &
p, info, itn, normA, condA, normr, normAr, normx )

! check convergence:
select case (info)
case(4)
call me%set_status(istat = -1004, &
string = 'LSMR Error: The system appears to be ill-conditioned. istop =', i=info)
exit
case(5)
call me%set_status(istat = -1005, &
string = 'LSMR Error: The iteration limit was reached. istop =', i=info)
exit
case default
info = 0
end select

case default
error stop 'invalid sparsity_mode'
end select

! check for errors:
Expand Down
30 changes: 22 additions & 8 deletions test/sparse_test.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,22 @@ program sparse_test
integer,parameter :: max_iter = 100
real(wp),parameter :: tol = 1.0e-8_wp
logical,parameter :: verbose = .false.
integer,dimension(3),parameter :: icol = [1,2,2]
integer,dimension(3),parameter :: irow = [1,1,2]

integer :: f_evals

call go(2, 'LSQR')
call go(3, 'LUSOL')
call go(4, 'LSMR')

contains

subroutine go(sparsity_mode, mode_name)
!! run the tests for the specified sparsity mode

integer,intent(in) :: sparsity_mode
character(len=*),intent(in) :: mode_name !! name of the sparsity mode used

type(nlesolver_type) :: solver
real(wp) :: alpha
Expand All @@ -22,21 +38,17 @@ program sparse_test
integer :: istat !! Integer status code.
character(len=:),allocatable :: message !! Text status message
real(wp),dimension(n) :: x
integer :: f_evals
integer :: i
character(len=:),allocatable :: description
real(wp) :: fmin_tol

integer,dimension(3),parameter :: icol = [1,2,2]
integer,dimension(3),parameter :: irow = [1,1,2]

fmin_tol = 1.0e-2_wp ! don't need a tight tol for this
n_intervals = 2
alpha = 1.0_wp

write(*,*) ''
write(*,*) '***********************'
write(*,*) '* sparse_test *'
write(*,*) '* sparse_test : '//trim(mode_name)
write(*,*) '***********************'
write(*,*) ''
do i = 1, 8
Expand Down Expand Up @@ -94,6 +106,8 @@ program sparse_test
m = m, &
max_iter = max_iter, &
tol = tol, &
atol = tol, &
btol = tol, &
func = func, &
grad_sparse = grad_sparse, &
step_mode = step_mode,&
Expand All @@ -102,10 +116,10 @@ program sparse_test
n_intervals = n_intervals, &
fmin_tol = fmin_tol, &
verbose = verbose,&
sparsity_mode = 3,&
sparsity_mode = sparsity_mode,& ! lsmr
irow = irow,&
icol = icol,&
damp = 1.0_wp)
damp = 0.0_wp)
call solver%status(istat, message)
write(*,'(I3,1X,A)') istat, message
if (istat /= 0) error stop
Expand All @@ -119,7 +133,7 @@ program sparse_test

end do

contains
end subroutine go

subroutine func(me,x,f)
!! compute the function
Expand Down

0 comments on commit 0584391

Please sign in to comment.