### SLSQP Working Array Address Setup and SLSQPB Call

Source: https://jacobwilliams.github.io/slsqp/proc/slsqp.html

Sets up initial addresses in the working array for various subarrays, assigns the n1 field of sdat, and calls the SLSQPB subroutine with all required arguments.

```fortran
    ! prepare data for calling sqpbdy  -  initial addresses in w

    im = 1
    il = im + max(1,m)
    il = im + la
    ix = il + n1*n/2 + 1
    ir = ix + n
    is = ir + n + n + max(1,m)
    is = ir + n + n + la
    iu = is + n1
    iv = iu + n1
    iw = iv + n1

    sdat%n1 = n1

    call slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,&
                w(ir),w(il),w(ix),w(im),w(is),w(iu),w(iv),w(iw),&
                sdat%t,sdat%f0,sdat%h1,sdat%h2,sdat%h3,sdat%h4,&
                sdat%n1,sdat%n2,sdat%n3,sdat%t0,sdat%gs,sdat%tol,sdat%line,&
                sdat%alpha,sdat%iexact,sdat%incons,sdat%ireset,sdat%itermx,&
                ldat,alphamin,alphamax,tolf,toldf,toldx,&
                max_iter_ls,nnls_mode,infBnd)

    end subroutine slsqp
```

--------------------------------

### SLSQP solver main routine (Fortran)

Source: https://jacobwilliams.github.io/slsqp/sourcefile/slsqp_module.f90.html

The main solver routine for the SLSQP class. It declares local variables, initializes arrays, checks setup conditions, and enters the main solver loop. The routine uses reverse communication with the underlying SLSQP algorithm.

```fortran
    class(slsqp_solver),intent(inout)   :: me
    real(wp),dimension(:),intent(inout) :: x          !! **in:**  initial optimization variables,
                                                      !! **out:** solution.
    integer,intent(out)                 :: istat      !! status code (see `mode` in [[slsqp]]).
    integer,intent(out),optional        :: iterations !! number of iterations
    character(len=:),intent(out),allocatable,optional :: status_message
                                                      !! string status message
                                                      !! corresponding to `istat`

    ! local variables:
    real(wp),dimension(:),allocatable   :: c        !! constraint vector -- `dimension(max(1,me%m))`
    real(wp),dimension(:,:),allocatable :: a        !! a matrix for [[slsqp]] -- `dimension(max(1,me%m),me%n+1)`
    real(wp),dimension(:),allocatable   :: g        !! g matrix for [[slsqp]] -- `dimension(me%n+1)`
    real(wp),dimension(:),allocatable   :: cvec     !! constraint vector -- `dimension(me%m)`
    real(wp),dimension(:),allocatable   :: dfdx     !! objective function partials -- `dimension(me%n)`
    real(wp),dimension(:,:),allocatable :: dcdx     !! constraint partials -- `dimension(me%m,me%n)`
    real(wp),dimension(:),allocatable   :: delta    !! perturbation step size to approximate gradient -- `dimension(me%n)`
    real(wp),dimension(:),allocatable   :: cvecr    !! right function value to approximate constraints vector's gradient -- `dimension(me%m)`
    real(wp),dimension(:),allocatable   :: cvecl    !! left function value to approximate constraints vector's gradient -- `dimension(me%m)`
    real(wp)                            :: f        !! objective function
    integer                             :: i        !! iteration counter
    integer                             :: mode     !! reverse communication flag for [[slsqp]]
    integer                             :: la       !! input to [[slsqp]]
    integer                             :: iter     !! in/out for [[slsqp]]
    real(wp)                            :: acc      !! in/out for [[slsqp]]
    integer                             :: ig       !! loop index to approximate gradient
    real(wp)                            :: fr       !! right function value to approximate objective function's gradient
    real(wp)                            :: fl       !! left function value to approximate objective function's gradient
    real(wp)                            :: fact     !! denominator factor for finite difference approximation

    !initialize:
    allocate(c(max(1,me%m))       )
    allocate(a(max(1,me%m),me%n+1))
    allocate(g(me%n+1)            )
    allocate(cvec(me%m)           )
    allocate(dfdx(me%n)           )
    allocate(dcdx(me%m,me%n)      )
    allocate(delta(me%n)          )
    allocate(cvecr(me%m)          )
    allocate(cvecl(me%m)          )
    i    = 0
    iter = me%max_iter
    la   = max(1,me%m)
    mode = 0
    a    = zero
    g    = zero
    c    = zero
    if (present(iterations)) iterations = 0
    call me%linmin%destroy()
    call me%slsqpb%destroy()

    !check setup:
    if (size(x)/=me%n) then
        istat = -100
        call me%report_message(mode_to_status_message(istat))
        if (present(status_message)) status_message = mode_to_status_message(istat)
        return
    end if

    !linesearch:
    select case(me%linesearch_mode)
    case(1)     !inexact (armijo-type linesearch)
        acc = abs(me%acc)
    case(2)     !exact
        acc = -abs(me%acc)
    case default
        istat = -101
        call me%report_message(mode_to_status_message(istat))
        if (present(status_message)) status_message = mode_to_status_message(istat)
        return
    end select

    !make sure the functions have been associated:
    if (.not. associated(me%f)) then
        istat = -102
        call me%report_message(mode_to_status_message(istat))
        if (present(status_message)) status_message = mode_to_status_message(istat)
        return
    end if
    if ((me%gradient_mode==0).and.(.not. associated(me%g))) then
        istat = -103
        call me%report_message(mode_to_status_message(istat))
        if (present(status_message)) status_message = mode_to_status_message(istat)
        return
    end if
    if (me%gradient_mode<0 .or. me%gradient_mode>3) then
        istat = -104
        call me%report_message(mode_to_status_message(istat))
        if (present(status_message)) status_message = mode_to_status_message(istat)
        return
    end if
    if (me%gradient_mode/=0 .and. me%gradient_delta<=epmach) then
        istat = -105
        call me%report_message(mode_to_status_message(istat))
        if (present(status_message)) status_message = mode_to_status_message(istat)
        return
    end if

    !main solver loop:
    do
```

