filterSD Subroutine

subroutine filterSD(n, m, x, al, f, fmin, cstype, bl, bu, ws, lws, v, nv, maxa, maxla, maxu, maxiu, kmax, maxg, rho, htol, rgtol, maxit, iprint, nout, ifail)

Copyright (C) 2010 Roger Fletcher

Current version dated 5 October 2011

THE ACCOMPANYING PROGRAM IS PROVIDED UNDER THE TERMS OF THE ECLIPSE PUBLIC LICENSE ("AGREEMENT"). ANY USE, REPRODUCTION OR DISTRIBUTION OF THE PROGRAM CONSTITUTES RECIPIENT'S ACCEPTANCE OF THIS AGREEMENT

Solves an NLP problem of the form: find a local solution x to

  minimize    f(x)
                    [  x   ]
  subject to  bl <= [      ] <= bu
                    [ c(x) ]

where f(x) is a given function of n variables x and c(x) is a vector of m given constraint functions. f(x) is minimized subject to lower and upper bounds bl and bu on x and c(x). f(x) and c(x) are defined by a user supplied subroutine 'functions'. The user is also required to supply a subroutine 'gradients' which calculates gradients of f(x) and c(x) with respect to x.

If linearized constraints and the trust region are incompatible, the code enters 'phase 1' in which an 'l1 feasibility problem' is solved. If this is unsuccessful in resolving the situation then the code exits with ifail=3 and returns a 'locally infeasible point' in x.

Parameter List (variables in a line starting with C must be set on entry)

n number of variables m number of general constraints x(n+m) x(1:n) stores the vector of variables. Initially an estimate of the solution must be set, replaced by the solution (if it exists) on exit. The rest of x is workspace al(n+m) stores Lagrange multipliers at the solution on exit. A positive multiplier indicates that the lower bound is active, and a negative multiplier indicates that the upper bound is active. Inactive constraints have a zero multiplier. f returns the value of f(x) when x is a feasible solution fmin set a strict lower bound on f(x) for feasible x (used to identify an unbounded NLP) cstype(m) character workspace: if ifail=3, cstype indicates constraints that are infeasible in the L1 solution. cstype(i)='A' if the lower bound on c/s i is infeasible, 'Z' if the upper bound is infeasible, else 'N' if feasible. bl(n+m) lower bounds on x and c(x) (use numbers no less than -ainfty (see below) and where possible supply realistic bounds on x) bu(n+m) upper bounds on x and c(x) (use numbers no greater than ainfty) ws() double precision workspace lws() integer workspace v(maxg) stores nv Ritz values (estimates of eigenvalues of reduced Hessian) supply the setting from a previous run of filterSD, or set nv=1 and v(1)=1.D0 in absence of other information nv number of values set in v maxa maximum number of entries in the Jacobian a() set by gradients maxla number of entries required for sparse matrix indices and pointers la(0:) to be set up in lws() (maxla>=maxa+m+3). Set maxla=1 if using dense matrix format maxu length of workspace user() passed through to user subroutines 'functions' and 'gradients' maxiu length of workspace iuser(*) passed through to user subroutines kmax maximum dimension of null space allowed for (kmax<=n) maxg maximum number of reduced gradient vectors stored by the limited memory method (typically 6 or 7) rho initial trust region radius (typically 1.D1) htol tolerance allowed in sum h of constraint feasibilities (e.g. 1.D-6) rgtol tolerance allowed in reduced gradient l2 norm (typically 1.D-4) maxit maximum number of major iterations allowed iprint verbosity of printing (0=none, 1=one line per iteration, 2=additional text information given) nout output channel for printing ifail returns failure indication as follows 0 = successful run 1 = unbounded NLP (f <= fmin at an htol-feasible point) 2 = bounds on x are inconsistent 3 = local minimum of feasibility problem and h > htol (nonlinear constraints are locally inconsistent) 4 = initial point x has h > ubd (reset ubd or x and re-enter) 5 = maxit major iterations have been carried out 6 = termination with rho <= htol 7 = not enough workspace in ws or lws (see message) 8 = insufficient space for filter (increase mxf and re-enter) >9 = unexpected fail in LCP solver (10 has been added to ifail)

