Skip to content

Commit

Permalink
added the LSMR solver as an option
Browse files Browse the repository at this point in the history 
Experimental
See #16
  • Loading branch information
@jacobwilliams
jacobwilliams committed Jan 25, 2024
1 parent 3367ca7 commit 89ffc2f
Show file tree
Hide file tree
Showing 3 changed files with 55 additions and 8 deletions.
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
32 changes: 32 additions & 0 deletions 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 Down Expand Up @@ -366,6 +368,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 Down Expand Up @@ -549,6 +552,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 :: localSize, itn !! for LSMR

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

Expand Down Expand Up @@ -746,6 +751,33 @@ 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

! TODO this should be an input
localSize = 0
!localSize = min(me%m, me%n)

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, 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

Check warning on line 772 in src/nlesolver_module.F90

View check run for this annotation

Codecov / codecov/patch

src/nlesolver_module.F90#L771-L772 

Added lines #L771 - L772 were not covered by tests
case(5)
call me%set_status(istat = -1005, &
string = 'LSMR Error: The iteration limit was reached. istop =', i=info)
exit

Check warning on line 776 in src/nlesolver_module.F90

View check run for this annotation

Codecov / codecov/patch

src/nlesolver_module.F90#L775-L776 

Added lines #L775 - L776 were not covered by tests
case default
info = 0
end select

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 89ffc2f

Please sign in to comment.