dual_range_space_quad_prog Subroutine

private subroutine dual_range_space_quad_prog(me, nf, nc, x, ix, xl, xu, cf, cfd, ic, ica, cl, cu, cg, cr, cz, g, go, h, s, mfp, kbf, kbc, idecf, eta2, eta9, eps7, eps9, umax, gmax, n, iterq)

dual range space quadratic programming method.

Note

This routine was formerly called plqdb1.

Type Bound

psqp_class

Arguments

Type IntentOptional Attributes Name
class(psqp_class), intent(inout) :: me
integer :: nf

number of variables.
integer :: nc

number of linear constraints.
real(kind=wp) :: x(*)

x(nf) vector of variables.
integer :: ix(*)

ix(nf) vector containing types of bounds.
real(kind=wp) :: xl(*)

xl(nf) vector containing lower bounds for variables.
real(kind=wp) :: xu(*)

xu(nf) vector containing upper bounds for variables.
real(kind=wp) :: cf(*)

cf(nf) vector containing values of the constraint functions.
real(kind=wp) :: cfd(*)

cfd(nc) vector containing increments of the constraint functions.
integer :: ic(*)

ic(nc) vector containing types of constraints.
integer :: ica(*)

ica(nf) vector containing indices of active constraints.
real(kind=wp) :: cl(*)

cl(nc) vector containing lower bounds for constraint functions.
real(kind=wp) :: cu(*)

cu(nc) vector containing upper bounds for constraint functions.
real(kind=wp) :: cg(*)

cg(nf*nc) matrix whose columns are normals of the linear constraints.
real(kind=wp) :: cr(*)

cr(nf*(nf+1)/2) triangular decomposition of kernel of the orthogonal projection.
real(kind=wp) :: cz(*)

cz(nf) vector of lagrange multipliers.
real(kind=wp) :: g(*)

g(nf) gradient of the lagrangian function.
real(kind=wp) :: go(*)

go(nf) saved gradient of the objective function.
real(kind=wp) :: h(*)

h(nf*(nf+1)/2) triangular decomposition or inversion of the hessian matrix approximation.
real(kind=wp) :: s(*)

s(nf) direction vector.
integer :: mfp

type of feasible point.

  • mfp=1-arbitrary feasible point.
  • mfp=2-optimum feasible point.
  • mfp=3-repeated solution.
integer :: kbf

specification of simple bounds.

  • kbf=0-no simple bounds.
  • kbf=1-one sided simple bounds.
  • kbf=2=two sided simple bounds.
integer :: kbc

specification of linear constraints.

  • kbc=0 - no linear constraints.
  • kbc=1 - one sided linear constraints.
  • kbc=2 - two sided linear constraints.
integer :: idecf

decomposition indicator.

  • idecf=0 - no decomposition.
  • idecf=1 - gill-murray decomposition.
  • idecf=9 - inversion.
  • idecf=10 - diagonal matrix.
real(kind=wp) :: eta2

tolerance for positive definiteness of the hessian matrix.
real(kind=wp) :: eta9

maximum for real numbers.
real(kind=wp) :: eps7

tolerance for linear independence of constraints.
real(kind=wp) :: eps9

tolerance for activity of constraints.
real(kind=wp) :: umax

maximum absolute value of a negative lagrange multiplier.
real(kind=wp) :: gmax

maximum absolute value of a partial derivative.
integer :: n

dimension of the manifold defined by active constraints.
integer :: iterq

type of feasible point.

  • iterq=1 - arbitrary feasible point.
  • iterq=2 - optimum feasible point.
  • iterq=-1 - feasible point does not exists.
  • iterq=-2 - optimum feasible point does not exists.

Calls

