diff --git a/ford.md b/ford.md index 677518a..9d44d66 100644 --- a/ford.md +++ b/ford.md @@ -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 diff --git a/src/nlesolver_module.F90 b/src/nlesolver_module.F90 index 174483a..ff2a20c 100755 --- a/src/nlesolver_module.F90 +++ b/src/nlesolver_module.F90 @@ -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: diff --git a/test/sparse_test.f90 b/test/sparse_test.f90 index 70124fa..baf4dbc 100644 --- a/test/sparse_test.f90 +++ b/test/sparse_test.f90 @@ -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