From ffd67ed66ddcc057e6a4f952c5cc0ae8b796b8b5 Mon Sep 17 00:00:00 2001 From: Pierre Payen Date: Mon, 16 Sep 2019 16:33:55 +0200 Subject: [PATCH 1/3] removed gfortran warning when comparing real number with equalities or inequalities --- src/slsqp_core.f90 | 12 ++++++------ src/slsqp_support.f90 | 4 ++-- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/slsqp_core.f90 b/src/slsqp_core.f90 index 5048afa..ccbe835 100644 --- a/src/slsqp_core.f90 +++ b/src/slsqp_core.f90 @@ -1413,7 +1413,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)<=0. ) then ! ****** end of secondary loop ****** do ip = 1 , nsetp @@ -1773,7 +1773,7 @@ subroutine h12(mode,lpivot,l1,m,u,iue,up,c,ice,icv,ncv) sm = sm + c(i3)*u(1,i) i3 = i3 + ice end do - if ( sm/=0 ) then + if ( abs(sm)>0 ) then sm = sm*b c(i2) = c(i2) + sm*up do i = l1 , m @@ -1831,7 +1831,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 +1883,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 +2004,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)<=0._wp ) 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)<=0._wp .or. abs(v-w)<=0._wp ) then v = u fv = fu end if 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 Date: Mon, 16 Sep 2019 16:47:25 +0200 Subject: [PATCH 2/3] Added test case to test different stopping criterion --- src/slsqp_module.f90 | 5 +- src/tests/slsqp_test_stopping_criterion.f90 | 132 ++++++++++++++++++++ 2 files changed, 136 insertions(+), 1 deletion(-) create mode 100644 src/tests/slsqp_test_stopping_criterion.f90 diff --git a/src/slsqp_module.f90 b/src/slsqp_module.f90 index 8cd478e..d81f725 100644 --- a/src/slsqp_module.f90 +++ b/src/slsqp_module.f90 @@ -27,7 +27,10 @@ module slsqp_module integer :: m = 0 !! number of constraints (\( m \ge 0 \)) 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) :: acc = zero !! accuracy tolerance + real(wp) :: tolf = zero !! accuracy tolerance over f: if |f| < tolf then stop + real(wp) :: toldf = zero !! accuracy tolerance over df: if |fn+1 - fn| < toldf then stop + real(wp) :: toldx = zero !! accuracy tolerance over xf: if |xn+1 - xn| < toldx then stop integer :: gradient_mode = 0 !! how the gradients are computed: !! diff --git a/src/tests/slsqp_test_stopping_criterion.f90 b/src/tests/slsqp_test_stopping_criterion.f90 new file mode 100644 index 0000000..b46cfcc --- /dev/null +++ b/src/tests/slsqp_test_stopping_criterion.f90 @@ -0,0 +1,132 @@ +!******************************************************************************* +!> 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 = 0 !! 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.0e-8_wp !! accuracy tolerance + real(wp),parameter :: tolf = 1.0e-8_wp !! accuracy tolerance over f: if |f| < tolf then stop + real(wp),parameter :: toldf = 1.0e-8_wp !! accuracy tolerance over df: if |fn+1 - fn| < toldf then stop + real(wp),parameter :: toldx = 1.0e-8_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) + + 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) + !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,'(*(A20,1X))') 'iteration', 'f','|df|','x(1)', 'x(2)', '|dx|', 'c(1)' + end if + + !write the iteration data: + write(output_unit,'(I20,1X,(*(F20.16,1X)))') iter,f,abs(f-flast),x,sqrt(sum((x-xlast)**2)),c + + xlast = x + flast = f + clast = c + + end subroutine report_iteration + + end program +!******************************************************************************* From 729684b3399b0c01111fb9469f90ff2aad78e7b9 Mon Sep 17 00:00:00 2001 From: pirpyn Date: Mon, 16 Sep 2019 23:20:18 +0200 Subject: [PATCH 3/3] added 3 stopping criterion in the slsqp core and module file --- src/slsqp_core.f90 | 21 ++++++++++------ src/slsqp_module.f90 | 16 +++++++++--- src/tests/slsqp_test_stopping_criterion.f90 | 27 +++++++++++++++------ 3 files changed, 44 insertions(+), 20 deletions(-) diff --git a/src/slsqp_core.f90 b/src/slsqp_core.f90 index ccbe835..6a941df 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) + jw,l_jw,sdat,ldat,alphamin,alphamax,tolf,toldf,toldx) implicit none @@ -198,7 +198,9 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & 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 |fn+1 - fn| < toldf then stop + real(wp),intent(in) :: toldx !! stopping criterion if ||xn+1 - xn|| < toldx then stop integer :: il , im , ir , is , iu , iv , iw , ix , mineq, n1 ! check length of working arrays @@ -239,7 +241,7 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & 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 !******************************************************************************* @@ -254,7 +256,7 @@ 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,& 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 @@ -308,7 +310,9 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& 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 |fn+1 - fn| < toldf then stop + real(wp),intent(in) :: toldx !! stopping criterion if ||xn+1 - xn|| < toldx then stop integer :: i, j, k if ( mode<0 ) then @@ -421,7 +425,8 @@ 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)