proc~~dual_range_space_quad_prog~~CallsGraph proc~dual_range_space_quad_prog psqp_class%dual_range_space_quad_prog proc~determine_new_active_linear_constr determine_new_active_linear_constr proc~dual_range_space_quad_prog->proc~determine_new_active_linear_constr proc~determine_new_active_simple_bound determine_new_active_simple_bound proc~dual_range_space_quad_prog->proc~determine_new_active_simple_bound proc~mxdpgb mxdpgb proc~dual_range_space_quad_prog->proc~mxdpgb proc~mxdpgf mxdpgf proc~dual_range_space_quad_prog->proc~mxdpgf proc~mxdprb mxdprb proc~dual_range_space_quad_prog->proc~mxdprb proc~mxdsmm mxdsmm proc~dual_range_space_quad_prog->proc~mxdsmm proc~mxvcop mxvcop proc~dual_range_space_quad_prog->proc~mxvcop proc~mxvdir mxvdir proc~dual_range_space_quad_prog->proc~mxvdir proc~mxvina mxvina proc~dual_range_space_quad_prog->proc~mxvina proc~mxvinv mxvinv proc~dual_range_space_quad_prog->proc~mxvinv proc~mxvneg mxvneg proc~dual_range_space_quad_prog->proc~mxvneg proc~ops_after_constr_deletion psqp_class%ops_after_constr_deletion proc~dual_range_space_quad_prog->proc~ops_after_constr_deletion proc~update_tri_decomp_general update_tri_decomp_general proc~dual_range_space_quad_prog->proc~update_tri_decomp_general proc~mxvdot mxvdot proc~determine_new_active_linear_constr->proc~mxvdot proc~update_tri_decomp_orthogonal update_tri_decomp_orthogonal proc~ops_after_constr_deletion->proc~update_tri_decomp_orthogonal proc~update_tri_decomp_general->proc~mxdpgb proc~update_tri_decomp_general->proc~mxdprb proc~update_tri_decomp_general->proc~mxdsmm proc~update_tri_decomp_general->proc~mxvcop proc~mxdsmv mxdsmv proc~update_tri_decomp_general->proc~mxdsmv proc~update_tri_decomp_general->proc~mxvdot proc~mxvset mxvset proc~update_tri_decomp_general->proc~mxvset proc~update_tri_decomp_orthogonal->proc~mxvcop proc~update_tri_decomp_orthogonal->proc~mxvset proc~mxvort mxvort proc~update_tri_decomp_orthogonal->proc~mxvort proc~mxvrot mxvrot proc~update_tri_decomp_orthogonal->proc~mxvrot

Called by

proc~~dual_range_space_quad_prog~~CalledByGraph proc~dual_range_space_quad_prog psqp_class%dual_range_space_quad_prog proc~psqp psqp_class%psqp proc~psqp->proc~dual_range_space_quad_prog proc~psqpn psqp_class%psqpn proc~psqpn->proc~psqp

