psqp Subroutine

   subroutine psqp(me, nf, nb, nc, x, ix, xl, xu, cf, ic, cl, cu, cg, cfo, cfd, gc, ica, &
                   cr, cz, cp, gf, g, h, s, xo, go, xmax, tolx, tolc, tolg, rpf, &
                   cmax, gmax, f, mit, mfv, met, mec, iprnt, iterm)

      class(psqp_class), intent(inout) :: me
      real(wp) :: f         !! value of the objective function.
      real(wp) :: cmax      !! maximum constraint violation.
      real(wp) :: gmax      !! maximum partial derivative of the lagrangian function.
      real(wp) :: rpf       !! penalty coefficient.
      real(wp) :: tolx      !! tolerance for change of variables.
      real(wp) :: tolc      !! tolerance for constraint violations.
      real(wp) :: tolg      !! tolerance for the gradient of the lagrangian function.
      real(wp) :: told      !!
      real(wp) :: tols      !!
      real(wp) :: xmax      !! maximum stepsize.
      integer :: iprnt      !! print specification.
                            !!
                            !! * iprnt=0      - no print.
                            !! * abs(iprnt)=1 - print of final results.
                            !! * abs(iprnt)=2 - print of final results and iterations.
                            !! * iprnt>0      - basic final results.
                            !! * iprnt<0      - extended final results.
      integer :: iterm      !! variable that indicates the cause of termination.
                            !!
                            !! * iterm=1-if abs(x-xo) was less than or equal to tolx in
                            !!   mtesx (usually two) subsequent iterations.
                            !! * iterm=2-if abs(f-fo) was less than or equal to tolf in
                            !!   mtesf (usually two) subsequent iterations.
                            !! * iterm=3-if f is less than or equal to tolb.
                            !! * iterm=4-if gmax is less than or equal to tolg.
                            !! * iterm=11-if nit exceeded mit.
                            !! * iterm=12-if nfv exceeded mfv.
                            !! * iterm=13-if nfg exceeded mfg.
                            !! * iterm<0-if the method failed.
                            !! * if iterm=-6, then the termination criterion has not been
                            !!   satisfied, but the point obtained if usually acceptable.
      integer :: met        !! variable metric update used.
                            !!
                            !! * met=1 - the bfgs update.
                            !! * met=2 - the hoshino update.
      integer :: met1       !!
      integer :: mec        !! correction if the negative curvature occurs.
                            !!
                            !! * mec=1 - correction suppressed.
                            !! * mec=2 - powell's correction.
      integer :: mes        !!
      integer :: mfv        !! maximum number of function evaluations.
      integer :: mit        !! maximum number of iterations.
      integer :: nb         !! choice of simple bounds.
                            !!
                            !! * nb=0 - simple bounds suppressed.
                            !! * nb>0 - simple bounds accepted.
      integer :: nc         !! number of linear constraints.
      integer :: nf         !! number of variables.
      real(wp) :: cf(*)     !! cf(nc+1)  vector containing values of the constraint functions.
      real(wp) :: cg(*)     !! cg(nf*nc)  matrix whose columns are normals of the linear constraints.
      real(wp) :: cfo(*)    !! cfo(nc)  vector containing saved values of the constraint functions.
      real(wp) :: cfd(*)    !! cfd(nc)  vector containing increments of the constraint functions.
      real(wp) :: cl(*)     !! cl(nc)  vector containing lower bounds for constraint functions.
      real(wp) :: cu(*)     !! cu(nc)  vector containing upper bounds for constraint functions.
      real(wp) :: cp(*)     !!
      real(wp) :: cr(*)     !! cr(nf*(nf+1)/2)  triangular decomposition of kernel
                            !! of the orthogonal projection.
      real(wp) :: cz(*)     !! cz(nf)  vector of lagrange multipliers.
      real(wp) :: g(*)      !! g(nf)  gradient of the objective function.
      real(wp) :: gc(*)     !! gc(nf)  gradient of the selected constraint function.
      real(wp) :: gf(*)     !! gf(nf)  gradient of the model function.
      real(wp) :: go(*)     !! go(nf)  gradients difference.
      real(wp) :: h(*)      !! h(nf*(nf+1)/2)  triangular decomposition or inversion of
                            !! the hessian matrix approximation.
      real(wp) :: s(*)      !! s(nf)  direction vector.
      real(wp) :: x(*)      !! x(nf)  vector of variables.
      real(wp) :: xl(*)     !! xl(nf)  vector containing lower bounds for variables.
      real(wp) :: xu(*)     !! xu(nf)  vector containing upper bounds for variables.
      real(wp) :: xo(*)     !! xo(nf)  vectors of variables difference.
      integer,intent(inout) :: ic(*) !! ic(nc)  vector containing types of constraints.
                                     !!
                                     !! * ic(kc)=0 - constraint cf(kc) is not used.
                                     !! * ic(kc)=1 - lower constraint cl(kc)<=cf(kc).
                                     !! * ic(kc)=2 - upper constraint cf(kc)<=cu(kc).
                                     !! * ic(kc)=3 - two side constraint cl(kc)<=cf(kc)<=cu(kc).
                                     !! * ic(kc)=5 - equality constraint cf(kc)==cl(kc).
      integer :: ica(*)     !! ica(nf)  vector containing indices of active constraints.
      integer,intent(inout) :: ix(*) !! ix(nf)  vector containing types of bounds.
                                     !!
                                     !! * ix(i)=0 - variable x(i) is unbounded.
                                     !! * ix(i)=1 - lover bound xl(i)<=x(i).
                                     !! * ix(i)=2 - upper bound x(i)<=xu(i).
                                     !! * ix(i)=3 - two side bound xl(i)<=x(i)<=xu(i).
                                     !! * ix(i)=5 - variable x(i) is fixed.

      real(wp) :: alf1, alf2, cmaxo, dmax, eps7, eps9, eta0, &
                  eta2, eta9, fmax, fmin, fo, gnorm, p, po, &
                  r, rmax, rmin, ro, snorm, tolb, tolf, &
                  umax, rp, fp, pp, ff, fc
      integer :: i, idecf, iext, irest, iterd, iterh, iterq, &
                 iters, kbc, kbf, kc, kd, kit, ld, mred, mtesf, &
                 mtesx, n, k, ntesx, iest, inits, kters, maxst, &
                 isys, mfp, nred, ipom, lds

      if (abs(iprnt) > 1) write (6, '(1x,"entry to psqp :")')

      ! initiation
      kbf = 0
      kbc = 0
      if (nb > 0) kbf = 2
      if (nc > 0) kbc = 2
      me%nres = 0
      me%ndec = 0
      me%nrem = 0
      me%nadd = 0
      me%nit = 0
      me%nfv = 0
      me%nfg = 0
      me%nfh = 0
      isys = 0
      iest = 0
      iext = 0
      mtesx = 2
      mtesf = 2
      inits = 1
      iterm = 0
      iters = 0
      iterd = 0
      iterq = 0
      mred = 20
      irest = 1
      iters = 2
      kters = 5
      idecf = 1
      eta0 = 1.0e-15_wp
      eta2 = 1.0e-15_wp
      eta9 = huge(1.0_wp) !1.0e60_wp
      eps7 = 1.0e-15_wp
      eps9 = 1.0e-8_wp
      alf1 = 1.0e-10_wp
      alf2 = 1.0e10_wp
      fmax = huge(1.0_wp) !1.0e60_wp
      fmin = -fmax
      tolb = -fmax
      dmax = eta9
      tolf = 1.0e-16_wp
      if (xmax <= 0.0_wp) xmax = 1.0e+16_wp
      if (tolx <= 0.0_wp) tolx = 1.0e-16_wp
      if (tolg <= 0.0_wp) tolg = 1.0e-6_wp
      if (tolc <= 0.0_wp) tolc = 1.0e-6_wp
      told = 1.0e-8_wp
      tols = 1.0e-4_wp
      if (rpf <= 0.0_wp) rpf = 1.0e-4_wp
      if (met <= 0) met = 1
      met1 = 2
      if (mec <= 0) mec = 2
      mes = 1
      if (mit <= 0) mit = 1000
      if (mfv <= 0) mfv = 2000
      kd = 1
      ld = -1
      kit = 0
      call mxvset(nc, 0.0_wp, cp)
      ! initial operations with simple bounds
      if (kbf > 0) then
         do i = 1, nf
            if ((ix(i) == 3 .or. ix(i) == 4) .and. xu(i) <= xl(i)) then
               xu(i) = xl(i)
               ix(i) = 5
            elseif (ix(i) == 5 .or. ix(i) == 6) then
               xl(i) = x(i)
               xu(i) = x(i)
               ix(i) = 5
            end if
            if (ix(i) == 1 .or. ix(i) == 3) x(i) = max(x(i), xl(i))
            if (ix(i) == 2 .or. ix(i) == 3) x(i) = min(x(i), xu(i))
         end do
      end if
      ! initial operations with general constraints
      if (kbc > 0) then
         k = 0
         do kc = 1, nc
            if ((ic(kc) == 3 .or. ic(kc) == 4) .and. cu(kc) <= cl(kc)) then
               cu(kc) = cl(kc)
               ic(kc) = 5
            elseif (ic(kc) == 5 .or. ic(kc) == 6) then
               cu(kc) = cl(kc)
               ic(kc) = 5
            end if
            k = k + nf
         end do
      end if
      if (kbf > 0) then
         do i = 1, nf
            if (ix(i) >= 5) ix(i) = -ix(i)
            if (ix(i) <= 0) then
            elseif ((ix(i) == 1 .or. ix(i) == 3) .and. x(i) <= xl(i)) then
               x(i) = xl(i)
            elseif ((ix(i) == 2 .or. ix(i) == 3) .and. x(i) >= xu(i)) then
               x(i) = xu(i)
            end if
            call test_simple_bound(x, ix, xl, xu, eps9, i)
            if (ix(i) > 10) ix(i) = 10 - ix(i)
         end do
      end if
      fo = fmin
      gmax = eta9
      dmax = eta9

      main: do

         lds = ld
         call me%compute_obj_and_dobj(nf, x, gf, gf, ff, f, kd, ld, iext)
         ld = lds
         call me%compute_con_and_dcon(nf, nc, x, fc, cf, cl, cu, ic, gc, cg, cmax, kd, ld)
         cf(nc + 1) = f
         ! JW : seems to start with iter 0, so add 1 to iterations for report output
         if (associated(me%report)) call me%report(me%nit+1,x(1:nf),j=f,f=cf(1:nc))    
         if (abs(iprnt) > 1) &
            write (6, '(1x,"nit=",i9,2x,"nfv=",i9,2x,"nfg=",i9,2x,"f=",g13.6,2x,"c=",e8.1,2x,"g=",e8.1)') &
            me%nit, me%nfv, me%nfg, f, cmax, gmax
         ! start of the iteration with tests for termination.
         if (iterm < 0) exit main
         if (iters /= 0) then
            if (f <= tolb) then
               iterm = 3
               exit main
            end if
            if (dmax <= tolx) then
               iterm = 1
               ntesx = ntesx + 1
               if (ntesx >= mtesx) exit main
            else
               ntesx = 0
            end if
         end if
         if (me%nit >= mit) then
            iterm = 11
            exit main
         end if
         if (me%nfv >= mfv) then
            iterm = 12
            exit main
         end if
         iterm = 0
         me%nit = me%nit + 1

         restart: do

            ! restart
            n = nf
            if (irest > 0) then
               call mxdsmi(n, h)
               ld = min(ld, 1)
               idecf = 1
               if (kit < me%nit) then
                  me%nres = me%nres + 1
                  kit = me%nit
               else
                  iterm = -10
                  if (iters < 0) iterm = iters - 5
                  exit main
               end if
            end if
            ! direction determination using a quadratic programming procedure
            call mxvcop(nc + 1, cf, cfo)
            mfp = 2
            ipom = 0
            dir_loop: do
               call me%dual_range_space_quad_prog(nf, nc, x, ix, xl, xu, cf, cfd, ic, ica, &
                                                  cl, cu, cg, cr, cz, g, gf, h, s, mfp, kbf, &
                                                  kbc, idecf, eta2, eta9, eps7, &
                                                  eps9, umax, gmax, n, iterq)
               if (iterq < 0) then
                  if (ipom < 10) then
                     ipom = ipom + 1
                     call transform_incompatible_qp_subproblem(nc, cf, ic, cl, cu, kbc)
                     cycle dir_loop
                  end if
                  iterd = iterq - 10
               else
                  ipom = 0
                  iterd = 1
                  gmax = mxvmax(nf, g)
                  gnorm = sqrt(mxvdot(nf, g, g))
                  snorm = sqrt(mxvdot(nf, s, s))
               end if
               exit dir_loop
            end do dir_loop
            if (iterd < 0) iterm = iterd
            if (iterm == 0) then
               call mxvcop(nc + 1, cfo, cf)
               ! test for sufficient descent
               p = mxvdot(nf, g, s)
               irest = 1
               if (snorm <= 0.0_wp) then
               elseif (p + told*gnorm*snorm <= 0.0_wp) then
                  irest = 0
               end if
               if (irest /= 0) cycle restart
               nred = 0
               rmin = alf1*gnorm/snorm
               rmax = min(alf2*gnorm/snorm, xmax/snorm)
               if (gmax <= tolg .and. cmax <= tolc) then
                  iterm = 4
                  exit main
               end if
               call compute_new_penalty_parameters(nf, n, nc, ica, cz, cp)
               call mxvina(nc, ic)
               call compute_augmented_lagrangian(nf, n, nc, cf, ic, ica, cl, cu, cz, rpf, fc, f)
               ! preparation of line search
               ro = 0.0_wp
               fo = f
               po = p
               cmaxo = cmax
               call mxvcop(nf, x, xo)
               call mxvcop(nf, g, go)
               call mxvcop(nf, gf, cr)
               call mxvcop(nc + 1, cf, cfo)

               line_search: do

                  ! line search without directional derivatives
                  call me%extended_line_search(r, ro, rp, f, fo, fp, po, pp, fmin, fmax, &
                                               rmin, rmax, tols, kd, ld, me%nit, kit, nred, &
                                               mred, maxst, iest, inits, iters, kters, mes, isys)
                  if (isys == 0) then
                     kd = 1
                     ! decision after unsuccessful line search
                     if (iters <= 0) then
                        r = 0.0_wp
                        f = fo
                        p = po
                        call mxvcop(nf, xo, x)
                        call mxvcop(nf, cr, gf)
                        call mxvcop(nc + 1, cfo, cf)
                        irest = 1
                        ld = kd
                        cycle restart
                     end if
                     ! computation of the value and the gradient of the objective
                     ! function together with the values and the gradients of the
                     ! approximated functions
                     if (kd > ld) then
                        lds = ld
                        call me%compute_obj_and_dobj(nf, x, gf, gf, ff, f, kd, ld, iext)
                        ld = lds
                        call me%compute_con_and_dcon(nf, nc, x, fc, cf, cl, cu, ic, gc, &
                                                     cg, cmax, kd, ld)
                     end if
                     ! preparation of variable metric update
                     call mxvcop(nf, gf, g)
                     call dual_range_space_qp(nf, n, x, xo, ica, cg, cz, g, go, r, f, fo, p, po, &
                                              cmax, cmaxo, dmax, kd, ld, iters)
                     ! variable metric update
                     call bfgs_variable_metric_update(n, h, g, s, xo, go, r, po, me%nit, &
                                                      kit, iterh, met, met1, mec)
                     ! if (mer>0.and.iterh>0) irest=1
                     cycle main   ! end of the iteration
                  else
                     ! go to (11174,11172) isys+1
                     call mxvdir(nf, r, s, xo, x)
                     lds = ld
                     call me%compute_obj_and_dobj(nf, x, gf, g, ff, f, kd, ld, iext)
                     ld = lds
                     call me%compute_con_and_dcon(nf, nc, x, fc, cf, cl, cu, ic, gc, cg, cmax, kd, ld)
                     cf(nc + 1) = f
                     call compute_augmented_lagrangian(nf, n, nc, cf, ic, ica, cl, cu, cz, rpf, fc, f)
                  end if

               end do line_search

            end if

            exit restart
         end do restart

         exit main
      end do main

      if (iprnt > 1 .or. iprnt < 0) write (6, '(1x,"exit from psqp :")')
      if (iprnt /= 0) &
         write (6, '(1x,"nit=",i4,2x,"nfv=",i4,2x,"nfg=",i4,2x,"f=",g13.6,2x,"c=",e8.1,2x,"g=",e8.1,2x,"iterm=",i3)') &
                  me%nit, me%nfv, me%nfg, f, cmax, gmax, iterm
      if (iprnt < 0) write (6, '(1x,"x=",5(g14.7,1x):/(3x,5(g14.7,1x)))') (x(i), i=1, nf)

   end subroutine psqp