User Routines

The user must supply two subroutines to calculate f(x), c(x) and their first derivatives as follows

     subroutine functions(n,m,x,f,c,user,iuser)
     implicit double precision (a-h,o-z)
     dimension x(*),c(*),user(*),iuser(*)
     !...
     !Statements to calculate f(x) and the m-vector c(x). The user is
     !responsible for ensuring that any failures such as IEEE errors
     !(overflow, NaN's etc.) are trapped and not returned to filterSD.
     !The same holds for gradients.
     !...
     return
     end

     subroutine gradients(n,m,x,a,user,iuser)
     implicit double precision(a-h,o-z)
     dimension x(*),a(*),user(*),iuser(*)
     !...
     !Statements to calculate gradients of f(x) and c(x) and set in a(*).
     !The column vector grad(f) must be followed by the column vectors
     !grad(c_i), i=1,2,...,m, in the one dimensional array a(*). Either a
     !dense or sparse data structure may be used. If using the sparse data
     !structure, only stucturally non-zero entries are set. Pointers etc. for
     !the data structure are set once and for all in lws as described below.
     !The user may assume that a call of 'gradients' immediately follows one
     ! of 'functions' with the same vector x.)
     !...
     return
     end

The user must also supply a driver routine which calls filterSD. This must set parameters and common blocks of filterSD as appropriate. Space for x,al,bl,bu,ws,lws,v and cstype must be assigned. If using the sparse data structure for setting gradients, indices and pointers la(0:maxla-1) must be set in the driver in lws, immediately following any user workspace in lws(1:maxiu). No changes in this data structure are allowed during the operation of filterSD. More details of the format of a() and la() are given in the file sparseA.f For dense format just set maxla=1 and set lws(maxiu+1) to the 'stride' (>=n) used in setting the columns of grad(f) and grad(c_i). For efficiency, constant entries in the gradients may be set in the driver. However two copies of the gradients are kept by filterSD . These reside in ws(maxu+1:maxu+maxa) and ws(maxu+maxa+1:maxu+2*maxa). Any constant entries must be set in both copies.

Arguments

    subroutine filterSD(n,m,x,al,f,fmin,cstype,bl,bu,ws,lws,v,nv, &
                        maxa,maxla,maxu,maxiu,kmax,maxg,rho,htol,rgtol,maxit,iprint, &
                        nout,ifail)

      implicit double precision (a-h,o-z)
      dimension x(*),al(*),bl(*),bu(*),ws(*),lws(*),v(*)
      character cstype(*)

!  Common blocks
!  =============
!      common/wsc/kk,ll,kkk,lll,mxws,mxlws
!  The user must specify the length of the workspace arrays ws(*) and lws(*)
!  in mxws and mxlws respectively. It may not be easy to specify a-priori how
!  large these arrays should be. Set a suitable large estimate, and filterSD
!  will prompt if larger values are required. As a guide, ws(*) contains
!  first user workspace, then workspace for filterSD, then workspace for glcpd,
!  and finally workspace for denseL.f or schurQR.f. lws(*) contains user
!  workspace, then maxla+n+m+mlp locations for filterSD and additional locations
!  for denseL.f or schurQR.f.
!     common/defaultc/ainfty,ubd,mlp,mxf
!  Default values of some control parameters are set here. ainfty is used to
!  represent infinity. ubd provides an upper bound on the allowed constraint
!  violation. mlp is the maximum length of arrays used in degeneracy control.
!  mxf is the maximum length of filter arrays. Default values are 1.D20, 1.D4,
!  50,  50 respectively.
!     common/ngrc/mxgr
!  The user can limit the time spent in each call of the LCP solver by setting
!  an upper limit on the number of gradient calls in mxgr (default=1000000)
!     common/mxm1c/mxm1
!  When using denseL.f, mxm1 must be set to the maximum number of general
!  constraints allowed in the active set. mxm1=min(m+1,n) is always sufficient
!     common/epsc/eps,tol,emin
!     common/repc/sgnf,nrep,npiv,nres
!  These common blocks provide default parameters that control glcpd
!     common/statsc/dnorm,h,hJt,hJ,ipeq,k,itn,nft,ngt
!  This common block returns information about the outcome of filterSD.
!  dnorm=final step length, h=final c/s violation, hJt=ditto for 'N' c/s,
!  hJ=ditto for 'A' and 'Z' c/s, ipeq=number of active equations, k=number of
!  free variables, itn=number of iterations, nft=total number of function
!  calls, ngt=total number of gradient calls

      common/wsc/kk,ll,kkk,lll,mxws,mxlws
      common/defaultc/ainfty,ubd,mlp,mxf
      common/functc/fxd,alc,m_,iph,last1,next1,nx,nx1, &
        nal,nal1,naal,naal1,nxd,nxd1,ncx,ncx1,ncxd,ncxd1,nla1

      m_=m
      nm=n+m