Source Code

   subroutine dual_range_space_quad_prog(me, nf, nc, x, ix, xl, xu, cf, cfd, &
                                         ic, ica, cl, cu, cg, cr, cz, g, go, h, s, &
                                         mfp, kbf, kbc, idecf, &
                                         eta2, eta9, eps7, eps9, umax, gmax, n, iterq)

      class(psqp_class), intent(inout) :: me
      integer :: nf     !! number of variables.
      integer :: nc     !! number of linear constraints.
      integer :: ix(*)  !! ix(nf)  vector containing types of bounds.
      integer :: ic(*)  !! ic(nc)  vector containing types of constraints.
      integer :: ica(*) !! ica(nf)  vector containing indices of active constraints.
      integer :: mfp    !! type of feasible point.
                        !!
                        !! * mfp=1-arbitrary feasible point.
                        !! * mfp=2-optimum feasible point.
                        !! * mfp=3-repeated solution.
      integer :: kbf    !! specification of simple bounds.
                        !!
                        !! * kbf=0-no simple bounds.
                        !! * kbf=1-one sided simple bounds.
                        !! * kbf=2=two sided simple bounds.
      integer :: kbc    !! specification of linear constraints.
                        !!
                        !! * kbc=0 - no linear constraints.
                        !! * kbc=1 - one sided linear constraints.
                        !! * kbc=2 - two sided linear constraints.
      integer :: idecf  !! decomposition indicator.
                        !!
                        !! * idecf=0  - no decomposition.
                        !! * idecf=1  - gill-murray decomposition.
                        !! * idecf=9  - inversion.
                        !! * idecf=10 - diagonal matrix.
      integer :: n      !! dimension of the manifold defined by active constraints.
      integer :: iterq  !! type of feasible point.
                        !!
                        !! * iterq=1  - arbitrary feasible point.
                        !! * iterq=2  - optimum feasible point.
                        !! * iterq=-1 - feasible point does not exists.
                        !! * iterq=-2 - optimum feasible point does not exists.
      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) :: cf(*) !! cf(nf)  vector containing 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) :: cg(*) !! cg(nf*nc)  matrix whose columns are normals of the linear constraints.
      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 lagrangian function.
      real(wp) :: go(*) !! go(nf)  saved gradient of the objective function.
      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) :: eta2  !! tolerance for positive definiteness of the hessian matrix.
      real(wp) :: eta9  !! maximum for real numbers.
      real(wp) :: eps7  !! tolerance for linear independence of constraints.
      real(wp) :: eps9  !! tolerance for activity of constraints.
      real(wp) :: umax  !! maximum absolute value of a negative lagrange multiplier.
      real(wp) :: gmax  !! maximum absolute value of a partial derivative.

      real(wp) :: con, temp, step, step1, step2, dmax, par, snorm
      integer :: nca, ncr, i, j, k, iold, jold, inew, jnew, knew, &
                 inf, ier, krem, kc, nred

      con = eta9
      if (idecf < 0) idecf = 1
      if (idecf == 0) then
         ! gill-murray decomposition
         temp = eta2
         call mxdpgf(nf, h, inf, temp, step)
         me%ndec = me%ndec + 1
         idecf = 1
      end if
      if (idecf >= 2 .and. idecf <= 8) then
         iterq = -10
         return
      end if

      ! initiation

      nred = 0
      jold = 0
      jnew = 0
      iterq = 0
      dmax = 0.0_wp
      if (mfp /= 3) then
         n = nf
         nca = 0
         ncr = 0
         if (kbf > 0) call mxvina(nf, ix)
         if (kbc > 0) call mxvina(nc, ic)
      end if

      outer: do

         ! direction determination

         call mxvneg(nf, go, s)
         do j = 1, nca
            kc = ica(j)
            if (kc > 0) then
               call mxvdir(nf, cz(j), cg((kc - 1)*nf + 1), s, s)
            else
               k = -kc
               s(k) = s(k) + cz(j)
            end if
         end do
         call mxvcop(nf, s, g)
         if (idecf == 1) then
            call mxdpgb(nf, h, s, 0)
         else
            call mxdsmm(nf, h, g, s)
         end if
         if (iterq /= 3) then
            ! check of feasibility
            inew = 0
            par = 0.0_wp
            call determine_new_active_linear_constr(nf, nc, cf, cfd, ic, cl, cu, &
                                                    cg, s, eps9, par, kbc, inew, knew)
            call determine_new_active_simple_bound(nf, ix, x, xl, xu, s, kbf, inew, &
                                                   knew, eps9, par)
            if (inew == 0) then
               ! solution achieved
               call mxvneg(nf, g, g)
               iterq = 2
               return
            else
               snorm = 0.0_wp
            end if

            inner: do

               ier = 0

               ! stepsize determination

               call update_tri_decomp_general(nf, n, ica, cg, cr, h, s, g, eps7, gmax, umax, &
                                              idecf, inew, me%nadd, ier, 1)
               call mxdprb(nca, cr, g, -1)
               if (knew < 0) call mxvneg(nca, g, g)

               ! primal stepsize

               if (ier /= 0) then
                  step1 = con
               else
                  step1 = -par/umax
               end if

               ! dual stepsize

               iold = 0
               step2 = con
               do j = 1, nca
                  kc = ica(j)
                  if (kc >= 0) then
                     k = ic(kc)
                  else
                     i = -kc
                     k = ix(i)
                  end if
                  if (k <= -5) then
                  elseif ((k == -1 .or. k == -3.) .and. g(j) <= 0.0_wp) then
                  elseif (.not. ((k == -2 .or. k == -4.) .and. g(j) >= 0.0_wp)) then
                     temp = cz(j)/g(j)
                     if (step2 > temp) then
                        iold = j
                        step2 = temp
                     end if
                  end if
               end do

               ! final stepsize

               step = min(step1, step2)
               if (step >= con) then
                  ! feasible solution does not exist
                  iterq = -1
                  return
               end if

               ! new lagrange multipliers

               dmax = step
               call mxvdir(nca, -step, g, cz, cz)
               snorm = snorm + sign(1, knew)*step
               par = par - (step/step1)*par
               if (step == step1) then
                  if (n <= 0) then
                     ! impossible situation
                     iterq = -5
                     return
                  end if

                  ! constraint addition

                  if (ier == 0) then
                     n = n - 1
                     nca = nca + 1
                     ncr = ncr + nca
                     cz(nca) = snorm
                  end if
                  if (inew > 0) then
                     kc = inew
                     call mxvinv(ic, kc, knew)
                  elseif (abs(knew) == 1) then
                     i = -inew
                     call mxvinv(ix, i, knew)
                  else
                     i = -inew
                     if (knew > 0) ix(i) = -3
                     if (knew < 0) ix(i) = -4
                  end if
                  nred = nred + 1
                  me%nadd = me%nadd + 1
                  jnew = inew
                  jold = 0
                  cycle outer
               end if

               ! constraint deletion

               do j = iold, nca - 1
                  cz(j) = cz(j + 1)
               end do
               call me%ops_after_constr_deletion(nf, nc, ix, ic, ica, cr, ic, g, n, iold, krem, ier)
               ncr = ncr - nca
               nca = nca - 1
               jold = iold
               jnew = 0
               if (kbc > 0) call mxvina(nc, ic)
               if (kbf > 0) call mxvina(nf, ix)
               do j = 1, nca
                  kc = ica(j)
                  if (kc > 0) then
                     ic(kc) = -ic(kc)
                  else
                     kc = -kc
                     ix(kc) = -ix(kc)
                  end if
               end do

            end do inner

         end if

         exit outer
      end do outer

   end subroutine dual_range_space_quad_prog