--------------------------------

### initialize

Source: https://jacobwilliams.github.io/slsqp/module/slsqp_module.html

Initializes the SLSQP solver with problem dimensions and settings.

```APIDOC
## initialize

### Description
Initializes the SLSQP solver with problem dimensions and settings.

### Method
initialize

### Parameters
- **n** (integer) - Required - Number of optimization variables (n>0)
- **m** (integer) - Required - Number of constraints (m≥0)
- **meq** (integer) - Required - Number of equality constraints (m≥meq≥0)
- **max_iter** (integer) - Required - Maximum number of iterations
- **acc** (real) - Required - Accuracy tolerance
- **tolf** (real) - Optional - Accuracy tolerance over f: if |f|<tolf then stop
- **toldf** (real) - Optional - Accuracy tolerance over df: if |fn+1−fn|<toldf then stop
- **toldx** (real) - Optional - Accuracy tolerance over dx: if |xn+1−xn|<toldx then stop
- **gradient_mode** (integer) - Optional - How the gradients are computed
- **gradient_delta** (real) - Optional - Perturbation step size to approximate gradients by finite differences
- **alphamin** (real) - Optional - Min α for line search
- **alphamax** (real) - Optional - Max α for line search
- **iprint** (integer) - Optional - Unit number of status printing (0 for no printing)
- **xl** (real array) - Optional - Lower bound on x
- **xu** (real array) - Optional - Upper bound on x
- **linesearch_mode** (integer) - Optional - Linesearch mode
- **nnls_mode** (integer) - Optional - Which NNLS method to use
- **max_iter_ls** (integer) - Optional - Max iterations in the least squares problem
- **infinite_bound** (real) - Optional - "Infinity" for the upper and lower bounds

### Response
No return value.
```

--------------------------------

### Build SLSQP with fpm

Source: https://jacobwilliams.github.io/slsqp/index.html

Build the library in release mode using the Fortran Package Manager.

```bash
fpm build --release
```

--------------------------------

### initialize_slsqp subroutine implementation

Source: https://jacobwilliams.github.io/slsqp/proc/initialize_slsqp.html

Fortran implementation of the initialize_slsqp subroutine. It checks bounds vector sizes, meq, m, n, and bound consistency, then sets optional parameters (linesearch_mode, alphamin, alphamax, tolf, toldf, toldx, max_iter_ls, nnls_mode, gradient_mode, gradient_delta, infinite_bound), allocates work arrays, and initializes solver fields. Reports errors via report_message and returns early on invalid inputs.

