Skip to content

Commit

Permalink
initial implementation of sparse Broyden option
Browse files Browse the repository at this point in the history
Fixes #13
jacobwilliams committed Jan 23, 2024
1 parent dd905e9 commit b519525
Showing 3 changed files with 84 additions and 46 deletions.
2 changes: 2 additions & 0 deletions ford.md
Original file line number Diff line number Diff line change
@@ -15,6 +15,8 @@ display: public
private
source: true
graph: true
externalize: true
external: fmin = https://jacobwilliams.github.io/fmin
extra_mods: iso_fortran_env: https://gcc.gnu.org/onlinedocs/gfortran/ISO_005fFORTRAN_005fENV.html
fmin_module: https://jacobwilliams.github.io/fmin
lusol_ez_module: https://jacobwilliams.github.io/lusol
124 changes: 82 additions & 42 deletions src/nlesolver_module.F90
Original file line number Diff line number Diff line change
@@ -494,12 +494,6 @@ subroutine initialize_nlesolver_variables(me,&
! LUSOL method
if (present(nout)) me%lusol_method = lusol_method

! now now, some options are not available for sparse mode
!...see if this is possible in sparse mode... TODO
if (me%use_broyden) then
call me%set_status(istat = -16, string = 'Error: broyden method not available for sparse solver')
return
end if
end if
end if

@@ -521,30 +515,34 @@ subroutine nlesolver_solver(me,x)
class(nlesolver_type),intent(inout) :: me
real(wp),dimension(:),intent(inout) :: x

real(wp),dimension(:) ,allocatable :: fvec !! function vector
real(wp),dimension(:,:),allocatable :: fjac !! jacobian matrix [dense]
real(wp),dimension(:), allocatable :: fjac_sparse !! jacobian matrix [sparse]
real(wp),dimension(:) ,allocatable :: rhs !! linear system right-hand side
real(wp),dimension(:) ,allocatable :: p !! search direction
real(wp),dimension(:) ,allocatable :: xold !! previous value of `x`
real(wp),dimension(:) ,allocatable :: prev_fvec !! previous function vector
real(wp),dimension(:,:),allocatable :: prev_fjac !! previous jacobian matrix
real(wp),dimension(:,:),allocatable :: delf !! used for Broyden (rank 2 for `matmul`)
real(wp),dimension(:,:),allocatable :: delx !! used for Broyden (rank 2 for `matmul`)
logical :: user_stop !! user stop button flag
integer :: info !! status flag from the [[linear_solver]]
integer :: iter !! iteration counter
real(wp) :: f !! magnitude of `fvec`
real(wp) :: fold !! previous value of `f`
integer :: n_uphill !! number of steps taken in the "uphill" direction
!! (where `f` is increasing)
real(wp) :: delxmag2 !! used for Broyden
logical :: recompute_jac !! if using Broyden, and we want to call the user
!! jacobian routine instead
integer :: broyden_counter !! number of times the broyden update has been used
integer :: alloc_stat !! allocation status flag
real(wp),dimension(:) ,allocatable :: fvec !! function vector
real(wp),dimension(:,:),allocatable :: fjac !! jacobian matrix [dense]
real(wp),dimension(:), allocatable :: fjac_sparse !! jacobian matrix [sparse]
real(wp),dimension(:), allocatable :: prev_fjac_sparse !! previous jacobian matrix [sparse]
real(wp),dimension(:) ,allocatable :: rhs !! linear system right-hand side
real(wp),dimension(:) ,allocatable :: p !! search direction
real(wp),dimension(:) ,allocatable :: xold !! previous value of `x`
real(wp),dimension(:) ,allocatable :: prev_fvec !! previous function vector
real(wp),dimension(:,:),allocatable :: prev_fjac !! previous jacobian matrix
real(wp),dimension(:,:),allocatable :: delf !! used for Broyden (rank 2 for `matmul`)
real(wp),dimension(:,:),allocatable :: delx !! used for Broyden (rank 2 for `matmul`)
logical :: user_stop !! user stop button flag
integer :: info !! status flag from the [[linear_solver]]
integer :: iter !! iteration counter
real(wp) :: f !! magnitude of `fvec`
real(wp) :: fold !! previous value of `f`
integer :: n_uphill !! number of steps taken in the "uphill" direction
!! (where `f` is increasing)
real(wp) :: delxmag2 !! used for Broyden
logical :: recompute_jac !! if using Broyden, and we want to call the user
!! jacobian routine instead
integer :: broyden_counter !! number of times the broyden update has been used
integer :: alloc_stat !! allocation status flag
type(lsqr_solver_ez) :: sparse_solver !! sparse LSQR solver class
type(lusol_settings) :: lusol_options
integer :: i !! counter
integer,dimension(:),allocatable :: idx, index_array !! for sparse indexing
character(len=10) :: i_str !! string version of `i` for row string

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

@@ -577,6 +575,7 @@ subroutine nlesolver_solver(me,x)
else
! sparse
if (alloc_stat==0) allocate(fjac_sparse (me%n_nonzeros) , stat=alloc_stat)
if (me%use_broyden .and. alloc_stat==0) allocate(prev_fjac_sparse(me%n_nonzeros) , stat=alloc_stat)
end if

if (alloc_stat==0) allocate(rhs (me%m) , stat=alloc_stat)
@@ -585,7 +584,10 @@ subroutine nlesolver_solver(me,x)
if (me%use_broyden) then
! only need these for broyden:
if (alloc_stat==0) allocate(prev_fvec(me%m) , stat=alloc_stat)
if (alloc_stat==0) allocate(prev_fjac(me%m,me%n) , stat=alloc_stat)
if (me%sparsity_mode/=1 .and. alloc_stat==0) then
allocate(prev_fjac(me%m,me%n) , stat=alloc_stat)
index_array = [(i, i=1,me%n_nonzeros)] ! an array to index the sparse jacobian elements
end if
if (alloc_stat==0) allocate(delf (me%m,1) , stat=alloc_stat)
if (alloc_stat==0) allocate(delx (me%n,1) , stat=alloc_stat)
end if
@@ -626,7 +628,13 @@ subroutine nlesolver_solver(me,x)
if (me%use_broyden .and. .not. recompute_jac) then
if (iter==1) then
! always compute Jacobian on the first iteration
call me%grad(x,fjac)
!call me%grad(x,fjac)
select case (me%sparsity_mode)
case (1)
call me%grad(x,fjac)
case (2:)
call me%grad_sparse(x,fjac_sparse)
end select
broyden_counter = 0
else
! and use Broyden update to estimate Jacobian
@@ -635,23 +643,55 @@ subroutine nlesolver_solver(me,x)
! note: fvec was computed the last iteration
delx(:,1) = x - xold
delf(:,1) = fvec - prev_fvec
delxmag2 = dot_product(delx(:,1),delx(:,1))

if (delxmag2 < eps) then
call me%set_status(istat = -8, &
string = 'Error: Divide by zero when computing Broyden update')
exit
end if
if (me%sparsity_mode==1) then ! dense

delxmag2 = dot_product(delx(:,1),delx(:,1))
if (delxmag2 < eps) then
call me%set_status(istat = -8, &
string = 'Error: Divide by zero when computing Broyden update')
exit
end if

! Jacobian estimate:
! This is the equation from wikipedia : https://en.wikipedia.org/wiki/Broyden%27s_method
fjac = prev_fjac + matmul((delf-matmul(prev_fjac,delx)), transpose(delx)) / delxmag2

! Jacobian estimate:
fjac = prev_fjac + &
matmul((delf-matmul(prev_fjac,delx)),&
transpose(delx))/delxmag2
! see also: eqn 4 in Schubert, which is different ...

else ! using a sparse option

! Just a row-by-row version of the fjac equation above.
! see L.K. Schubert, 1970
do i = 1, me%m

idx = pack(index_array, mask = me%irow==i) ! the nonzero elements in this row
if (size(idx)==0) cycle ! no non-zeros in this row

associate (dx => delx(me%icol(idx),:)) ! nonzero x vec for this row
delxmag2 = dot_product(dx(:,1),dx(:,1)) ! only those x's for this row
if (delxmag2 < eps) then
write(i_str,'(I10)') i
call me%set_status(istat = -8, &
string = 'Error: Divide by zero when computing sparse Broyden update for row '//&
trim(adjustl(i_str)))
exit
end if
fjac_sparse(idx) = prev_fjac_sparse(idx) + &
matmul((delf(i,1)-matmul(prev_fjac_sparse(idx),dx)), transpose(dx)) / delxmag2
end associate

end do

end if
broyden_counter = broyden_counter + 1
end if
prev_fjac = fjac
prev_fvec = fvec
if (me%sparsity_mode==1) then
prev_fjac = fjac
else
prev_fjac_sparse = fjac_sparse
end if
else
! compute the jacobian:
select case (me%sparsity_mode)
@@ -702,7 +742,7 @@ subroutine nlesolver_solver(me,x)
! 1 => TRP: Threshold Rook Pivoting.
! 2 => TCP: Threshold Complete Pivoting.
call solve(me%n,me%m,me%n_nonzeros,me%irow,me%icol,fjac_sparse,rhs,p,info,&
settings=lusol_options)
settings=lusol_options)
end select

! check for errors:
4 changes: 0 additions & 4 deletions test/sparse_test.f90
Original file line number Diff line number Diff line change
@@ -48,7 +48,6 @@ program sparse_test
f_evals = 0
description = 'Constant alpha'
case(2)
cycle ! not yet supported
step_mode = 1
use_broyden = .true.
f_evals = 0
@@ -59,7 +58,6 @@ program sparse_test
f_evals = 0
description = 'Backtracking line search'
case(4)
cycle ! not yet supported
step_mode = 2
use_broyden = .true.
f_evals = 0
@@ -70,7 +68,6 @@ program sparse_test
f_evals = 0
description = 'Exact line search'
case(6)
cycle ! not yet supported
step_mode = 3
use_broyden = .true.
f_evals = 0
@@ -81,7 +78,6 @@ program sparse_test
f_evals = 0
description = 'Fixed point search'
case(8)
cycle ! not yet supported
step_mode = 4
use_broyden = .true.
f_evals = 0

0 comments on commit b519525

Please sign in to comment.