Skip to content

Commit 955acb5

Browse files
zaikunzhangzaikunzhang
committed
250819.145512.CST fix a bug in preproc.f90 due to the co-exisitence of rhobeg_in and rhobeg_old, rhoend_in and rhoend_old; they were not initialized correctly
1 parent 3468967 commit 955acb5

File tree

4 files changed

+42
-33
lines changed

4 files changed

+42
-33
lines changed

‎fortran/common/preproc.f90‎

Lines changed: 9 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ module preproc_mod
66
!
77
! Started: July 2020
88
!
9-
! Last Modified: Sun 17 Aug 2025 07:47:41 PM CST
9+
! Last Modified: Tue 19 Aug 2025 01:55:20 PM CST
1010
!--------------------------------------------------------------------------------------------------!
1111

1212
! N.B.:
@@ -90,10 +90,9 @@ subroutine preproc(solver, n, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
9090
real(RP) :: gamma2_in
9191
real(RP) :: rhobeg_default
9292
real(RP) :: rhobeg_in
93-
real(RP) :: rhobeg_old
9493
real(RP) :: rhoend_default
9594
real(RP) :: rhoend_in
96-
real(RP) :: x0_old(n)
95+
real(RP) :: x0_in(n)
9796

9897
! Preconditions
9998
if (DEBUGGING) then
@@ -284,6 +283,9 @@ subroutine preproc(solver, n, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
284283

285284
! Validate RHOBEG and RHOEND
286285

286+
rhobeg_in = rhobeg
287+
rhoend_in = rhoend
288+
287289
! Revise the default values for RHOBEG/RHOEND according to the solver.
288290
if (lower(solver) == 'bobyqa') then
289291
rhobeg_default = max(EPS, min(RHOBEG_DFT, minval(xu - xl) / 4.0_RP))
@@ -297,7 +299,6 @@ subroutine preproc(solver, n, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
297299
! Do NOT merge the IF below into the ELSEIF above! Otherwise, XU and XL may be accessed even if
298300
! the solver is not BOBYQA, because the logical evaluation is not short-circuit.
299301
if (rhobeg > minval(xu - xl) / TWO) then
300-
rhobeg_in = rhobeg
301302
! Do NOT make this revision if RHOBEG not positive or not finite, because otherwise RHOBEG
302303
! will get a huge value when XU or XL contains huge values that indicate unbounded variables.
303304
rhobeg = minval(xu - xl) / 4.0_RP ! Here, we do not take RHOBEG_DEFAULT.
@@ -308,7 +309,6 @@ subroutine preproc(solver, n, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
308309
end if
309310

310311
if (.not. (is_finite(rhobeg) .and. rhobeg > 0)) then ! RHOBEG = NaN falls into this case.
311-
rhobeg_in = rhobeg
312312
! Take RHOEND into account if it has a valid value. We do not do this if the solver is BOBYQA,
313313
! which requires that RHOBEG <= (XU-XL)/2.
314314
if (is_finite(rhoend) .and. rhoend > 0 .and. lower(solver) /= 'bobyqa') then
@@ -321,7 +321,6 @@ subroutine preproc(solver, n, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
321321
end if
322322

323323
if (.not. (is_finite(rhoend) .and. rhoend >= 0 .and. rhoend <= rhobeg)) then ! RHOEND = NaN falls into this case.
324-
rhoend_in = rhoend
325324
rhoend = max(EPS, min((RHOEND_DFT / RHOBEG_DFT) * rhobeg, rhoend_default))
326325
call warning(solver, 'Invalid RHOEND: '//num2str(rhoend_in)// &
327326
& '; we should have '//num2str(rhobeg)//' = RHOBEG >= RHOEND >= 0; it is set to '//num2str(rhoend))
@@ -334,7 +333,7 @@ subroutine preproc(solver, n, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
334333
if (lower(solver) == 'bobyqa') then
335334
! Revise X0 if allowed and needed.
336335
if (.not. honour_x0) then
337-
x0_old = x0 ! Recorded to see whether X0 is really revised.
336+
x0_in = x0 ! Recorded to see whether X0 is really revised.
338337
! N.B.: The following revision is valid only if XL <= X0 <= XU and RHOBEG <= MINVAL(XU-XL)/2,
339338
! which should hold at this point due to the revision of RHOBEG and moderation of X0.
340339
! The cases below are mutually exclusive in precise arithmetic as MINVAL(XU-XL) >= 2*RHOBEG.
@@ -358,7 +357,7 @@ subroutine preproc(solver, n, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
358357
!!x0(ubx) = xu(ubx);
359358
!!x0(ubx_minus) = xu(ubx_minus) - rhobeg;
360359

361-
if (any(abs(x0_old - x0) > 0)) then
360+
if (any(abs(x0_in - x0) > 0)) then
362361
call warning(solver, 'X0 is revised so that the distance between X0 and the inactive bounds is at least RHOBEG = '// &
363362
& num2str(rhobeg)//'; revise RHOBEG or set HONOUR_X0 to .TRUE. if you prefer to keep X0 unchanged')
364363
end if
@@ -367,14 +366,13 @@ subroutine preproc(solver, n, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
367366
! Revise RHOBEG if needed.
368367
! N.B.: If X0 has been revised above (i.e., HONOUR_X0 is FALSE), then the following revision
369368
! is unnecessary in precise arithmetic. However, it may still be needed due to rounding errors.
370-
rhobeg_old = rhobeg
371369
lbx = (is_finite(xl) .and. x0 - xl <= EPS * max(ONE, abs(xl))) ! X0 essentially equals XL
372370
ubx = (is_finite(xu) .and. x0 - xu >= -EPS * max(ONE, abs(xu))) ! X0 essentially equals XU
373371
x0(trueloc(lbx)) = xl(trueloc(lbx))
374372
x0(trueloc(ubx)) = xu(trueloc(ubx))
375373
rhobeg = max(EPS, minval([rhobeg, x0(falseloc(lbx)) - xl(falseloc(lbx)), xu(falseloc(ubx)) - x0(falseloc(ubx))]))
376-
if (rhobeg_old - rhobeg > EPS * max(ONE, rhobeg_old)) then
377-
rhoend = max(EPS, min((rhoend / rhobeg_old) * rhobeg, rhoend)) ! We do not revise RHOEND unless RHOBEG is truly revised.
374+
if (rhobeg_in - rhobeg > EPS * max(ONE, rhobeg_in)) then
375+
rhoend = max(EPS, min((rhoend / rhobeg_in) * rhobeg, rhoend)) ! We do not revise RHOEND unless RHOBEG is truly revised.
378376
if (has_rhobeg) then
379377
call warning(solver, 'RHOBEG is revised from '//num2str(rhobeg_in)//' to '//num2str(rhobeg)// &
380378
& ' and RHOEND from '//num2str(rhoend_in)//' to '//num2str(rhoend)// &

‎fortran/tests/makefiles/Makefile.common‎

Lines changed: 25 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -496,16 +496,16 @@ gtest_i2_r4_d1_tst_c gtest_i4_r4_d1_tst_c gtest_i8_r4_d1_tst_c gtest_i2_r4_d0_ts
496496
FCF := $(GFORT) -Wno-function-elimination -ffpe-trap=zero $(GFF)
497497

498498
gtest_i2_r8_d1_tst_c gtest_i4_r8_d1_tst_c gtest_i8_r8_d1_tst_c gtest_i2_r8_d0_tst_c gtest_i4_r8_d0_tst_c gtest_i8_r8_d0_tst_c: \
499-
FCO := $(GFORT) $(FFLAGSO) -Wno-function-elimination -ffpe-trap=zero,overflow#,invalid,underflow,denorm
499+
FCO := $(GFORT) $(FFLAGSO) -Wno-function-elimination -ffpe-trap=zero,overflow,invalid#,underflow,denorm
500500
gtest_i2_r8_d1_tst_c gtest_i4_r8_d1_tst_c gtest_i8_r8_d1_tst_c gtest_i2_r8_d0_tst_c gtest_i4_r8_d0_tst_c gtest_i8_r8_d0_tst_c: \
501-
FCG := $(GFORT) $(FFLAGSG) -ffpe-trap=zero,overflow#,invalid,underflow,denorm
501+
FCG := $(GFORT) $(FFLAGSG) -ffpe-trap=zero,overflow,invalid#,underflow,denorm
502502
gtest_i2_r8_d1_tst_c gtest_i4_r8_d1_tst_c gtest_i8_r8_d1_tst_c gtest_i2_r8_d0_tst_c gtest_i4_r8_d0_tst_c gtest_i8_r8_d0_tst_c: \
503-
FCF := $(GFORT) -Wno-function-elimination -ffpe-trap=zero $(GFF)
503+
FCF := $(GFORT) -Wno-function-elimination -ffpe-trap=zero,invalid $(GFF)
504504

505505
gtest_i2_r16_d1_tst_c gtest_i4_r16_d1_tst_c gtest_i8_r16_d1_tst_c gtest_i2_r16_d0_tst_c gtest_i4_r16_d0_tst_c gtest_i8_r16_d0_tst_c: \
506-
FCO := $(GFORT) $(FFLAGSO) -Wno-function-elimination -ffpe-trap=zero,overflow#,invalid,underflow,denorm
506+
FCO := $(GFORT) $(FFLAGSO) -Wno-function-elimination -ffpe-trap=zero,overflow,invalid#,underflow,denorm
507507
gtest_i2_r16_d1_tst_c gtest_i4_r16_d1_tst_c gtest_i8_r16_d1_tst_c gtest_i2_r16_d0_tst_c gtest_i4_r16_d0_tst_c gtest_i8_r16_d0_tst_c: \
508-
FCG := $(GFORT) $(FFLAGSG) -ffpe-trap=zero,overflow#invalid,underflow,denorm
508+
FCG := $(GFORT) $(FFLAGSG) -ffpe-trap=zero,overflow,invalid#,underflow,denorm
509509
gtest_i2_r16_d1_tst_c gtest_i4_r16_d1_tst_c gtest_i8_r16_d1_tst_c gtest_i2_r16_d0_tst_c gtest_i4_r16_d0_tst_c gtest_i8_r16_d0_tst_c: \
510510
FCF := $(GFORT) -Wno-function-elimination -ffpe-trap=zero,invalid $(GFF)
511511

@@ -554,11 +554,14 @@ itest_i2_r8_d1_tst_c itest_i4_r8_d1_tst_c itest_i8_r8_d1_tst_c itest_i2_r8_d0_ts
554554
FCF := $(IFORT) -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags -fp-trap=divzero,overflow $(IFF)
555555

556556
itest_i2_r16_d1_tst_c itest_i4_r16_d1_tst_c itest_i8_r16_d1_tst_c itest_i2_r16_d0_tst_c itest_i4_r16_d0_tst_c itest_i8_r16_d0_tst_c: \
557-
FCO := $(IFORT) $(FFLAGSO)
557+
FCO := $(IFORT) $(FFLAGSO) -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags \
558+
-fp-trap=divzero,overflow#invalid,underflow,denormal
558559
itest_i2_r16_d1_tst_c itest_i4_r16_d1_tst_c itest_i8_r16_d1_tst_c itest_i2_r16_d0_tst_c itest_i4_r16_d0_tst_c itest_i8_r16_d0_tst_c: \
559-
FCG := $(IFORT) $(FFLAGSG) -check all
560+
FCG := $(IFORT) $(FFLAGSG) -check all -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags \
561+
-fp-trap=divzero,overflow#,invalid,underflow,denormal
560562
itest_i2_r16_d1_tst_c itest_i4_r16_d1_tst_c itest_i8_r16_d1_tst_c itest_i2_r16_d0_tst_c itest_i4_r16_d0_tst_c itest_i8_r16_d0_tst_c: \
561-
FCF := $(IFORT) $(IFF)
563+
FCF := $(IFORT) -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags -fp-trap=divzero,overflow $(IFF)
564+
562565

563566
# NAG nagfor
564567
# In massive tests, we skip the useful -mtrace option (print memory allocation trace), as its output is enormous.
@@ -728,22 +731,30 @@ xtest_i2_r4_d1_tst_c xtest_i4_r4_d1_tst_c xtest_i8_r4_d1_tst_c xtest_i2_r4_d0_ts
728731
FCF := $(XFORT) -fp-trap=divzero $(XFF)
729732

730733
xtest_i2_r8_d1_tst_c xtest_i4_r8_d1_tst_c xtest_i8_r8_d1_tst_c xtest_i2_r8_d0_tst_c xtest_i4_r8_d0_tst_c xtest_i8_r8_d0_tst_c: \
731-
FCO := $(XFORT) $(FFLAGSO) -ftrapuv -fp-trap=divzero#,overflow#,invalid,underflow,denormal
734+
FCO := $(XFORT) $(FFLAGSO) -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags \
735+
-fp-trap=divzero,overflow#invalid,underflow,denormal
732736
#-no-ftz
733737
xtest_i2_r8_d1_tst_c xtest_i4_r8_d1_tst_c xtest_i8_r8_d1_tst_c xtest_i2_r8_d0_tst_c xtest_i4_r8_d0_tst_c xtest_i8_r8_d0_tst_c: \
734-
FCG := $(XFORT) $(FFLAGSG) -check all,nouninit -ftrapuv -fp-trap=divzero#,overflow#,invalid,underflow,denormal
738+
FCG := $(XFORT) $(FFLAGSG) -check all,nouninit -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags \
739+
-fp-trap=divzero,overflow#invalid,underflow,denormal
735740
#-no-ftz
736741
xtest_i2_r8_d1_tst_c xtest_i4_r8_d1_tst_c xtest_i8_r8_d1_tst_c xtest_i2_r8_d0_tst_c xtest_i4_r8_d0_tst_c xtest_i8_r8_d0_tst_c: \
737-
FCF := $(XFORT) -ftrapuv -fp-trap=divzero $(XFF)
742+
FCF := $(XFORT) -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags \
743+
-fp-trap=divzero,overflow $(XFF)
744+
738745

739746
xtest_i2_r16_d1_tst_c xtest_i4_r16_d1_tst_c xtest_i8_r16_d1_tst_c xtest_i2_r16_d0_tst_c xtest_i4_r16_d0_tst_c xtest_i8_r16_d0_tst_c: \
740-
FCO := $(XFORT) $(FFLAGSO) -ftrapuv -fp-trap=divzero#,overflow#,invalid,underflow,denormal
747+
FCO := $(XFORT) $(FFLAGSO) -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags \
748+
-fp-trap=divzero,overflow#invalid,underflow,denormal
741749
#-no-ftz
742750
xtest_i2_r16_d1_tst_c xtest_i4_r16_d1_tst_c xtest_i8_r16_d1_tst_c xtest_i2_r16_d0_tst_c xtest_i4_r16_d0_tst_c xtest_i8_r16_d0_tst_c: \
743-
FCG := $(XFORT) $(FFLAGSG) -check all,nouninit -ftrapuv -fp-trap=divzero#,overflow#invalid,underflow,denormal
751+
FCG := $(XFORT) $(FFLAGSG) -check all,nouninit -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags \
752+
-fp-trap=divzero,overflow#invalid,underflow,denormal
744753
#-no-ftz
745754
xtest_i2_r16_d1_tst_c xtest_i4_r16_d1_tst_c xtest_i8_r16_d1_tst_c xtest_i2_r16_d0_tst_c xtest_i4_r16_d0_tst_c xtest_i8_r16_d0_tst_c: \
746-
FCF := $(XFORT) -ftrapuv -fp-trap=divzero $(XFF)
755+
FCF := $(XFORT) -ftrapuv -fpe0 -fpe-all=0 -assume ieee_fpe_flags \
756+
-fp-trap=divzero,overflow $(XFF)
757+
747758

748759
# G95
749760
# G95 supports f2003, but not higher.

‎matlab/interfaces/private/preprima.m‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1983,7 +1983,7 @@
19831983
if isfield(options, 'honour_x0') && ~options.honour_x0
19841984
% N.B.: The following code is valid only if lb <= x0 <= ub and rhobeg <= min(ub-lb)/2, which
19851985
% hold after `pre_options` and `project` are invoked.
1986-
x0_old = x0;
1986+
x0_in = x0;
19871987
lbx = (x0 <= lb + 0.5*options.rhobeg);
19881988
lbx_plus = (x0 > lb + 0.5*options.rhobeg) & (x0 < lb + options.rhobeg);
19891989
ubx_minus = (x0 < ub - 0.5*options.rhobeg) & (x0 > ub - options.rhobeg);
@@ -1992,7 +1992,7 @@
19921992
x0(lbx_plus) = lb(lbx_plus) + options.rhobeg;
19931993
x0(ubx_minus) = ub(ubx_minus) - options.rhobeg;
19941994
x0(ubx) = ub(ubx);
1995-
if any(abs(x0_old-x0) > 0)
1995+
if any(abs(x0_in-x0) > 0)
19961996
wid = sprintf('%s:ReviseX0', invoker);
19971997
wmsg = sprintf('%s: x0 is revised so that the distance between x0 and the inactive bounds is at least rhobeg; set options.honour_x0 to true if you prefer to keep x0.', invoker);
19981998
warning(wid, '%s', wmsg);
@@ -2003,14 +2003,14 @@
20032003
% Revise rhobeg if needed.
20042004
% N.B.: If x0 has been revised above (i.e., options.honour_x0 is false), then the following revision
20052005
% is unnecessary in precise arithmetic. However, it may still be needed due to rounding errors.
2006-
rhobeg_old = options.rhobeg;
2006+
rhobeg_in = options.rhobeg;
20072007
rho_ratio = options.rhoend / options.rhobeg;
20082008
lbx = (lb > -inf & x0 - lb <= eps*max(abs(lb), 1)); % x0 essentially equals lb
20092009
ubx = (ub < inf & x0 - ub >= -eps*max(abs(ub), 1)); % x0 essentially equals ub
20102010
x0(lbx) = lb(lbx);
20112011
x0(ubx) = ub(ubx);
20122012
options.rhobeg = max(eps, min([options.rhobeg; x0(~lbx) - lb(~lbx); ub(~ubx) - x0(~ubx)]));
2013-
if rhobeg_old - options.rhobeg > eps*max(1, rhobeg_old)
2013+
if rhobeg_in - options.rhobeg > eps*max(1, rhobeg_in)
20142014
options.rhoend = max(eps, min(rho_ratio*options.rhobeg, options.rhoend)); % We do not revise rhoend unless rhobeg is revised
20152015
if ismember('rhobeg', user_options_fields) || ismember('rhoend', user_options_fields)
20162016
wid = sprintf('%s:ReviseRhobeg', invoker);

‎pyprima/src/pyprima/common/preproc.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -205,20 +205,20 @@ def preproc(solver, num_vars, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
205205
# least RHOBEG. If HONOUR_X0 == TRUE, revise RHOBEG if needed; otherwise, revise HONOUR_X0 if needed.
206206
if present(honour_x0):
207207
if honour_x0:
208-
rhobeg_old = rhobeg;
208+
rhobeg_in = rhobeg;
209209
lbx = np.isfinite(xl) & (x0 - xl <= EPS * np.maximum(1, np.abs(xl))) # X0 essentially equals XL
210210
ubx = np.isfinite(xu) & (x0 - xu >= -EPS * np.maximum(1, np.abs(xu))) # X0 essentially equals XU
211211
x0[lbx] = xl[lbx]
212212
x0[ubx] = xu[ubx]
213213
rhobeg = max(EPS, np.min([rhobeg, x0[~lbx] - xl[~lbx], xu[~ubx] - x0[~ubx]]))
214-
if rhobeg_old - rhobeg > EPS * max(1, rhobeg_old):
214+
if rhobeg_in - rhobeg > EPS * max(1, rhobeg_in):
215215
rhoend = max(EPS, min(0.1 * rhobeg, rhoend)) # We do not revise RHOEND unless RHOBEG is truly revised.
216216
if has_rhobeg:
217217
warn(f'{solver}: RHOBEG is revised to {rhobeg} and RHOEND to at most 0.1*RHOBEG so that the distance between X0 and the inactive bounds is at least RHOBEG')
218218
else:
219219
rhoend = np.minimum(rhoend, rhobeg) # This may update RHOEND slightly.
220220
else:
221-
x0_old = x0 # Recorded to see whether X0 is really revised.
221+
x0_in = x0 # Recorded to see whether X0 is really revised.
222222
# N.B.: The following revision is valid only if XL <= X0 <= XU and RHOBEG <= MINVAL(XU-XL)/2,
223223
# which should hold at this point due to the revision of RHOBEG and moderation of X0.
224224
# The cases below are mutually exclusive in precise arithmetic as MINVAL(XU-XL) >= 2*RHOBEG.
@@ -231,7 +231,7 @@ def preproc(solver, num_vars, iprint, maxfun, maxhist, ftarget, rhobeg, rhoend,
231231
x0[ubx] = xu[ubx]
232232
x0[ubx_minus] = xu[ubx_minus] - rhobeg
233233

234-
if (any(np.abs(x0_old - x0) > 0)):
234+
if (any(np.abs(x0_in - x0) > 0)):
235235
warn(f'{solver}: X0 is revised so that the distance between X0 and the inactive bounds is at least RHOBEG set HONOUR_X0 to .TRUE. if you prefer to keep X0 unchanged')
236236

237237
# Validate CTOL (it can be 0)

0 commit comments

Comments
 (0)