```fortran
    if (size(xl)/=size(xu) .or. size(xl)/=n) then
        call me%report_message('error: invalid upper or lower bound vector size')
        call me%report_message('  size(xl) =',ival=size(xl))
        call me%report_message('  size(xu) =',ival=size(xu))
        call me%report_message('  n        =',ival=n)
    else if (meq<0 .or. meq>m) then
        call me%report_message('error: invalid meq value:', ival=meq)
    else if (m<0) then
        call me%report_message('error: invalid m value:', ival=m)
    else if (n<1) then
        call me%report_message('error: invalid n value:', ival=n)
    else if (any(xl>xu .and. .not. ieee_is_nan(xl) .and. .not. ieee_is_nan(xu))) then
        call me%report_message('error: lower bounds must be <= upper bounds.')
        do i=1,n
            if (xl(i)>xu(i) .and. .not. ieee_is_nan(xl(i)) .and. .not. ieee_is_nan(xu(i))) then
                call me%report_message('  xl(i)>xu(i) for variable',ival=i)
            end if
        end do
    else

        if (present(linesearch_mode)) then
            !two linesearch modes:
            select case (linesearch_mode)
            case(1)     !inexact
                me%linesearch_mode = linesearch_mode
            case(2)     !exact
                me%linesearch_mode = linesearch_mode
            case default
                call me%report_message('error: invalid linesearch_mode (must be 1 or 2): ',&
                                        ival=linesearch_mode)
                call me%destroy()
                return
            end select
        end if

        !optional linesearch bounds:
        if (present(alphamin)) me%alphamin = alphamin
        if (present(alphamax)) me%alphamax = alphamax

        !verify valid values for alphamin and alphamax: 0<alphamin<alphamax<=1
        if (me%alphamin<=zero .or. me%alphamax<=zero .or. &
            me%alphamax<=me%alphamin .or. &
            me%alphamin>=one .or. me%alphamax>one) then

            call me%report_message('error: invalid values for alphamin or alphamax.')
            call me%report_message('  alphamin =',rval=me%alphamin)
            call me%report_message('  alphamax =',rval=me%alphamax)
            call me%destroy()
            return

        end if

        if (present(tolf))  me%tolf  = tolf
        if (present(toldf)) me%toldf = toldf
        if (present(toldx)) me%toldx = toldx

        if (present(max_iter_ls)) me%max_iter_ls = max_iter_ls
        if (present(nnls_mode)) then
            select case (nnls_mode)
            case(1:2)
                me%nnls_mode = nnls_mode
            case default
                call me%report_message('error: invalid value for nnls_mode. defaulting to 1.')
                me%nnls_mode = 1
            end select
        end if

        status_ok = .true.
        me%n = n
        me%m = m
        me%meq = meq
        me%max_iter = max_iter
        me%acc = acc
        me%f => f
        me%g => g
        if (present(report)) me%report => report

        allocate(me%xl(n)); me%xl = xl
        allocate(me%xu(n)); me%xu = xu

        !work arrays:
        n1 = n+1
        mineq = m - meq + 2*n1

        me%l_w = n1*(n1+1) + meq*(n1+1) + mineq*(n1+1) + &   !for lsq
                 (n1-meq+1)*(mineq+2) + 2*mineq        + &   !for lsi
                 (n1+mineq)*(n1-meq) + 2*meq + n1      + &   !for lsei
                  n1*n/2 + 2*m + 3*n +3*n1 + 1               !for slsqpb

        allocate(me%w(me%l_w))
        me%w = zero

        if (present(gradient_mode)) then
            me%gradient_mode = gradient_mode
            if (present(gradient_delta)) then
                me%gradient_delta = gradient_delta
            end if
        end if

        if (present(infinite_bound)) then
            me%infinite_bound = abs(infinite_bound)
        else
            me%infinite_bound = huge(one)
        end if

    end if

    end subroutine initialize_slsqp
```

--------------------------------

### SLSQP lsq subroutine declaration and variable setup

Source: https://jacobwilliams.github.io/slsqp/sourcefile/slsqp_core.f90.html

Declares the lsq subroutine interface, including input arrays l, g, a, b, xl, xu, and working array w, along with output arrays x and y. Defines the mode flag meanings and initializes internal variables for the optimization.