!  set real storage map for ws
!  first maxu locations are user storage for functions and gradients
!  vectors required by funct: two slots of length maxa for a(*)
      last1=maxu+1
      next1=last1+maxa
!  slot of length n for x
      nx1=next1+maxa
      nx=nx1-1
!  slot of length m for lambda
      nal=nx+n
      nal1=nal+1
!  slot of length n for Ak.al, or for storing x
      naal=nal+m
      naal1=naal+1
!  slot of length n for x at x+d
      nxd=naal+n
      nxd1=nxd+1
!  slot of length m for c at x
      ncx=nxd+n
      ncx1=ncx+1
!  slot of length m for c at x+d
      ncxd=ncx+m
      ncxd1=ncxd+1
!  local storage for filter_SD
!  slot of length n for d
      id1=ncxd1+m
!  slot of length n+m for dl
      idl1=id1+n
!  slot of length n+m for du
      idu1=idl1+nm
!  slot of length n for g
      ig1=idu1+nm
!  slot of length n+m for e
      ie1=ig1+n
!  slot of length mlp for alp
      ialp1=ie1+nm
!  slot of length mxf for filh
      ifilh1=ialp1+mlp
!  slot of length mxf for filf
      ifilf1=ifilh1+mxf
!  total length of ws so far
      kk=ifilf1+mxf-1

!  set integer storage map for lws
!  first maxiu locations are user storage for functions and gradients
!  storage of length maxla for la(0:*)
      nla1=maxiu+1
!  local storage for filter_SD
!  slot of length n+m for ls
      ils1=nla1+maxla
      ils=ils1-1
!  slot of length mlp for lp
      ilp1=ils1+nm
!  total length of lws so far
      ll=ilp1+mlp-1

      do i=1,n
        if (bl(i)>bu(i)) then
          if (iprint>1) write(nout,*)'simple bounds infeasible'
          ifail=2
          return
        end if
        ws(nx+i)=min(max(bl(i),x(i)),bu(i))
      end do

!  note x and al are just used as workspace: the true values are those in ws

      call filter_SD(n,f,fmin,cstype,bl,bu,ws,lws,v,nv, &
        maxa,kmax,maxg, &
        ws(id1),ws(idl1),ws(idu1),ws(ig1),x,al,ws(ie1),lws(ils1), &
        ws(ialp1),lws(ilp1),ws(ifilh1),ws(ifilf1),rho,htol,rgtol, &
        maxit,iprint,nout,ifail)

      if (ifail>=7) return
!  scatter ws(nx.. and ws(nal.. and bound multipliers into x and al
      do i=1,n
        al(i)=0.D0
      end do
      do i=1,m
        al(n+i)=ws(nal+i)
      end do
      do j=1,n
        i=abs(lws(ils+j))
        if (i<=n) then
          if (ws(nx+i)==bl(i)) then
            al(i)=x(i)
          else if (ws(nx+i)==bu(i)) then
            al(i)=-x(i)
          end if
        end if
      end do
      do i=1,n
        x(i)=ws(nx+i)
      end do

      return
1     format(A,15I5)
    end subroutine filterSD