diff --git a/.travis.yml b/.travis.yml index a09fca6..d8bf71d 100644 --- a/.travis.yml +++ b/.travis.yml @@ -24,6 +24,7 @@ script: - ./bin/slsqp_test - ./bin/slsqp_test_2 - ./bin/slsqp_test_3 + - ./bin/slsqp_test_71 after_success: - git config --global user.name "TRAVIS-CI-for-$(git --no-pager show -s --format='%cn' $TRAVIS_COMMIT)" diff --git a/LICENSE b/LICENSE index f21ce4e..ea2b6aa 100644 --- a/LICENSE +++ b/LICENSE @@ -1,7 +1,7 @@ Modern Fortran Edition of the SLSQP Optimizer https://github.com/jacobwilliams/slsqp -Copyright (c) 2016-2018, Jacob Williams +Copyright (c) 2016-2019, Jacob Williams All rights reserved. Redistribution and use in source and binary forms, with or without modification, diff --git a/README.md b/README.md index 14edbcb..bd22dfa 100644 --- a/README.md +++ b/README.md @@ -38,6 +38,8 @@ A [FoBiS](https://github.com/szaghi/FoBiS) configuration file (`slsqp.fobis`) is To generate the documentation using [ford](https://github.com/cmacmackin/ford), run: ```FoBis.py rule --execute makedoc -f slsqp.fobis``` + To run the test programs, run: ```FoBis.py rule --execute tests -f slsqp.fobis``` + ### Development [![Build Status](https://img.shields.io/travis/jacobwilliams/slsqp/master.svg?style=plastic)](https://travis-ci.org/jacobwilliams/slsqp) diff --git a/build.sh b/build.sh new file mode 100755 index 0000000..e96bfd8 --- /dev/null +++ b/build.sh @@ -0,0 +1,52 @@ +#!/bin/bash + +# +# Simple build script. +# +# Requires: FoBiS and Ford +# + +MODCODE='slsqp_module.f90' # module file name +LIBOUT='libslsqp.a' # name of library +DOCDIR='./doc/' # build directory for documentation +SRCDIR='./src/' # library source directory +TESTSRCDIR='./src/tests/' # unit test source directory +BINDIR='./bin/' # build directory for unit tests +LIBDIR='./lib/' # build directory for library +FORDMD='slsqp.md' # FORD config file name + +#compiler flags: + +#FCOMPILER='gnu' #Set compiler to gfortran +#FCOMPILERFLAGS='-c -O2 -std=f2008' + +FCOMPILER='intel' #Set compiler to ifort +FCOMPILERFLAGS='-c -O2' + +#build using FoBiS: + +if hash FoBiS.py 2>/dev/null; then + + echo "Building library..." + + FoBiS.py build -compiler ${FCOMPILER} -cflags "${FCOMPILERFLAGS}" -dbld ${LIBDIR} -s ${SRCDIR} -dmod ./ -dobj ./ -t ${MODCODE} -o ${LIBOUT} -mklib static -colors + + echo "Building test programs..." + + FoBiS.py build -compiler ${FCOMPILER} -cflags "${FCOMPILERFLAGS}" -dbld ${BINDIR} -s ${TESTSRCDIR} -dmod ./ -dobj ./ -colors -libs ${LIBDIR}${LIBOUT} --include ${LIBDIR} + +else + echo "FoBiS.py not found! Cannot build library. Install using: sudo pip install FoBiS.py" +fi + +# build the documentation using FORD: + +if hash ford 2>/dev/null; then + + echo "Building documentation..." + + ford ${FORDMD} + +else + echo "Ford not found! Cannot build documentation. Install using: sudo pip install ford" +fi diff --git a/slsqp.code-workspace b/slsqp.code-workspace new file mode 100644 index 0000000..489f23c --- /dev/null +++ b/slsqp.code-workspace @@ -0,0 +1,13 @@ +{ + "folders": [ + { + "path": "." + } + ], + "settings": { + "files.trimTrailingWhitespace": true, + "editor.insertSpaces": true, + "editor.tabSize": 4, + "editor.trimAutoWhitespace": true + } +} \ No newline at end of file diff --git a/slsqp.fobis b/slsqp.fobis index cba6803..e995ddf 100644 --- a/slsqp.fobis +++ b/slsqp.fobis @@ -142,3 +142,15 @@ template = template-tests [rule-makedoc] help = Rule for building documentation from source files rule_1 = ford $FORD_FILE + +[rule-tests] +help = Rule for running the test programs +rule_1 = (cd bin + GLOBIGNORE='*.*' + ls slsqp_test_* | sed 's/^\([^0-9]*\)\([0-9]*\)/\1 \2/' | sort -k2,2n | tr -d ' ' | + while read TEST; do + echo "" + echo "Running ${TEST}" + "./${TEST}" + done) + diff --git a/src/slsqp_core.f90 b/src/slsqp_core.f90 index 5048afa..e539822 100644 --- a/src/slsqp_core.f90 +++ b/src/slsqp_core.f90 @@ -110,7 +110,7 @@ module slsqp_core !@note `f`, `c`, `g`, `a` must all be set by the user before each call. subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & - jw,l_jw,sdat,ldat,alphamin,alphamax) + sdat,ldat,alphamin,alphamax,tolf,toldf,toldx) implicit none @@ -167,9 +167,8 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & !! * ** 7 **: rank-deficient equality constraint subproblem [[hfti]] !! * ** 8 **: positive directional derivative for linesearch !! * ** 9 **: more than `iter` iterations in sqp - !! * ** >=10 **: working space `w` or `jw` too small, + !! * ** >=10 **: working space `w` too small, !! `w` should be enlarged to `l_w=mode/1000`, - !! `jw` should be enlarged to `l_jw=mode-1000*l_w` integer,intent(in) :: l_w !! the length of `w`, which should be at least: !! !! * `(3*n1+m)*(n1+1)` **for lsq** @@ -178,8 +177,6 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & !! * `+ n1*n/2 + 2*m + 3*n + 3*n1 + 1` **for slsqpb** !! !! with `mineq = m - meq + 2*n1` & `n1 = n+1` - integer,intent(in) :: l_jw !! the length of `jw` which should be at least - !! `mineq = m - meq + 2*(n+1)`. real(wp),dimension(l_w),intent(inout) :: w !! `w()` is a one dimensional working space. !! the first `m+n+n*n1/2` elements of `w` must not be !! changed between subsequent calls of [[slsqp]]. @@ -193,11 +190,15 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & !! contain the multipliers associated with all !! constraints of the quadratic program finding !! the search direction to the solution `x*` - integer,dimension(l_jw),intent(inout) :: jw !! `jw()` is a one dimensional integer working space type(slsqpb_data),intent(inout) :: sdat !! data for [[slsqpb]]. 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) :: 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 :: il , im , ir , is , iu , iv , iw , ix , mineq, n1 @@ -207,15 +208,17 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & il = (3*n1+m)*(n1+1) + (n1-meq+1)*(mineq+2) + 2*mineq + (n1+mineq)& *(n1-meq) + 2*meq + n1*n/2 + 2*m + 3*n + 4*n1 + 1 im = max(mineq,n1-meq) - if ( l_wn) then ! note: calling lsq when meq>n is corrupting the ! memory in some way, so just catch this here. mode = 2 + iter = 0 return end if @@ -235,11 +238,11 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & 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),jw,& + 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) + ldat,alphamin,alphamax,tolf,toldf,toldx) end subroutine slsqp !******************************************************************************* @@ -251,10 +254,10 @@ end subroutine slsqp ! l1 - line search, positive definite bfgs update subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& - r,l,x0,mu,s,u,v,w,iw,& + 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) + alphamin,alphamax,tolf,toldf,toldx) implicit none integer,intent(in) :: m @@ -286,7 +289,6 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& !! * `+(n1+mineq)*(n1-meq) + 2*meq + n1` for [[lsei]] !! !! with `mineq = m - meq + 2*n1` & `n1 = n+1` - integer,dimension(*),intent(inout) :: iw real(wp),intent(inout) :: t real(wp),intent(inout) :: f0 real(wp),intent(inout) :: h1 @@ -306,8 +308,13 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& 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) :: 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 :: i, j, k @@ -421,11 +428,7 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& 100 ireset = ireset + 1 if ( ireset>5 ) then ! check relaxed convergence in case of positive directional derivative - if ( (abs(f-f0) +! Check for convergence. + + function check_convergence(n,f,f0,x,x0,s,h3,acc,tolf,toldf,toldx,& + converged,not_converged) result(mode) + + implicit none + + integer,intent(in) :: n + real(wp),intent(in) :: f + real(wp),intent(in) :: f0 + real(wp),dimension(:),intent(in) :: x + real(wp),dimension(:),intent(in) :: x0 + real(wp),dimension(:),intent(in) :: s + real(wp),intent(in) :: h3 + real(wp),intent(in) :: acc + real(wp),intent(in) :: tolf + real(wp),intent(in) :: toldf + real(wp),intent(in) :: toldx + integer,intent(in) :: converged !! mode value if converged + integer,intent(in) :: not_converged !! mode value if not converged + integer :: mode + + logical :: ok ! temp variable + real(wp),dimension(n) :: xmx0 + + if (h3=zero) ok = abs(f)=zero) ok = abs(f-f0)=zero) then + xmx0 = x-x0 ! to avoid array temporary warning + ok = dnrm2(n,xmx0,1)zero) then @@ -1181,7 +1229,7 @@ end function diff !### History ! * Jacob Williams, refactored into modern Fortran, Jan. 2016. - subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,index,mode) + subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,mode) implicit none @@ -1200,14 +1248,6 @@ subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,index,mode) !! the dual solution vector. `w` will satisfy `w(i) = 0` !! for all `i` in set `p` and `w(i) <= 0` for all `i` in set `z`. real(wp),dimension(m),intent(inout) :: zz !! an m-array of working space. - integer,dimension(n),intent(out) :: index !! an integer working array. - !! on exit the contents of this array define the sets - !! `p` and `z` as follows: - !! - !! * `index(1:nsetp) = set p`. - !! * `index(iz1:iz2) = set z`. - !! - !! where: `iz1 = nsetp + 1 = npp1`, `iz2 = n` integer,intent(out) :: mode !! this is a success-failure flag with the following meanings: !! !! * ***1*** the solution has been computed successfully. @@ -1217,6 +1257,14 @@ subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,index,mode) integer :: i,ii,ip,iter,itmax,iz,iz1,iz2,izmax,j,jj,jz,l,npp1,nsetp,rtnkey real(wp) :: alpha,asave,cc,sm,ss,t,temp,unorm,up,wmax,ztest real(wp),dimension(1) :: dummy + integer,dimension(n) :: index !! an integer working array. + !! the contents of this array define the sets + !! `p` and `z` as follows: + !! + !! * `index(1:nsetp) = set p`. + !! * `index(iz1:iz2) = set z`. + !! + !! where: `iz1 = nsetp + 1 = npp1`, `iz2 = n` real(wp),parameter :: factor = 0.01_wp @@ -1228,12 +1276,9 @@ subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,index,mode) iter = 0 itmax = 3*n - ! initialize the arrays index() and x(). - - do i = 1 , n - x(i) = zero - index(i) = i - end do + ! initialize the arrays index(1:n) and x(1:n). + x = zero + index = [(i, i=1,n)] iz2 = n iz1 = 1 @@ -1377,12 +1422,7 @@ subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,index,mode) jj = index(ip) zz(ip) = zz(ip)/a(ip,jj) end do - ! if ( rtnkey==1 ) then !.....original - ! !continue - ! else if ( rtnkey/=2 ) then - ! return - ! end if - if (rtnkey/=1 .and. rtnkey/=2) return !......replaced with + if (rtnkey/=1 .and. rtnkey/=2) return ! ****** secondary loop begins here ****** @@ -1413,7 +1453,7 @@ subroutine nnls(a,mda,m,n,b,x,rnorm,w,zz,index,mode) ! if all new constrained coeffs are feasible then alpha will ! still = 2. if so exit from secondary loop to main loop. - if ( alpha==two ) then + if ( abs(alpha-two)<=zero ) then ! ****** end of secondary loop ****** do ip = 1 , nsetp @@ -1502,7 +1542,7 @@ end subroutine nnls !### History ! * Jacob Williams, refactored into modern Fortran, Jan. 2016. - subroutine hfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip) + subroutine hfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g) implicit none @@ -1524,8 +1564,6 @@ subroutine hfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip) !! defined by the `j-th` column vector of the array `b`. real(wp),dimension(n),intent(inout) :: h !! array of working space real(wp),dimension(n),intent(inout) :: g !! array of working space - integer,dimension(n),intent(inout) :: ip !! integer array of working space - !! recording permutation indices of column vectors real(wp),dimension(mdb,nb),intent(inout) :: b !! if `nb = 0` the subroutine will make no reference !! to the array `b`. if `nb > 0` the array `b` must !! initially contain the `m x nb` matrix `b` of the @@ -1535,6 +1573,8 @@ subroutine hfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip) integer :: i, ii, ip1, j, jb, jj, k, kp1, l, ldiag, lmax real(wp) :: hmax, sm, tmp logical :: need_lmax + integer,dimension(n) :: ip !! integer array of working space + !! recording permutation indices of column vectors real(wp),parameter :: factor = 0.001_wp @@ -1559,14 +1599,14 @@ subroutine hfti(a,mda,m,n,b,mdb,nb,tau,krank,rnorm,h,g,ip) h(l) = h(l) - a(j-1,l)**2 if ( h(l)>h(lmax) ) lmax = l end do - if ( diff(hmax+factor*h(lmax),hmax)>0 ) need_lmax = .false. + if ( diff(hmax+factor*h(lmax),hmax)>zero ) need_lmax = .false. end if if (need_lmax) then ! compute squared column lengths and find lmax lmax = j do l = j , n - h(l) = 0. + h(l) = zero do i = j , m h(l) = h(l) + a(i,l)**2 end do @@ -1742,17 +1782,17 @@ subroutine h12(mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv) do j = l1 , m cl = max(abs(u(1,j)),cl) end do - if ( cl<=0 ) return + if ( cl<=zero ) return clinv = one/cl sm = (u(1,lpivot)*clinv)**2 do j = l1 , m sm = sm + (u(1,j)*clinv)**2 end do cl = cl*sqrt(sm) - if ( u(1,lpivot)>0 ) cl = -cl + if ( u(1,lpivot)>zero ) cl = -cl up = u(1,lpivot) - cl u(1,lpivot) = cl - else if ( cl<=0 ) then + else if ( cl<=zero ) then return end if @@ -1760,7 +1800,7 @@ subroutine h12(mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv) ! apply the transformation i+u*(u**t)/b to c. b = up*u(1,lpivot) ! b must be nonpositive here. - if ( b<0 ) then + if ( bzero ) then sm = sm*b c(i2) = c(i2) + sm*up do i = l1 , m @@ -1831,7 +1871,7 @@ subroutine g1(a,b,c,s,sig) s = c*xr sig = abs(a)*yr else - if ( b/=zero ) then + if ( abs(b)>zero ) then xr = a/b yr = sqrt(one+xr**2) s = sign(one/yr,b) @@ -1883,7 +1923,7 @@ subroutine ldl(n,a,z,sigma,w) integer :: i , ij , j real(wp) :: t , v , u , tp , beta , alpha , delta , gamma - if ( sigma/=zero ) then + if ( abs(sigma)>zero ) then ij = 1 t = one/sigma if ( sigma<=zero ) then @@ -2004,12 +2044,12 @@ real(wp) function linmin(mode,ax,bx,f,tol,& if ( fu>fx ) then if ( u=x ) b = u - if ( fu<=fw .or. w==x ) then + if ( fu<=fw .or. abs(w-x)<=zero ) then v = w fv = fw w = u fw = fu - else if ( fu<=fv .or. v==x .or. v==w ) then + else if ( fu<=fv .or. abs(v-x)<=zero .or. abs(v-w)<=zero ) then v = u fv = fu end if diff --git a/src/slsqp_module.f90 b/src/slsqp_module.f90 index 8cd478e..faa8ae9 100644 --- a/src/slsqp_module.f90 +++ b/src/slsqp_module.f90 @@ -28,6 +28,10 @@ module slsqp_module 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 xf: if \( |x_{n+1} - x_n| < toldx \) then stop integer :: gradient_mode = 0 !! how the gradients are computed: !! @@ -50,9 +54,6 @@ module slsqp_module integer :: l_w = 0 !! size of `w` real(wp),dimension(:),allocatable :: w !! real work array - integer :: l_jw = 0 !! size of `jw` - integer,dimension(:),allocatable :: jw !! integer 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 @@ -137,7 +138,7 @@ end subroutine stop_iterations 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) + gradient_mode,gradient_delta,tolf,toldf,toldx) implicit none @@ -168,6 +169,9 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& !! 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 :: n1,mineq,i @@ -228,6 +232,10 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& end if + if (present(tolf)) me%tolf = tolf + if (present(toldf)) me%toldf = toldf + if (present(toldx)) me%toldx = toldx + status_ok = .true. me%n = n me%m = m @@ -253,9 +261,6 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& allocate(me%w(me%l_w)) me%w = zero - me%l_jw = mineq - allocate(me%jw(me%l_jw)) - me%jw = 0 if (present(gradient_mode)) then me%gradient_mode = gradient_mode if (present(gradient_delta)) then @@ -263,6 +268,7 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& end if end if end if + end subroutine initialize_slsqp !******************************************************************************* @@ -292,30 +298,41 @@ subroutine slsqp_wrapper(me,x,istat,iterations,status_message) !! **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) :: f !! objective function - real(wp),dimension(max(1,me%m)) :: c !! constraint vector - real(wp),dimension(max(1,me%m),me%n+1) :: a !! a matrix for [[slsqp]] - real(wp),dimension(me%n+1) :: g !! g matrix for [[slsqp]] - real(wp),dimension(me%m) :: cvec !! constraint vector - real(wp),dimension(me%n) :: dfdx !! objective function partials - real(wp),dimension(me%m,me%n) :: dcdx !! constraint partials - 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),dimension(me%n) :: delta !! perturbation step size 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),dimension(me%m) :: cvecr !! right function value to approximate constraints vector's gradient - real(wp),dimension(me%m) :: cvecl !! left function value to approximate constraints vector's gradient - real(wp) :: fact !! denominator factor for finite difference approximation + 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) @@ -379,7 +396,7 @@ subroutine slsqp_wrapper(me,x,istat,iterations,status_message) if (mode==0 .or. mode==1) then !function evaluation (f&c) call me%f(x,f,cvec) - c(1:me%m) = cvec + c(1:me%m) = cvec end if if (mode==0 .or. mode==-1) then !gradient evaluation (g&a) @@ -426,8 +443,9 @@ subroutine slsqp_wrapper(me,x,istat,iterations,status_message) !main routine: call slsqp(me%m,me%meq,la,me%n,x,me%xl,me%xu,& f,c,g,a,acc,iter,mode,& - me%w,me%l_w,me%jw,me%l_jw,& - me%slsqpb,me%linmin,me%alphamin,me%alphamax) + me%w,me%l_w,& + me%slsqpb,me%linmin,me%alphamin,me%alphamax,& + me%tolf,me%toldf,me%toldx) if (mode==1 .or. mode==-1) then !continue to next call diff --git a/src/slsqp_support.f90 b/src/slsqp_support.f90 index 2203d2c..d23fd7c 100644 --- a/src/slsqp_support.f90 +++ b/src/slsqp_support.f90 @@ -42,7 +42,7 @@ subroutine daxpy(n,da,dx,incx,dy,incy) integer :: i , incx , incy , ix , iy , m , mp1 , n if ( n<=0 ) return - if ( da==zero ) return + if ( abs(da)<=zero ) return if ( incx==1 .and. incy==1 ) then ! code for both increments equal to 1 @@ -243,7 +243,7 @@ real(wp) function dnrm2(n,x,incx) ! auxiliary routine: ! call dlassq( n, x, incx, scale, ssq ) do ix = 1 , 1 + (n-1)*incx , incx - if ( x(ix)/=zero ) then + if ( abs(x(ix))>zero ) then absxi = abs(x(ix)) if ( scale author: Jacob Williams +! +! Test for the [[slsqp_module]]. +! +! This is problem 71 from the Hock-Schittkowsky test suite: +! +!``` +! min x1*x4*(x1 + x2 + x3) + x3 +! s.t. x1*x2*x3*x4 - x5 - 25 = 0 +! x1**2 + x2**2 + x3**2 + x4**2 - 40 = 0 +! 1 <= x1,x2,x3,x4 <= 5 +! 0 <= x5 +! +! Starting point: +! x = (1, 5, 5, 1, -24) +! +! Optimal solution: +! x = (1.00000000, 4.74299963, 3.82114998, 1.37940829, 0) +!``` + + program slsqp_test_71 + + use slsqp_module + use slsqp_kinds + + implicit none + + integer,parameter :: n = 5 !! number of optimization variables + integer,parameter :: m = 2 !! total number of constraints + integer,parameter :: meq = 2 !! number of equality constraints + integer,parameter :: max_iter = 100 !! maximum number of allowed iterations + real(wp),dimension(n),parameter :: xl = [1.0_wp, 1.0_wp, 1.0_wp, 1.0_wp, 0.0_wp] !! lower bounds + real(wp),dimension(n),parameter :: xu = [5.0_wp, 5.0_wp, 5.0_wp, 5.0_wp, 1.0e10_wp] !! upper bounds + real(wp),parameter :: acc = 1.0e-8_wp !! tolerance + integer,parameter :: linesearch_mode = 1 !! use inexact linesearch. + + type(slsqp_solver) :: solver !! instantiate an slsqp solver + real(wp),dimension(n) :: x !! optimization variable vector + integer :: istat !! for solver status check + logical :: status_ok !! for initialization status check + integer :: iterations !! number of iterations by the solver + + x = [1.0_wp, 5.0_wp, 5.0_wp, 1.0_wp, -24.0_wp] !initial guess + + call solver%initialize(n,m,meq,max_iter,acc,func,grad,& + xl,xu,linesearch_mode=linesearch_mode,& + status_ok=status_ok,& + report=report_iteration) + + if (status_ok) then + call solver%optimize(x,istat,iterations) + write(*,*) '' + write(*,*) 'solution :', x + write(*,*) 'istat :', istat + write(*,*) 'iterations :', iterations + write(*,*) '' + else + error stop 'error calling slsqp.' + end if + + contains + + subroutine func(me,x,f,c) + + implicit none + class(slsqp_solver),intent(inout) :: me + real(wp),dimension(:),intent(in) :: x !! optimization variable vector + real(wp),intent(out) :: f !! value of the objective function + real(wp),dimension(:),intent(out) :: c !! the constraint vector `dimension(m)`, + !! equality constraints (if any) first. + + f = x(1)*x(4)*(x(1)+x(2)+x(3)) + x(3) + + c(1) = x(1)*x(2)*x(3)*x(4) - x(5) - 25.0_wp + c(2) = x(1)**2 + x(2)**2 + x(3)**2 + x(4)**2 - 40.0_wp + + end subroutine func + + subroutine grad(me,x,g,a) + !! gradients for [[func]]. + implicit none + class(slsqp_solver),intent(inout) :: me + real(wp),dimension(:),intent(in) :: x !! optimization variable vector + real(wp),dimension(:),intent(out) :: g !! objective function partials w.r.t x `dimension(n)` + real(wp),dimension(:,:),intent(out) :: a !! gradient matrix of constraints w.r.t. x `dimension(m,n)` + + g(1) = x(4)*(2.0_wp*x(1)+x(2)+x(3)) + g(2) = x(1)*x(4) + g(3) = x(1)*x(4) + 1.0_wp + g(4) = x(1)*(x(1)+x(2)+x(3)) + g(5) = 0.0_wp + + a = 0.0_wp + a(1,1) = x(2)*x(3)*x(4) + a(1,2) = x(1)*x(3)*x(4) + a(1,3) = x(1)*x(2)*x(4) + a(1,4) = x(1)*x(2)*x(3) + a(1,5) = -1.0_wp + a(2,1) = 2.0_wp*x(1) + a(2,2) = 2.0_wp*x(2) + a(2,3) = 2.0_wp*x(3) + a(2,4) = 2.0_wp*x(4) + + end subroutine grad + + subroutine report_iteration(me,iter,x,f,c) + use, intrinsic :: iso_fortran_env, only: output_unit + !! report an iteration (print to the console). + implicit none + class(slsqp_solver),intent(inout) :: me + integer,intent(in) :: iter + real(wp),dimension(:),intent(in) :: x + real(wp),intent(in) :: f + real(wp),dimension(:),intent(in) :: c + + !write a header: + if (iter==0) then + write(output_unit,'(*(A20,1X))') 'iteration', 'x(1)', 'x(2)', 'x(3)', 'x(4)', 'x(5)', 'f', 'c(1)', 'c(2)' + end if + + !write the iteration data: + write(output_unit,'(I20,1X,(*(F20.16,1X)))') iter,x,f,c + + end subroutine report_iteration + + end program slsqp_test_71 +!******************************************************************************* \ No newline at end of file diff --git a/src/tests/slsqp_test_stopping_criterion.f90 b/src/tests/slsqp_test_stopping_criterion.f90 new file mode 100644 index 0000000..822a539 --- /dev/null +++ b/src/tests/slsqp_test_stopping_criterion.f90 @@ -0,0 +1,143 @@ +!******************************************************************************* +!> author: Pierre Payen +! +! Test for the [[slsqp_module]]. + + module report_variable + + use slsqp_kinds + + implicit none + + real(wp) :: flast + real(wp),dimension(2) :: xlast + real(wp),dimension(1) :: clast + + end module + + program slsqp_test_criterion + + use slsqp_module + use slsqp_kinds + use report_variable + + implicit none + + integer,parameter :: n = 2 !! number of optimization variables + integer,parameter :: m = 1 !! total number of constraints + integer,parameter :: meq = 1 !! number of equality constraints + integer,parameter :: max_iter = 100 !! maximum number of allowed iterations + real(wp),dimension(n),parameter :: xl = [-1.0_wp, -1.0_wp] !! lower bounds + real(wp),dimension(n),parameter :: xu = [ 1.0_wp, 1.0_wp] !! upper bounds + real(wp),parameter :: acc = 1.e-8_wp !! accuracy tolerance + real(wp),parameter :: tolf = 0.e+0_wp !! accuracy tolerance over f: if |f| < tolf then stop + real(wp),parameter :: toldf = 0.e+0_wp !! accuracy tolerance over df: if |fn+1 - fn| < toldf then stop. + !! It's different from acc in the event of positive derivative + real(wp),parameter :: toldx = 0.e+0_wp !! accuracy tolerance over dx: if |xn+1 - xn| < toldx then stop + integer,parameter :: linesearch_mode = 1 !! use inexact linesearch. + + type(slsqp_solver) :: solver !! instantiate an slsqp solver + real(wp),dimension(n) :: x !! optimization variable vector + integer :: istat !! for solver status check + logical :: status_ok !! for initialization status check + integer :: iterations !! number of iterations by the solver + + x = [0.1_wp, 0.1_wp] !initial guess + xlast = x + call rosenbrock_func(solver,xlast,flast,clast) + + write(*,*) 'The stopping criterion are' + write(*,*) 'acc ',acc + write(*,*) 'tolf ',tolf + write(*,*) 'toldf',toldf + write(*,*) 'toldx',toldx + + call solver%initialize(n,m,meq,max_iter,acc,rosenbrock_func,rosenbrock_grad,& + xl,xu,linesearch_mode=linesearch_mode,status_ok=status_ok,& + report=report_iteration,tolf=tolf,toldf=toldf,toldx=toldx) + !alphamin=0.1_wp, alphamax=0.5_wp) !to limit search steps + + if (status_ok) then + call solver%optimize(x,istat,iterations) + write(*,*) '' + write(*,*) 'solution :', x + write(*,*) 'istat :', istat + write(*,*) 'iterations :', iterations + write(*,*) '' + else + error stop 'error calling slsqp.' + end if + + !solution: x1 = 0.78641515097183889 + ! x2 = 0.61769831659541152 + ! f = 4.5674808719160388E-002 + ! c = 2.8654301154062978E-012 + + contains + + subroutine rosenbrock_func(me,x,f,c) + !! Rosenbrock function + !! + !! Minimize the Rosenbrock function: \( f(x) = 100 (x_2 - x_1)^2 + (1 - x_1)^2 \), + !! subject to the inequality constraint: \( x_1^2 + x_2^2 \le 1 \). + !! + !! see: http://www.mathworks.com/help/optim/ug/example-nonlinear-constrained-minimization.html + implicit none + class(slsqp_solver),intent(inout) :: me + real(wp),dimension(:),intent(in) :: x !! optimization variable vector + real(wp),intent(out) :: f !! value of the objective function + real(wp),dimension(:),intent(out) :: c !! the constraint vector `dimension(m)`, + !! equality constraints (if any) first. + + f = 100.0_wp*(x(2) - x(1)**2)**2 + (1.0_wp - x(1))**2 !objective function + c(1) = 1.0_wp - x(1)**2 - x(2)**2 !equality constraint (>=0) + + end subroutine rosenbrock_func + + subroutine rosenbrock_grad(me,x,g,a) + !! gradients for [[rosenbrock_func]]. + implicit none + class(slsqp_solver),intent(inout) :: me + real(wp),dimension(:),intent(in) :: x !! optimization variable vector + real(wp),dimension(:),intent(out) :: g !! objective function partials w.r.t x `dimension(n)` + real(wp),dimension(:,:),intent(out) :: a !! gradient matrix of constraints w.r.t. x `dimension(m,n)` + + g(1) = -400.0_wp*(x(2)-x(1)**2)*x(1) - 2.0_wp*(1.0_wp-x(1)) !df/x1 + g(2) = 200.0_wp*(x(2) - x(1)**2) !df/x2 + + a(1,1) = -2.0_wp*x(1) ! dc/dx1 + a(1,2) = -2.0_wp*x(2) ! dc/dx2 + + end subroutine rosenbrock_grad + + subroutine report_iteration(me,iter,x,f,c) + use, intrinsic :: iso_fortran_env, only: output_unit + use report_variable + !! report an iteration (print to the console). + implicit none + class(slsqp_solver),intent(inout) :: me + integer,intent(in) :: iter + real(wp),dimension(:),intent(in) :: x + real(wp),intent(in) :: f + real(wp),dimension(:),intent(in) :: c + + !write a header: + if (iter==0) then + write(output_unit,'(*(A17,1X))') & + 'iteration', 'f','|fn+1-fn|', '||xn+1-xn||', 'c(1)',& + 'f