```fortran
    subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode,infBnd)

    implicit none

    integer,intent(in)        :: m
    integer,intent(in)        :: n
    integer,intent(in)        :: meq
    integer,intent(in)        :: nl
    integer,intent(in)        :: la
    real(wp),dimension(n)     :: x  !! stores the n-dimensional solution vector
    real(wp),dimension(m+n+n) :: y  !! stores the vector of lagrange multipliers of dimension
                                    !! m+n+n (constraints+lower+upper bounds)
    integer :: mode                 !! is a success-failure flag with the following meanings:
                                    !!
                                    !! * **1:** successful computation,
                                    !! * **2:** error return because of wrong dimensions (`n<1`),
                                    !! * **3:** iteration count exceeded by [[nnls]],
                                    !! * **4:** inequality constraints incompatible,
                                    !! * **5:** matrix `e` is not of full rank,
                                    !! * **6:** matrix `c` is not of full rank,
                                    !! * **7:** rank defect in [[hfti]]
    integer,intent(in) :: max_iter_ls !! maximum number of iterations in the [[nnls]] problem
    integer,intent(in) :: nnls_mode !! which NNLS method to use
    real(wp),intent(in) :: infbnd !! "infinity" for the upper and lower bounds.

    real(wp),dimension(nl)   :: l
    real(wp),dimension(n)    :: g
    real(wp),dimension(la,n) :: a
    real(wp),dimension(la)   :: b
    real(wp),dimension(*)    :: w
    real(wp),dimension(n)    :: xl
    real(wp),dimension(n)    :: xu
    real(wp) :: diag , xnorm
    integer :: i , ic , id , ie , iff , ig , ih , il , ip , &
               iw , i1 , i2 , i3 , i4 , mineq , &
               m1 , n1 , n2 , n3, num_unbounded, j

      n1 = n + 1
      mineq = m - meq
      m1 = mineq + n + n

      !  determine whether to solve problem
      !  with inconsistent linerarization (n2=1)
      !  or not (n2=0)

      n2 = n1*n/2 + 1
      if ( n2==nl ) then
         n2 = 0
      else
         n2 = 1
      end if
      n3 = n - n2

      !  recover matrix e and vector f from l and g

      i2 = 1
      i3 = 1
      i4 = 1
      ie = 1
      iff = n*n + 1
      do i = 1 , n3
         i1 = n1 - i
         diag = sqrt(l(i2))
         w(i3) = zero
         call dcopy(i1,w(i3),0,w(i3),1)
         call dcopy(i1-n2,l(i2),1,w(i3),n)
         call dscal(i1-n2,diag,w(i3),n)
         w(i3) = diag
         w(iff-1+i) = (g(i)-ddot(i-1,w(i4),1,w(iff),1))/diag
         i2 = i2 + i1 - n2
         i3 = i3 + n1
         i4 = i4 + n
      end do
      if ( n2==1 ) then
         w(i3) = l(nl)
         w(i4) = zero
         call dcopy(n3,w(i4),0,w(i4),1)
         w(iff-1+n) = zero
      end if
      call dscal(n,-one,w(iff),1)

      ic = iff + n
      id = ic + meq*n

      if ( meq>0 ) then

         !  recover matrix c from upper part of a

         do i = 1 , meq
            call dcopy(n,a(i,1),la,w(ic-1+i),meq)
         end do

         !  recover vector d from upper part of b

         call dcopy(meq,b(1),1,w(id),1)
         call dscal(meq,-one,w(id),1)

      end if

      ig = id + meq

      if ( mineq>0 ) then
         ! recover matrix g from lower part of a
         ! The matrix G(mineq+2*n,m1) is stored at w(ig)
         ! Not all rows will be filled if some of the upper/lower
         ! bounds are unbounded.
         do i = 1 , mineq
            call dcopy(n,a(meq+i,1),la,w(ig-1+i),m1)
         end do
      end if

      ih = ig + m1*n
      iw = ih + mineq + 2*n
```

--------------------------------

### Define the slsqp_solver derived type

Source: https://jacobwilliams.github.io/slsqp/type/slsqp_solver.html

This is the full definition of the `slsqp_solver` type, the main class for interfacing with the SLSQP solver. It includes all components for problem setup (number of variables, constraints, tolerances), gradient computation modes, line search settings, and work arrays. It also declares the public procedures `initialize`, `destroy`, `optimize`, and `abort`.

