From 58aa1a8a95468ee4a77b49355bfd90bf201c7a56 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 4 Jun 2023 10:34:04 -0500 Subject: [PATCH 1/5] port some updates from SciPy check constraint violation when exiting due to too many resets check for NaNs do not exit with convergence on steps where SQP is inconsistent Fixes #35 --- src/slsqp_core.f90 | 57 +++++++++++++++++++++++++++++++++------------- 1 file changed, 41 insertions(+), 16 deletions(-) diff --git a/src/slsqp_core.f90 b/src/slsqp_core.f90 index 373300b..977bbc0 100644 --- a/src/slsqp_core.f90 +++ b/src/slsqp_core.f90 @@ -8,7 +8,8 @@ module slsqp_core use slsqp_kinds use slsqp_support - use bvls_module, only: bvls_wrapper + use bvls_module, only: bvls_wrapper + use, intrinsic :: ieee_arithmetic, only: ieee_is_nan, ieee_value, ieee_quiet_nan implicit none @@ -325,6 +326,9 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& integer,intent(in) :: nnls_mode !! which NNLS method to use integer :: i, j, k + logical :: inconsistent_linearization !! if the SQP problem is inconsistent + + inconsistent_linearization = .false. ! initialize if ( mode<0 ) then @@ -378,9 +382,15 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& call dscal(n,h4,u,1) call daxpy(n,one-h4,v,1,u,1) end if - call ldl(n,l,u,+one/h1,v) - call ldl(n,l,v,-one/h2,u) - + if (h1==zero .or. h2==zero) then + ! Singular update: reset hessian. + ! [ JW : this is based on a SciPy update ] + call reset_bfgs_matrix() + if ( ireset>5 ) return + else + call ldl(n,l,u,+one/h1,v) + call ldl(n,l,v,-one/h2,u) + end if ! end of main iteration else if ( mode==0 ) then @@ -423,14 +433,14 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& select case (iexact) case(0) if ( h1<=h3/ten .or. line>10 ) then - call convergence_check() + call convergence_check(acc,0,-1) else alpha = min(max(h3/(two*(h3-h1)),alphamin),alphamax) call inexact_linesearch() end if case(1) call exact_linesearch() - if ( line==3 ) call convergence_check() + if ( line==3 ) call convergence_check(acc,0,-1) end select return @@ -455,10 +465,13 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& ! augmented problem for inconsistent linearization + inconsistent_linearization = .false. ! initialize if ( mode==6 ) then if ( n==meq ) mode = 4 end if if ( mode==4 ) then + ! Will reject this iteration if the SQP problem is inconsistent, + inconsistent_linearization = .true. do j = 1 , m if ( j<=meq ) then a(j,n1) = -c(j) @@ -519,7 +532,9 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& ! check convergence mode = 0 - if ( h15 ) then ! check relaxed convergence in case of positive directional derivative - mode = check_convergence(n,f,f0,x,x0,s,h3,tol,tolf,toldf,toldx,0,8) - return + ! [ JW: reuse this routine so that h3 is recomputed. + ! this is based on a SciPy update to SLSQP ] + call convergence_check(tol,0,8) + ! the caller should return in this case else l(1) = zero call dcopy(n2,l(1),0,l,1) @@ -595,17 +613,22 @@ subroutine exact_linesearch() ! 400 end if end subroutine exact_linesearch - subroutine convergence_check() ! 500 + subroutine convergence_check(tolerance,converged,not_converged) ! 500 + real(wp),intent(in) :: tolerance !! tolerance + integer,intent(in) :: converged !! mode value if converged + integer,intent(in) :: not_converged !! mode value if not converged + h3 = zero do j = 1 , m - if ( j<=meq ) then + if (j<=meq) then h1 = c(j) else h1 = zero end if h3 = h3 + max(-c(j),h1) end do - mode = check_convergence(n,f,f0,x,x0,s,h3,acc,tolf,toldf,toldx,0,-1) + mode = check_convergence(n,f,f0,x,x0,s,h3,tolerance,tolf,toldf,toldx,& + converged,not_converged,inconsistent_linearization) end subroutine convergence_check end subroutine slsqpb @@ -616,7 +639,7 @@ end subroutine slsqpb ! Check for convergence. pure function check_convergence(n,f,f0,x,x0,s,h3,acc,tolf,toldf,toldx,& - converged,not_converged) result(mode) + converged,not_converged,inconsistent_linearization) result(mode) implicit none @@ -633,12 +656,13 @@ pure function check_convergence(n,f,f0,x,x0,s,h3,acc,tolf,toldf,toldx,& real(wp),intent(in) :: toldx integer,intent(in) :: converged !! mode value if converged integer,intent(in) :: not_converged !! mode value if not converged + logical,intent(in) :: inconsistent_linearization !! if the SQP problem is inconsistent (will return `not_converged`) integer :: mode logical :: ok ! temp variable real(wp),dimension(n) :: xmx0 - if (h3>=acc) then + if (h3>=acc .or. inconsistent_linearization .or. ieee_is_nan(f)) then mode = not_converged else @@ -1097,7 +1121,7 @@ subroutine lsi(e,f,g,h,le,me,lg,mg,n,x,xnorm,w,mode,max_iter_ls,nnls_mode) mode = 5 do i = 1 , mg do j = 1 , n - if ( abs(e(j,j))zero) then ! compute solution of primal problem fac=one-ddot(m,h,1,w(iy),1) + if (ieee_is_nan(fac)) return if (diff(one+fac,one)>zero) then mode=1 fac=one/fac From 58ce1eb414eee61eede5823e811678e87cd64a54 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 4 Jun 2023 13:54:44 -0500 Subject: [PATCH 2/5] better handling of unbounded variables Can now specify an infinite_bound value to represent a bound that doesn't exist. Also, NaN can be used. Based on SciPy updates here: * https://github.com/scipy/scipy/commit/db343b85053f365dc9a747923328d4e062375d34 * https://github.com/scipy/scipy/commit/bac03893d698657f4d0db7c93e88f4ee07e4438d Fixes #37 --- src/slsqp_core.f90 | 119 ++++++++++++++++--------- src/slsqp_module.f90 | 33 +++++-- test/slsqp_test_stopping_criterion.f90 | 6 ++ 3 files changed, 108 insertions(+), 50 deletions(-) diff --git a/src/slsqp_core.f90 b/src/slsqp_core.f90 index 977bbc0..68575e5 100644 --- a/src/slsqp_core.f90 +++ b/src/slsqp_core.f90 @@ -113,7 +113,7 @@ module slsqp_core subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & sdat,ldat,alphamin,alphamax,tolf,toldf,toldx,& - max_iter_ls,nnls_mode) + max_iter_ls,nnls_mode,infinite_bound) implicit none @@ -204,10 +204,21 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & real(wp),intent(in) :: toldx !! stopping criterion if \( ||x_{n+1} - x_n|| < toldx \) then stop. integer,intent(in) :: max_iter_ls !! maximum number of iterations in the [[nnls]] problem integer,intent(in) :: nnls_mode !! which NNLS method to use + real(wp),intent(in) :: infinite_bound !! "infinity" for the upper and lower bounds. + !! if `xl<=-infinite_bound` or `xu>=infinite_bound` + !! then these bounds are considered nonexistant. + !! If `infinite_bound=0` then `huge()` is used for this. integer :: il , im , ir , is , iu , iv , iw , ix , mineq, n1 + real(wp) :: infBnd !! local copy of `infinite_bound` - ! check length of working arrays + if (infinite_bound==zero) then + infBnd = huge(one) + else + infBnd = abs(infinite_bound) + end if + + ! check length of working arrays n1 = n + 1 mineq = m - meq + n1 + n1 il = (3*n1+m)*(n1+1) + (n1-meq+1)*(mineq+2) + 2*mineq + (n1+mineq)& @@ -227,7 +238,7 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & return end if - ! prepare data for calling sqpbdy - initial addresses in w + ! prepare data for calling sqpbdy - initial addresses in w im = 1 il = im + max(1,m) @@ -248,7 +259,7 @@ subroutine slsqp(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,w,l_w, & 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,tolf,toldf,toldx,& - max_iter_ls,nnls_mode) + max_iter_ls,nnls_mode,infBnd) end subroutine slsqp !******************************************************************************* @@ -264,7 +275,7 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& t,f0,h1,h2,h3,h4,n1,n2,n3,t0,gs,tol,line,& alpha,iexact,incons,ireset,itermx,ldat,& alphamin,alphamax,tolf,toldf,toldx,& - max_iter_ls,nnls_mode) + max_iter_ls,nnls_mode,infBnd) implicit none integer,intent(in) :: m @@ -324,6 +335,7 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& real(wp),intent(in) :: toldx !! stopping criterion if \( ||x_{n+1} - x_n|| < toldx \) then stop integer,intent(in) :: max_iter_ls !! maximum number of iterations in the [[nnls]] problem integer,intent(in) :: nnls_mode !! which NNLS method to use + real(wp),intent(in) :: infBnd !! "infinity" for the upper and lower bounds. integer :: i, j, k logical :: inconsistent_linearization !! if the SQP problem is inconsistent @@ -461,7 +473,7 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& call daxpy(n,-one,x,1,u,1) call daxpy(n,-one,x,1,v,1) h4 = one - call lsq(m,meq,n,n3,la,l,g,a,c,u,v,s,r,w,mode,max_iter_ls,nnls_mode) + call lsq(m,meq,n,n3,la,l,g,a,c,u,v,s,r,w,mode,max_iter_ls,nnls_mode,infBnd) ! augmented problem for inconsistent linearization @@ -489,7 +501,7 @@ subroutine slsqpb(m,meq,la,n,x,xl,xu,f,c,g,a,acc,iter,mode,& v(n1) = one incons = 0 do - call lsq(m,meq,n1,n3,la,l,g,a,c,u,v,s,r,w,mode,max_iter_ls,nnls_mode) + call lsq(m,meq,n1,n3,la,l,g,a,c,u,v,s,r,w,mode,max_iter_ls,nnls_mode,infBnd) h4 = one - s(n1) if ( mode==4 ) then l(n3) = ten*l(n3) @@ -595,7 +607,7 @@ subroutine inexact_linesearch() ! 300 call dscal(n,alpha,s,1) call dcopy(n,x0,1,x,1) call daxpy(n,one,s,1,x,1) - call enforce_bounds(x,xl,xu) ! ensure that x doesn't violate bounds + call enforce_bounds(x,xl,xu,infBnd) ! ensure that x doesn't violate bounds mode = 1 end subroutine inexact_linesearch @@ -721,7 +733,7 @@ end function check_convergence ! * coded dieter kraft, april 1987 ! * revised march 1989 - subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode) + subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode,infBnd) implicit none @@ -744,6 +756,7 @@ subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode) !! * **7:** rank defect in [[hfti]] integer,intent(in) :: max_iter_ls !! maximum number of iterations in the [[nnls]] problem integer,intent(in) :: nnls_mode !! which NNLS method to use + real(wp),intent(in) :: infbnd !! "infinity" for the upper and lower bounds. real(wp),dimension(nl) :: l real(wp),dimension(n) :: g @@ -755,7 +768,7 @@ subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode) real(wp) :: diag , xnorm integer :: i , ic , id , ie , if , ig , ih , il , im , ip , & iu , iw , i1 , i2 , i3 , i4 , mineq , & - m1 , n1 , n2 , n3 + m1 , n1 , n2 , n3, num_unbounded, j n1 = n + 1 mineq = m - meq @@ -822,57 +835,75 @@ subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode) ig = id + meq if ( mineq>0 ) then - ! recover matrix g from lower part of a + ! recover matrix g from lower part of a + ! The matrix G(mineq+2*n,m1) is stored at w(ig) + ! Not all rows will be filled if some of the upper/lower + ! bounds are unbounded. do i = 1 , mineq call dcopy(n,a(meq+i,1),la,w(ig-1+i),m1) end do end if - ! augment matrix g by +i and -i - - ip = ig + mineq - do i = 1 , n - w(ip-1+i) = zero - call dcopy(n,w(ip-1+i),0,w(ip-1+i),m1) - end do - w(ip) = one - call dcopy(n,w(ip),0,w(ip),m1+1) - - im = ip + n - do i = 1 , n - w(im-1+i) = zero - call dcopy(n,w(im-1+i),0,w(im-1+i),m1) - end do - w(im) = -one - call dcopy(n,w(im),0,w(im),m1+1) - ih = ig + m1*n + iw = ih + mineq + 2*n if ( mineq>0 ) then ! recover h from lower part of b + ! The vector H(mineq+2*n) is stored at w(ih) call dcopy(mineq,b(meq+1),1,w(ih),1) call dscal(mineq,-one,w(ih),1) end if + ! augment matrix g by +i and -i, and, ! augment vector h by xl and xu + ! NaN or infBnd value indicates no bound + ip = ig + mineq il = ih + mineq - call dcopy(n,xl,1,w(il),1) - iu = il + n - call dcopy(n,xu,1,w(iu),1) - call dscal(n,-one,w(iu),1) + num_unbounded = 0 + + do i=1,n + if (ieee_is_nan(xl(i)) .or. xl(i)<=-infbnd) then + num_unbounded = num_unbounded + 1 + else + w(il) = xl(i) + do j=1,n + w(ip + m1*(j-1)) = zero + end do + w(ip + m1*(i-1)) = one + ip = ip + 1 + il = il + 1 + end if + end do - iw = iu + n + do i=1,n + if (ieee_is_nan(xu(i)) .or. xu(i)>=infbnd) then + num_unbounded = num_unbounded + 1 + else + w(il) = -xu(i) + do j=1,n + w(ip + m1*(j-1)) = zero + end do + w(ip + m1*(i-1)) = -one + ip = ip + 1 + il = il + 1 + end if + end do call lsei(w(ic),w(id),w(ie),w(if),w(ig),w(ih),max(1,meq),meq,n,n, & - m1,m1,n,x,xnorm,w(iw),mode,max_iter_ls,nnls_mode) + m1,m1-num_unbounded,n,x,xnorm,w(iw),mode,max_iter_ls,nnls_mode) if ( mode==1 ) then - ! restore lagrange multipliers + ! restore lagrange multipliers (only for user-defined variables) call dcopy(m,w(iw),1,y(1),1) - call dcopy(n3,w(iw+m),1,y(m+1),1) - call dcopy(n3,w(iw+m+n),1,y(m+n3+1),1) - call enforce_bounds(x,xl,xu) ! to ensure that bounds are not violated + if (n3 > 0) then + !set rest of the multipliers to nan (they are not used) + y(m+1) = ieee_value(one, ieee_quiet_nan) + do i=m+2,m+n3+n3 + y(i) = y(m+1) + end do + end if + call enforce_bounds(x,xl,xu,infbnd) ! to ensure that bounds are not violated end if end subroutine lsq @@ -2236,19 +2267,21 @@ end function linmin !******************************************************************************* !> -! enforce the bound constraints on x. +! enforce the bound constraints on `x`. - pure subroutine enforce_bounds(x,xl,xu) + subroutine enforce_bounds(x,xl,xu,infbnd) implicit none real(wp),dimension(:),intent(inout) :: x !! optimization variable vector real(wp),dimension(:),intent(in) :: xl !! lower bounds (must be same dimension as `x`) real(wp),dimension(:),intent(in) :: xu !! upper bounds (must be same dimension as `x`) + real(wp),intent(in) :: infbnd !! "infinity" for the upper and lower bounds. + !! Note that `NaN` may also be used to indicate no bound. - where (x-infbnd .and. .not. ieee_is_nan(xl)) x = xl - elsewhere (x>xu) + elsewhere (x>xu .and. xu=infinite_bound` + !! then these bounds are considered nonexistant. + contains private @@ -149,7 +154,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,tolf,toldf,toldx,& - max_iter_ls,nnls_mode) + max_iter_ls,nnls_mode,infinite_bound) implicit none @@ -161,8 +166,10 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& procedure(func) :: f !! problem function procedure(grad) :: g !! function to compute gradients (must be !! associated if `gradient_mode=0`) - real(wp),dimension(n),intent(in) :: xl !! lower bounds on `x` - real(wp),dimension(n),intent(in) :: xu !! upper bounds on `x` + real(wp),dimension(n),intent(in) :: xl !! lower bounds on `x`. + !! `xl(i)=NaN` (or `xl(i)<=-infinite_bound`) indicates to ignore `i`th bound + real(wp),dimension(n),intent(in) :: xu !! upper bounds on `x`. + !! `xu(i)=NaN` (or `xu(i)>=infinite_bound`) indicates to ignore `i`th bound real(wp),intent(in) :: acc !! accuracy logical,intent(out) :: status_ok !! will be false if there were errors integer,intent(in),optional :: linesearch_mode !! 1 = inexact (default), 2 = exact @@ -188,6 +195,10 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& !! !! 1. Use the original [[nnls]] !! 2. Use the newer [[bvls]] + real(wp),intent(in),optional :: infinite_bound !! "infinity" for the upper and lower bounds. + !! if `xl<=-infinite_bound` or `xu>=infinite_bound` + !! then these bounds are considered nonexistant. + !! If not present then `huge()` is used for this. integer :: n1,mineq,i @@ -207,10 +218,10 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& call me%report_message('error: invalid m value:', ival=m) else if (n<1) then call me%report_message('error: invalid n value:', ival=n) - else if (any(xl>xu)) then + else if (any(xl>xu .and. .not. ieee_is_nan(xl) .and. .not. ieee_is_nan(xu))) then call me%report_message('error: lower bounds must be <= upper bounds.') do i=1,n - if (xl(i)>xu(i)) then + if (xl(i)>xu(i) .and. .not. ieee_is_nan(xl(i)) .and. .not. ieee_is_nan(xu(i))) then call me%report_message(' xl(i)>xu(i) for variable',ival=i) end if end do @@ -294,6 +305,13 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& me%gradient_delta = gradient_delta end if end if + + if (present(infinite_bound)) then + me%infinite_bound = abs(infinite_bound) + else + me%infinite_bound = huge(one) + end if + end if end subroutine initialize_slsqp @@ -475,7 +493,8 @@ subroutine slsqp_wrapper(me,x,istat,iterations,status_message) me%w,me%l_w, & me%slsqpb,me%linmin,me%alphamin,me%alphamax,& me%tolf,me%toldf,me%toldx,& - me%max_iter_ls,me%nnls_mode) + me%max_iter_ls,me%nnls_mode,& + me%infinite_bound) if (mode==1 .or. mode==-1) then !continue to next call diff --git a/test/slsqp_test_stopping_criterion.f90 b/test/slsqp_test_stopping_criterion.f90 index dfe9a75..ba90f38 100644 --- a/test/slsqp_test_stopping_criterion.f90 +++ b/test/slsqp_test_stopping_criterion.f90 @@ -32,6 +32,12 @@ program slsqp_test_stopping_criterion real(wp),dimension(2) :: xlast real(wp),dimension(1) :: clast + write(*,*) '' + write(*,*) '-------------------------------' + write(*,*) ' slsqp_test_stopping_criterion ' + write(*,*) '-------------------------------' + write(*,*) '' + x = [0.1_wp, 0.1_wp] !initial guess xlast = x call rosenbrock_func(solver,xlast,flast,clast) From a905f5913f8f660c6918c7a63dc4966a265c135d Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 4 Jun 2023 16:22:41 -0500 Subject: [PATCH 3/5] added test and some code consolidation --- src/slsqp_core.f90 | 30 ++++++++++++++------------ test/slsqp_test_stopping_criterion.f90 | 12 +++++++++-- 2 files changed, 26 insertions(+), 16 deletions(-) diff --git a/src/slsqp_core.f90 b/src/slsqp_core.f90 index 68575e5..da466f1 100644 --- a/src/slsqp_core.f90 +++ b/src/slsqp_core.f90 @@ -866,13 +866,7 @@ subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode,infB if (ieee_is_nan(xl(i)) .or. xl(i)<=-infbnd) then num_unbounded = num_unbounded + 1 else - w(il) = xl(i) - do j=1,n - w(ip + m1*(j-1)) = zero - end do - w(ip + m1*(i-1)) = one - ip = ip + 1 - il = il + 1 + call update_w(xl(i), one) end if end do @@ -880,13 +874,7 @@ subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode,infB if (ieee_is_nan(xu(i)) .or. xu(i)>=infbnd) then num_unbounded = num_unbounded + 1 else - w(il) = -xu(i) - do j=1,n - w(ip + m1*(j-1)) = zero - end do - w(ip + m1*(i-1)) = -one - ip = ip + 1 - il = il + 1 + call update_w(xu(i), -one) end if end do @@ -906,6 +894,20 @@ subroutine lsq(m,meq,n,nl,la,l,g,a,b,xl,xu,x,y,w,mode,max_iter_ls,nnls_mode,infB call enforce_bounds(x,xl,xu,infbnd) ! to ensure that bounds are not violated end if + contains + + subroutine update_w(val, fact) + real(wp),intent(in) :: val !! xu(i) or xl(i) + real(wp),intent(in) :: fact !! -1 or 1 + w(il) = fact*val + do j=1,n + w(ip + m1*(j-1)) = zero + end do + w(ip + m1*(i-1)) = fact + ip = ip + 1 + il = il + 1 + end subroutine update_w + end subroutine lsq !******************************************************************************* diff --git a/test/slsqp_test_stopping_criterion.f90 b/test/slsqp_test_stopping_criterion.f90 index ba90f38..1d5d41e 100644 --- a/test/slsqp_test_stopping_criterion.f90 +++ b/test/slsqp_test_stopping_criterion.f90 @@ -7,6 +7,7 @@ program slsqp_test_stopping_criterion use slsqp_module use slsqp_kinds + use, intrinsic :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan implicit none @@ -14,8 +15,6 @@ program slsqp_test_stopping_criterion 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. @@ -31,6 +30,15 @@ program slsqp_test_stopping_criterion real(wp) :: flast real(wp),dimension(2) :: xlast real(wp),dimension(1) :: clast + real(wp),dimension(n) :: xl !! lower bounds + real(wp),dimension(n) :: xu !! upper bounds + real(wp) :: nan !! not a number + + nan = ieee_value(1.0_wp, ieee_quiet_nan) + + ! test with some missing bounds: + xl = [-1.0_wp, nan] + xu = [nan, 1.0_wp] write(*,*) '' write(*,*) '-------------------------------' From 07771f1b61acfb290adb7adeb242c96e3f81c424 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 4 Jun 2023 16:59:32 -0500 Subject: [PATCH 4/5] test update --- test/slsqp_test_stopping_criterion.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/slsqp_test_stopping_criterion.f90 b/test/slsqp_test_stopping_criterion.f90 index 1d5d41e..a483c85 100644 --- a/test/slsqp_test_stopping_criterion.f90 +++ b/test/slsqp_test_stopping_criterion.f90 @@ -68,6 +68,8 @@ program slsqp_test_stopping_criterion write(*,*) 'istat :', istat write(*,*) 'iterations :', iterations write(*,*) '' + if (any(abs(x - [0.78641515097183889_wp, 0.61769831659541152_wp]) > 1.0e-4_wp)) & + error stop 'Error: incorrect solution' else error stop 'error calling slsqp.' end if From e53bb20b62910a0dfbf15d34a660bbbf8cb084e4 Mon Sep 17 00:00:00 2001 From: Jacob Williams Date: Sun, 4 Jun 2023 17:05:12 -0500 Subject: [PATCH 5/5] minor commenting --- src/slsqp_module.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/slsqp_module.f90 b/src/slsqp_module.f90 index 8fc7652..58ebbae 100644 --- a/src/slsqp_module.f90 +++ b/src/slsqp_module.f90 @@ -173,7 +173,7 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& real(wp),intent(in) :: acc !! accuracy logical,intent(out) :: status_ok !! will be false if there were errors integer,intent(in),optional :: linesearch_mode !! 1 = inexact (default), 2 = exact - integer,intent(in),optional :: iprint !! unit number of status messages (default=output_unit) + integer,intent(in),optional :: iprint !! unit number of status messages (default=`output_unit`) procedure(iterfunc),optional :: report !! user-defined procedure that will be called once per iteration real(wp),intent(in),optional :: alphamin !! minimum alpha for linesearch [default 0.1] real(wp),intent(in),optional :: alphamax !! maximum alpha for linesearch [default 1.0] @@ -183,6 +183,8 @@ subroutine initialize_slsqp(me,n,m,meq,max_iter,acc,f,g,xl,xu,status_ok,& !! * 1 - approximate by basic backward differences !! * 2 - approximate by basic forward differences !! * 3 - approximate by basic central differences + !! + !! Note that modes 1-3 do not respect the variable bounds. real(wp),intent(in),optional :: gradient_delta !! perturbation step size (>epsilon) to compute the approximated !! gradient by finite differences (`gradient_mode` 1-3). !! note that this is an absolute step that does not respect