slsqpb Subroutine

    subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,&
                      r,l,x0,mu,s,u,v,w,&
                      t,f0,h1,h2,h3,h4,n1,n2,n3,t0,gs,tol,line,&
                      alpha,iexact,incons,ireset,itermx,ldat,&
                      alphamin,alphamax,tolf,toldf,toldx,&
                      max_iter_ls,nnls_mode,infBnd)
    implicit none

    integer,intent(in)                  :: m
    integer,intent(in)                  :: meq
    integer,intent(in)                  :: la
    integer,intent(in)                  :: n
    real(wp),dimension(n)               :: x
    real(wp),dimension(n)               :: xl
    real(wp),dimension(n)               :: xu
    real(wp)                            :: f
    real(wp),dimension(la)              :: c
    real(wp),dimension(n+1)             :: g
    real(wp),dimension(la,n+1)          :: a
    real(wp)                            :: acc
    integer,intent(inout)               :: iter     !! **in:**  maximum number of iterations.
                                                    !! **out:** actual number of iterations.
    integer,intent(inout)               :: mode
    real(wp),dimension(m+n+n+2)         :: r
    real(wp),dimension((n+1)*(n+2)/2)   :: l
    real(wp),dimension(n)               :: x0
    real(wp),dimension(la)              :: mu
    real(wp),dimension(n+1)             :: s
    real(wp),dimension(n+1)             :: u
    real(wp),dimension(n+1)             :: v
    real(wp),dimension(*),intent(inout) :: w   !! `dim(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]]
                                               !!
                                               !! with `mineq = m - meq + 2*n1` & `n1 = n+1`
    real(wp),intent(inout)              :: t
    real(wp),intent(inout)              :: f0
    real(wp),intent(inout)              :: h1
    real(wp),intent(inout)              :: h2
    real(wp),intent(inout)              :: h3
    real(wp),intent(inout)              :: h4
    integer,intent(inout)               :: n1
    integer,intent(inout)               :: n2
    integer,intent(inout)               :: n3
    real(wp),intent(inout)              :: t0
    real(wp),intent(inout)              :: gs
    real(wp),intent(inout)              :: tol
    integer,intent(inout)               :: line
    real(wp),intent(inout)              :: alpha
    integer,intent(inout)               :: iexact
    integer,intent(inout)               :: incons
    integer,intent(inout)               :: ireset
    integer,intent(inout)               :: itermx
    type(linmin_data),intent(inout)     :: ldat      !! data for [[linmin]].
    real(wp),intent(in)                 :: alphamin  !! min \( \alpha \) for line search
                                                     !! \( 0 < \alpha_{min} < \alpha_{max} \le 1 \)
    real(wp),intent(in)                 :: alphamax  !! max \( \alpha \) for line search
                                                     !! \( 0 < \alpha_{min} < \alpha_{max} \le 1 \)
    real(wp),intent(in)                 :: tolf      !! stopping criterion if \( |f| < tolf \) then stop.
    real(wp),intent(in)                 :: toldf     !! stopping criterion if \( |f_{n+1} - f_n| < toldf \) then stop
    real(wp),intent(in)                 :: toldx     !! stopping criterion if \( ||x_{n+1} - x_n|| < toldx \) then stop
    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.

    integer :: i, j, k
    logical :: inconsistent_linearization !! if the SQP problem is inconsistent

    inconsistent_linearization = .false. ! initialize

    if ( mode<0 ) then

        ! call jacobian at current x

        ! update cholesky-factors of hessian matrix by modified bfgs formula

        do i = 1 , n
            u(i) = g(i) - ddot(m,a(1,i),1,r,1) - v(i)
        end do

        ! l'*s

        k = 0
        do i = 1 , n
            h1 = zero
            k = k + 1
            do j = i + 1 , n
                k = k + 1
                h1 = h1 + l(k)*s(j)
            end do
            v(i) = s(i) + h1
        end do

        ! d*l'*s

        k = 1
        do i = 1 , n
            v(i) = l(k)*v(i)
            k = k + n1 - i
        end do

        ! l*d*l'*s

        do i = n , 1 , -1
            h1 = zero
            k = i
            do j = 1 , i - 1
                h1 = h1 + l(k)*v(j)
                k = k + n - j
            end do
            v(i) = v(i) + h1
        end do

        h1 = ddot(n,s,1,u,1)
        h2 = ddot(n,s,1,v,1)
        h3 = 0.2_wp*h2
        if ( h1<h3 ) then
            h4 = (h2-h3)/(h2-h1)
            h1 = h3
            call dscal(n,h4,u,1)
            call daxpy(n,one-h4,v,1,u,1)
        end if
        if (h1==zero .or. h2==zero) then
            ! Singular update: reset hessian.
            ! [ JW : this is based on a SciPy update ]
            call reset_bfgs_matrix()
            if ( ireset>5 ) return
        else
            call ldl(n,l,u,+one/h1,v)
            call ldl(n,l,v,-one/h2,u)
        end if
        ! end of main iteration

    else if ( mode==0 ) then

        itermx = iter
        if ( acc>=zero ) then
            iexact = 0
        else
            iexact = 1
        end if
        acc = abs(acc)
        tol = ten*acc
        iter = 0
        ireset = 0
        n1 = n + 1
        n2 = n1*n/2
        n3 = n2 + 1
        s(1) = zero
        mu(1) = zero
        call dcopy(n,s(1),0,s,1)
        call dcopy(m,mu(1),0,mu,1)

        call reset_bfgs_matrix()
        if ( ireset>5 ) return

    else

        ! call functions at current x

        t = f
        do j = 1 , m
            if ( j<=meq ) then
                h1 = c(j)
            else
                h1 = zero
            end if
            t = t + mu(j)*max(-c(j),h1)
        end do
        h1 = t - t0
        select case (iexact)
        case(0)
            if ( h1<=h3/ten .or. line>10 ) then
                call convergence_check(acc,0,-1)
            else
                alpha = min(max(h3/(two*(h3-h1)),alphamin),alphamax)
                call inexact_linesearch()
            end if
        case(1)
            call exact_linesearch()
            if ( line==3 ) call convergence_check(acc,0,-1)
        end select

        return

    end if

    do
        ! main iteration : search direction, steplength, ldl'-update

        iter = iter + 1
        mode = 9
        if ( iter>itermx ) return

        ! search direction as solution of qp - subproblem

        call dcopy(n,xl,1,u,1)
        call dcopy(n,xu,1,v,1)
        call daxpy(n,-one,x,1,u,1)
        call daxpy(n,-one,x,1,v,1)
        h4 = one
        call lsq(m,meq,n,n3,la,l,g,a,c,u,v,s,r,w,mode,max_iter_ls,nnls_mode,infBnd)

        ! augmented problem for inconsistent linearization

        inconsistent_linearization = .false. ! initialize
        if ( mode==6 ) then
            if ( n==meq ) mode = 4
        end if
        if ( mode==4 ) then
            ! Will reject this iteration if the SQP problem is inconsistent,
            inconsistent_linearization = .true.
            do j = 1 , m
                if ( j<=meq ) then
                    a(j,n1) = -c(j)
                else
                    a(j,n1) = max(-c(j),zero)
                end if
            end do
            s(1) = zero
            call dcopy(n,s(1),0,s,1)
            h3 = zero
            g(n1) = zero
            l(n3) = hun
            s(n1) = one
            u(n1) = zero
            v(n1) = one
            incons = 0
            do
                call lsq(m,meq,n1,n3,la,l,g,a,c,u,v,s,r,w,mode,max_iter_ls,nnls_mode,infBnd)
                h4 = one - s(n1)
                if ( mode==4 ) then
                    l(n3) = ten*l(n3)
                    incons = incons + 1
                    if ( incons<=5 ) cycle
                    return
                else if ( mode/=1 ) then
                    return
                else
                    exit
                end if
            end do

        else if ( mode/=1 ) then
            return
        end if

        ! update multipliers for l1-test

        do i = 1 , n
            v(i) = g(i) - ddot(m,a(1,i),1,r,1)
        end do
        f0 = f
        call dcopy(n,x,1,x0,1)
        gs = ddot(n,g,1,s,1)
        h1 = abs(gs)
        h2 = zero
        do j = 1 , m
            if ( j<=meq ) then
                h3 = c(j)
            else
                h3 = zero
            end if
            h2 = h2 + max(-c(j),h3)
            h3 = abs(r(j))
            mu(j) = max(h3,(mu(j)+h3)/two)
            h1 = h1 + h3*abs(c(j))
        end do

        ! check convergence

        mode = 0
        if ( h1<acc .and. h2<acc .and. &
             .not. inconsistent_linearization .and. &
             .not. ieee_is_nan(f)) return
        h1 = zero
        do j = 1 , m
            if ( j<=meq ) then
                h3 = c(j)
            else
                h3 = zero
            end if
            h1 = h1 + mu(j)*max(-c(j),h3)
        end do
        t0 = f + h1
        h3 = gs - h1*h4
        mode = 8
        if ( h3>=zero ) then
            call reset_bfgs_matrix()
            if ( ireset>5 ) return
        else
            exit
        end if

    end do

    ! line search with an l1-testfunction

    line = 0
    alpha = alphamax
    if ( iexact==1 ) then
        call exact_linesearch()
        if ( line==3 ) call convergence_check(acc,0,-1)
    else
        call inexact_linesearch()
    end if

    contains

        subroutine reset_bfgs_matrix()  ! 100
            !! reset BFGS matrix
            ireset = ireset + 1
            if ( ireset>5 ) then
                ! check relaxed convergence in case of positive directional derivative
                ! [ JW: reuse this routine so that h3 is recomputed.
                !   this is based on a SciPy update to SLSQP ]
                call convergence_check(tol,0,8)
                ! the caller should return in this case
            else
                l(1) = zero
                call dcopy(n2,l(1),0,l,1)
                j = 1
                do i = 1 , n
                    l(j) = one
                    j = j + n1 - i
                end do
            end if
        end subroutine reset_bfgs_matrix

        subroutine inexact_linesearch()  ! 300
            line = line + 1
            h3 = alpha*h3
            call dscal(n,alpha,s,1)
            call dcopy(n,x0,1,x,1)
            call daxpy(n,one,s,1,x,1)
            call enforce_bounds(x,xl,xu,infBnd)  ! ensure that x doesn't violate bounds
            mode = 1
        end subroutine inexact_linesearch

        subroutine exact_linesearch() ! 400
            if ( line/=3 ) then
                alpha = linmin(line,alphamin,alphamax,t,tol, &
                               ldat%a, ldat%b, ldat%d, ldat%e, ldat%p,   ldat%q,   &
                               ldat%r, ldat%u, ldat%v, ldat%w, ldat%x,   ldat%m,   &
                               ldat%fu,ldat%fv,ldat%fw,ldat%fx,ldat%tol1,ldat%tol2 )
                call dcopy(n,x0,1,x,1)
                call daxpy(n,alpha,s,1,x,1)
                mode = 1
            else
                call dscal(n,alpha,s,1)
            end if
        end subroutine exact_linesearch

        subroutine convergence_check(tolerance,converged,not_converged)  ! 500
            real(wp),intent(in) :: tolerance !! tolerance
            integer,intent(in) :: converged     !! mode value if converged
            integer,intent(in) :: not_converged !! mode value if not converged

            h3 = zero
            do j = 1 , m
                if (j<=meq) then
                    h1 = c(j)
                else
                    h1 = zero
                end if
                h3 = h3 + max(-c(j),h1)
            end do
            mode = check_convergence(n,f,f0,x,x0,s,h3,tolerance,tolf,toldf,toldx,&
                                     converged,not_converged,inconsistent_linearization)
        end subroutine convergence_check

    end subroutine slsqpb