```fortran
    type,public :: slsqp_solver

        !! The main class used to interface with the SLSQP solver.

        private

        integer  :: n           = 0        !! number of optimization variables (\( n > 0 \))
        integer  :: m           = 0        !! number of constraints (\( m \ge 0 \))
        integer  :: meq         = 0        !! number of equality constraints (\( m \ge m_{eq} \ge 0 \))
        integer  :: max_iter    = 0        !! maximum number of iterations
        real(wp) :: acc         = zero     !! accuracy tolerance
        real(wp) :: tolf        = -one     !! accuracy tolerance over f:  if \( |f| < tolf \) then stop
        real(wp) :: toldf       = -one     !! accuracy tolerance over df: if \( |f_{n+1} - f_n| < toldf \) then stop.
                                           !! It's different from `acc` in the case of positive derivative
        real(wp) :: toldx       = -one     !! accuracy tolerance over dx: if \( |x_{n+1} - x_n| < toldx \) then stop

        integer  :: gradient_mode = 0      !! how the gradients are computed:
                                           !!
                                           !! * 0 - use the user-supplied `g` subroutine. [default]
                                           !! * 1 - approximate by basic backward differences
                                           !! * 2 - approximate by basic forward differences
                                           !! * 3 - approximate by basic central differences
        real(wp) :: gradient_delta  = 1.0e-8_wp !! perturbation step size to approximate gradients
                                                !! by finite differences (`gradient_mode` 1-3).

        !these two were not in the original code:
        real(wp) :: alphamin = 0.1_wp   !! min \( \alpha \) for line search \( 0 < \alpha_{min} < \alpha_{max} \le 1 \)
        real(wp) :: alphamax = 1.0_wp   !! max \( \alpha \) for line search \( 0 < \alpha_{min} < \alpha_{max} \le 1 \)

        integer :: iprint = output_unit !! unit number of status printing (0 for no printing)

        real(wp),dimension(:),allocatable :: xl  !! lower bound on x
        real(wp),dimension(:),allocatable :: xu  !! upper bound on x

        integer :: l_w = 0 !! size of `w`
        real(wp),dimension(:),allocatable :: w !! real work array

        procedure(func),pointer     :: f          => null() !! problem function subroutine
        procedure(grad),pointer     :: g          => null() !! gradient subroutine
        procedure(iterfunc),pointer :: report     => null() !! for reporting an iteration
        procedure(msgfunc),pointer  :: report_msg => null() !! for reporting a string message from slsqp to the user

        integer :: linesearch_mode = 1  !! linesearch mode:
                                        !!
                                        !! * `1` = inexact (Armijo) linesearch,
                                        !! * `2` = exact linesearch.
        type(linmin_data) :: linmin !! data formerly within [[linmin]].
                                    !! Only used when `linesearch_mode=2`
        type(slsqpb_data) :: slsqpb  !! data formerly within [[slsqpb]].

        ! note: the following two maybe should be combined into a separate type
        ! along with the two methods...
        integer :: nnls_mode = 1 !! Which NNLS method to use:
                                 !!
                                 !! 1. Use the original [[nnls]]
                                 !! 2. Use the newer [[bvls]]
        integer :: max_iter_ls = 0  !! max iterations in the least squares problem.
                                    !! if `<=0`, defaults to `3*n`.
                                    !! (use by either [[nnls]] or [[bvls]])

        logical :: user_triggered_stop = .false.    !! if the `abort` method has been called
                                                    !! to stop the iterations

        real(wp) :: infinite_bound = huge(one)  !! "infinity" for the upper and lower bounds.
                                                !! if `xl<=-infinite_bound` or `xu>=infinite_bound`
                                                !! then these bounds are considered nonexistant.

    contains

        private

        procedure,public :: initialize => initialize_slsqp
        procedure,public :: destroy    => destroy_slsqp
        procedure,public :: optimize   => slsqp_wrapper
        procedure,public :: abort      => stop_iterations

        procedure :: report_message  !! for reporting messages to the user

    end type slsqp_solver
```

--------------------------------

### initialize_slsqp

Source: https://jacobwilliams.github.io/slsqp/proc/initialize_slsqp.html

Initializes the SLSQP solver with problem dimensions, objective and constraint functions, bounds, and optional configuration parameters.

```APIDOC
## initialize_slsqp

### Description
Initializes the SLSQP solver instance with problem dimensions, objective and constraint functions, bounds, and optional configuration parameters.

### Method
Subroutine call

### Signature
```fortran
subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,&
                            linesearch_mode,iprint,report,alphamin,alphamax,&
                            gradient_mode,gradient_delta,tolf,toldf,toldx,&
                            max_iter_ls,nnls_mode,infinite_bound,report_msg)
```

### Parameters
#### Required Parameters
- **me** (class(slsqp_solver)) - The solver instance to initialize.
- **n** (integer) - Number of variables, \( n \ge 1 \).
- **m** (integer) - Total number of constraints, \( m \ge 0 \).
- **meq** (integer) - Number of equality constraints, \( m_{eq} \ge 0 \).
- **max_iter** (integer) - Maximum number of iterations.
- **f** (procedure(func)) - Problem function.
- **g** (procedure(grad)) - Function to compute gradients (must be associated if `gradient_mode=0`).
- **xl** (real(wp), dimension(n)) - Lower bounds on `x`. `xl(i)=NaN` (or `xl(i)<=-infinite_bound`) indicates to ignore `i`th bound.
- **xu** (real(wp), dimension(n)) - Upper bounds on `x`. `xu(i)=NaN` (or `xu(i)>=infinite_bound`) indicates to ignore `i`th bound.
- **acc** (real(wp)) - Accuracy.
- **status_ok** (logical, intent(out)) - Will be false if there were errors.

#### Optional Parameters
- **linesearch_mode** (integer) - 1 = inexact (default), 2 = exact.
- **iprint** (integer) - Unit number of status messages (default=`output_unit`).
- **report** (procedure(iterfunc)) - User-defined procedure called once per iteration.
- **alphamin** (real(wp)) - Minimum alpha for linesearch [default 0.1].
- **alphamax** (real(wp)) - Maximum alpha for linesearch [default 1.0].
- **gradient_mode** (integer) - How gradients are computed: 0 = user-supplied `g` (default), 1 = backward differences, 2 = forward differences, 3 = central differences. Modes 1-3 do not respect variable bounds.
- **gradient_delta** (real(wp)) - Perturbation step size (>epsilon) for finite difference gradients (modes 1-3). Absolute step that does not respect bounds.
- **tolf** (real(wp)) - Stopping criterion: stop if \( |f| < tolf \).
- **toldf** (real(wp)) - Stopping criterion: stop if \( |f_{n+1} - f_n| < toldf \).
- **toldx** (real(wp)) - Stopping criterion: stop if \( ||x_{n+1} - x_n|| < toldx \).
- **max_iter_ls** (integer) - Maximum number of iterations in the [[nnls]] problem.
- **nnls_mode** (integer) - Which NNLS method to use: 1 = original [[nnls]], 2 = newer [[bvls]].
- **infinite_bound** (real(wp)) - "Infinity" for upper and lower bounds. If `xl<=-infinite_bound` or `xu>=infinite_bound`, bounds are considered nonexistent. Default is `huge()`.
- **report_msg** (procedure(msgfunc)) - User-defined procedure called with any warning or error messages.

### Behavior
- Sets `status_ok` to `.false.` initially.
- Calls `me%destroy()` to reset the solver.
- If `iprint` is present, sets `me%iprint`.
- If `report_msg` is present, associates `me%report_msg`.

### Notes
- The subroutine does not return a value; it initializes the solver object.
- The `g` procedure must be associated if `gradient_mode=0`.
- Finite difference gradient modes (1-3) do not respect variable bounds.
- The `infinite_bound` parameter defines what is considered an infinite bound; if not provided, `huge()` is used.
```

--------------------------------

### initialize

Source: https://jacobwilliams.github.io/slsqp/sourcefile/slsqp_module.f90.html

Initializes the SLSQP solver with problem dimensions, objective function, gradient function, bounds, and optional solver settings. This method must be called before solving.

```APIDOC
## initialize

### Description
Initializes the SLSQP solver with the problem definition and solver options. This method must be called before solving.

### Method
initialize

### Parameters
- **me** (class(slsqp_solver)) - Required - The solver instance to initialize.
- **n** (integer) - Required - The number of variables, n >= 1.
- **m** (integer) - Required - Total number of constraints, m >= 0.
- **meq** (integer) - Required - Number of equality constraints, meq >= 0.
- **max_iter** (integer) - Required - Maximum number of iterations.
- **f** (procedure(func)) - Required - Problem function.
- **g** (procedure(grad)) - Required - Function to compute gradients (must be associated if gradient_mode=0).
- **xl** (real(wp), dimension(n)) - Required - Lower bounds on x. xl(i)=NaN (or xl(i)<=-infinite_bound) indicates to ignore i-th bound.
- **xu** (real(wp), dimension(n)) - Required - Upper bounds on x. xu(i)=NaN (or xu(i)>=infinite_bound) indicates to ignore i-th bound.
- **acc** (real(wp)) - Required - Accuracy.
- **status_ok** (logical) - Required - Will be false if there were errors.
- **linesearch_mode** (integer, optional) - Optional - 1 = inexact (default), 2 = exact.
- **iprint** (integer, optional) - Optional - Unit number of status messages (default=output_unit).
- **report** (procedure(iterfunc), optional) - Optional - User-defined procedure that will be called once per iteration.
- **alphamin** (real(wp), optional) - Optional - Minimum alpha for linesearch [default 0.1].
- **alphamax** (real(wp), optional) - Optional - Maximum alpha for linesearch [default 1.0].
- **gradient_mode** (integer, optional) - Optional - How the gradients are to be computed: 0 - use the user-supplied g subroutine. [default]; 1 - approximate by basic backward differences; 2 - approximate by basic forward differences; 3 - approximate by basic central differences. Note that modes 1-3 do not respect the variable bounds.
- **gradient_delta** (real(wp), optional) - Optional - Perturbation step size (>epsilon) to compute the approximated gradient by finite differences (gradient_mode 1-3). Note that this is an absolute step that does not respect the xl or xu variable bounds.
- **tolf** (real(wp), optional) - Optional - Stopping criterion if |f| < tolf then stop.
- **toldf** (real(wp), optional) - Optional - Stopping criterion if |f_{n+1} - f_n| < toldf then stop.
- **toldx** (real(wp), optional) - Optional - Stopping criterion if ||x_{n+1} - x_n|| < toldx then stop.
- **max_iter_ls** (integer, optional) - Optional - Maximum number of iterations in the nnls problem.
- **nnls_mode** (integer, optional) - Optional - Which NNLS method to use: 1. Use the original nnls; 2. Use the newer bvls.
- **infinite_bound** (real(wp), optional) - Optional - "infinity" for the upper and lower bounds. If xl<=-infinite_bound or xu>=infinite_bound then these bounds are considered nonexistant. If not present then huge() is used for this.
- **report_msg** (procedure(msgfunc), optional) - Optional - User-defined procedure that will be called with any warning or error messages.

### Return Value
- **status_ok** (logical) - Indicates if initialization was successful (false if errors).
```

--------------------------------

### BVLS initialize subroutine

Source: https://jacobwilliams.github.io/slsqp/sourcefile/bvls_module.f90.html

Initializes the BVLS algorithm. Checks array sizes and sets error codes (IERR=1 for invalid dimensions, IERR=2 for inconsistent sizes, IERR=3 for invalid range). Sets iteration limit to 3*n if max_iter is not positive. Initializes index array and sets up the initial X values based on bounds, adjusting B for nonzero starting values.

```fortran
subroutine initialize()

   m=size(a,1)
   n=size(a,2)
   if (m  <=  0 .or. n  <=  0) then
      ierr = 1
      return
   end if

   ! Check array sizes for consistency and with M and N.
   if (size(x) < n) then
      ierr=2
      return
   end if
   if (size(b) < m) then
      ierr=2
      return
   end if
   if (size(bnd,1) /= 2) then
      ierr=2
      return
   end if
   if (size(bnd,2) < n) then
      ierr=2
      return
   end if
   if (size(w) < n) then
      ierr=2
      return
   end if
   if (size(index) < n) then
      ierr=2
      return
   end if

   ierr = 0
   if (max_iter<=0) then
      itmax = 3*n
   else
      itmax = max_iter
   end if
   iter=0
   ! Initialize the array index().
   do i=1,n
      index(i)=i
   end do

   iz2=n
   iz1=1
   nsetp=0
   npp1=1

   ! Begin:  Loop on IZ to initialize  X().
   iz=iz1
   do
      if (iz  >  iz2 ) exit
      j=index(iz)
      if ( bnd(1,j)   <= -huge(one)) then
         if (bnd(2,j) >=  huge(one)) then
            x(j) = zero
         else
            x(j) = min(zero,bnd(2,j))
         end if
     else  if ( bnd(2,j) >= huge(one)) then
        x(j) = max(zero,bnd(1,j))
     else
        range = bnd(2,j) - bnd(1,j)
        if ( range <= zero ) then
           ! Here X(J) is constrained to a single value.
           index(iz)=index(iz2)
           index(iz2)=j
           iz=iz-1
           iz2=iz2-1
           x(j)=bnd(1,j)
           w(j)=zero
         else  if ( range  >  zero) then
            !! The following statement sets X(J) to 0 if the constraint interval
            !! includes 0, and otherwise sets X(J) to the endpoint of the
            !! constraint interval that is closest to 0.
            x(j) = max(bnd(1,j), min(bnd(2,j),zero))
         else
            ierr = 3
            return
         end if
   end if
      if ( abs(x(j)) > zero ) then
         ! Change B() to reflect a nonzero starting value for X(J).
         b(1:m)=b(1:m)-a(1:m,j)*x(j)
      end if
      iz=iz+1
   end do

end subroutine initialize
```

--------------------------------

### optimize

Source: https://jacobwilliams.github.io/slsqp/module/slsqp_module.html

Runs the SLSQP optimization algorithm to find the optimal solution.

```APIDOC
## optimize

### Description
Runs the SLSQP optimization algorithm to find the optimal solution.

### Method
optimize

### Parameters
- **f** (procedure) - Required - Problem function subroutine
- **g** (procedure) - Required - Gradient subroutine
- **report** (procedure) - Optional - For reporting an iteration
- **report_msg** (procedure) - Optional - For reporting a string message from slsqp to the user

### Response
No return value.
```

--------------------------------

### slsqp_wrapper

Source: https://jacobwilliams.github.io/slsqp/src/slsqp_module.f90

Runs the SLSQP optimization algorithm using the provided objective function, gradient, and iteration callbacks. This is the main optimization routine.

```APIDOC
## slsqp_wrapper

### Description
Runs the SLSQP optimization algorithm using the provided objective function, gradient, and iteration callbacks. This is the main optimization routine.

### Method
optimize

### Signature
call me%optimize(func, grad, iterfunc)

### Parameters
- **func** (procedure) - Required - User-defined subroutine to compute objective and constraints.
- **grad** (procedure) - Required - User-defined subroutine to compute gradients.
- **iterfunc** (procedure) - Required - User-defined subroutine called at each iteration.

### Callback Interfaces
#### func
subroutine func(me, x, f, c)
- **me** (class(slsqp_solver)) - The solver instance.
- **x** (real array) - Optimization variable vector.
- **f** (real) - Objective function value.
- **c** (real array) - Constraint vector (equality constraints first).

#### grad
subroutine grad(me, x, g, a)
- **me** (class(slsqp_solver)) - The solver instance.
- **x** (real array) - Optimization variable vector.
- **g** (real array) - Objective function partials w.r.t x.
- **a** (real matrix) - Gradient matrix of constraints w.r.t x.

#### iterfunc
subroutine iterfunc(me, iter, x, f, c)
- **me** (class(slsqp_solver)) - The solver instance.
- **iter** (integer) - Iteration number.
- **x** (real array) - Optimization variable vector.
- **f** (real) - Objective function value.
- **c** (real array) - Constraint vector (equality constraints first).
```

--------------------------------

### Build SLSQP with single precision flag

Source: https://jacobwilliams.github.io/slsqp/index.html

Build a single precision version of the library by specifying the REAL32 preprocessor flag.

```bash
fpm build --profile release --flag "-DREAL32"
```

--------------------------------

### initialize_slsqp subroutine definition

Source: https://jacobwilliams.github.io/slsqp/proc/initialize_slsqp.html

Defines the initialize_slsqp subroutine, which initializes the SLSQP solver with problem dimensions, bounds, and optional settings. It sets status_ok to false and destroys the solver instance before applying optional parameters.

```fortran
    subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,&
                                linesearch_mode,iprint,report,alphamin,alphamax,&
                                gradient_mode,gradient_delta,tolf,toldf,toldx,&
                                max_iter_ls,nnls_mode,infinite_bound,report_msg)

    implicit none

    class(slsqp_solver),intent(inout) :: me
    integer,intent(in)                :: n               !! the number of variables, \( n \ge 1 \)
    integer,intent(in)                :: m               !! total number of constraints, \( m \ge 0 \)
    integer,intent(in)                :: meq             !! number of equality constraints, \( m_{eq} \ge 0 \)
    integer,intent(in)                :: max_iter        !! maximum number of iterations
    procedure(func)                   :: f               !! problem function
    procedure(grad)                   :: g               !! function to compute gradients (must be
                                                         !! associated if `gradient_mode=0`)
    real(wp),dimension(n),intent(in)  :: xl              !! lower bounds on `x`.
                                                         !! `xl(i)=NaN` (or `xl(i)<=-infinite_bound`) indicates to ignore `i`th bound
    real(wp),dimension(n),intent(in)  :: xu              !! upper bounds on `x`.
                                                         !! `xu(i)=NaN` (or `xu(i)>=infinite_bound`) indicates to ignore `i`th bound
    real(wp),intent(in)               :: acc             !! accuracy
    logical,intent(out)               :: status_ok       !! will be false if there were errors
    integer,intent(in),optional       :: linesearch_mode !! 1 = inexact (default), 2 = exact
    integer,intent(in),optional       :: iprint          !! unit number of status messages (default=`output_unit`)
    procedure(iterfunc),optional      :: report          !! user-defined procedure that will be called once per iteration
    real(wp),intent(in),optional      :: alphamin        !! minimum alpha for linesearch [default 0.1]
    real(wp),intent(in),optional      :: alphamax        !! maximum alpha for linesearch [default 1.0]
    integer,intent(in),optional       :: gradient_mode   !! how the gradients are to be computed:
                                                         !!
                                                         !! * 0 - use the user-supplied `g` subroutine. [default]
                                                         !! * 1 - approximate by basic backward differences
                                                         !! * 2 - approximate by basic forward differences
                                                         !! * 3 - approximate by basic central differences
                                                         !!
                                                         !! Note that modes 1-3 do not respect the variable bounds.
    real(wp),intent(in),optional      :: gradient_delta  !! perturbation step size (>epsilon) to compute the approximated
                                                         !! gradient by finite differences (`gradient_mode` 1-3).
                                                         !! note that this is an absolute step that does not respect
                                                         !! the `xl` or `xu` variable bounds.
    real(wp),intent(in),optional      :: tolf            !! stopping criterion if \( |f| < tolf \) then stop.
    real(wp),intent(in),optional      :: toldf           !! stopping criterion if \( |f_{n+1} - f_n| < toldf \) then stop
    real(wp),intent(in),optional      :: toldx           !! stopping criterion if \( ||x_{n+1} - x_n|| < toldx \) then stop
    integer,intent(in),optional       :: max_iter_ls     !! maximum number of iterations in the [[nnls]] problem
    integer,intent(in),optional       :: nnls_mode       !! Which NNLS method to use:
                                                         !!
                                                         !! 1. Use the original [[nnls]]
                                                         !! 2. Use the newer [[bvls]]
    real(wp),intent(in),optional :: infinite_bound !! "infinity" for the upper and lower bounds.
                                                   !! if `xl<=-infinite_bound` or `xu>=infinite_bound`
                                                   !! then these bounds are considered nonexistant.
                                                   !! If not present then `huge()` is used for this.
    procedure(msgfunc),optional  :: report_msg  !! user-defined procedure that will be called
                                                !! with any warning or error messages.

    integer :: n1,mineq,i

    status_ok = .false.
    call me%destroy()

    if (present(iprint)) me%iprint = iprint

    if (present(report_msg)) me%report_msg => report_msg

```