diff --git a/fortran/bobyqa/bobyqb.f90 b/fortran/bobyqa/bobyqb.f90 index 9d14c9d2a8..694a2be246 100644 --- a/fortran/bobyqa/bobyqb.f90 +++ b/fortran/bobyqa/bobyqb.f90 @@ -330,7 +330,7 @@ subroutine bobyqb(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm ! Set QRED to the reduction of the quadratic model when the move D is made from XOPT. QRED ! should be positive. If it is nonpositive due to rounding errors, we will not take this step. qred = -quadinc(d, xpt, gopt, pq, hq) ! QRED = Q(XOPT) - Q(XOPT + D) - trfail = (.not. qred > 1.0E-6 * rho**2) ! QRED is tiny/negative or NaN. + trfail = (.not. qred > 1.0E-6_RP * rho**2) ! QRED is tiny/negative or NaN. ! When D is short, make a choice between reducing RHO and improving the geometry depending ! on whether or not our work with the current RHO seems complete. RHO is reduced if the diff --git a/fortran/bobyqa/geometry.f90 b/fortran/bobyqa/geometry.f90 index 863d7a16fc..118df98a2d 100644 --- a/fortran/bobyqa/geometry.f90 +++ b/fortran/bobyqa/geometry.f90 @@ -491,9 +491,9 @@ function geostep(knew, kopt, bmat, delbar, sl, su, xpt, zmat) result(d) ! How to make this condition adaptive? A naive idea is to replace the thresholds to, ! e.g.,1.0E-2*RHOBEG. However, in a test on 20220517, this adaptation worsened the performance. In ! such a test, RHOBEG must take a value that is quite different from one. We tried RHOBEG = 0.9E-2. -!if (delbar > 1.0E-3) then -!if (delbar > 1.0E-1) then -if (delbar > 1.0E-2) then +!if (delbar > 1.0E-3_RP) then +!if (delbar > 1.0E-1_RP) then +if (delbar > 1.0E-2_RP) then return end if !--------------------------------------------------------------------------------------------------! diff --git a/fortran/bobyqa/rescue.f90 b/fortran/bobyqa/rescue.f90 index 382d47fb5b..14b1b9366c 100644 --- a/fortran/bobyqa/rescue.f90 +++ b/fortran/bobyqa/rescue.f90 @@ -499,7 +499,7 @@ subroutine rescue(calfun, solver, iprint, maxfun, delta, ftarget, xl, xu, kopt, ! Skipping an XNEW that is close but not identical to XPT(:, KPT) will cause discrepancy ! between [BMAT, ZMAT] and XPT, since the former has been updated, but it is not severe as ! the difference between XNEW and XPT(:, KPT) is tiny. - if (sum(abs(xnew - xpt(:, kpt))) <= 1.0E-2 * delta .or. .not. is_finite(sum(abs(xnew)))) then + if (sum(abs(xnew - xpt(:, kpt))) <= 1.0E-2_RP * delta .or. .not. is_finite(sum(abs(xnew)))) then cycle end if xpt(:, kpt) = xnew @@ -747,7 +747,7 @@ subroutine updateh_rsc(knew, beta, vlag_in, bmat, zmat, info) ! Apply Givens rotations to put zeros in the KNEW-th row of ZMAT. After this, ZMAT(KNEW, :) contains ! only one nonzero at ZMAT(KNEW, 1). Entries of ZMAT are treated as 0 if the moduli are quite small. do j = 2, npt - n - 1_IK - if (abs(zmat(knew, j)) > 1.0E-20 * maxval(abs(zmat))) then ! This threshold is by Powell + if (abs(zmat(knew, j)) > 1.0E-20_RP * maxval(abs(zmat))) then ! This threshold is by Powell grot = planerot(zmat(knew, [1_IK, j])) zmat(:, [1_IK, j]) = matprod(zmat(:, [1_IK, j]), transpose(grot)) end if diff --git a/fortran/bobyqa/update.f90 b/fortran/bobyqa/update.f90 index 7e49fec3f9..eb32d5f6c3 100644 --- a/fortran/bobyqa/update.f90 +++ b/fortran/bobyqa/update.f90 @@ -153,7 +153,7 @@ subroutine updateh(knew, kopt, d, xpt, bmat, zmat, info) ! Apply Givens rotations to put zeros in the KNEW-th row of ZMAT. After this, ZMAT(KNEW, :) contains ! only one nonzero at ZMAT(KNEW, 1). Entries of ZMAT are treated as 0 if the moduli are quite small. do j = 2, npt - n - 1_IK - if (abs(zmat(knew, j)) > 1.0E-20 * maxval(abs(zmat))) then ! This threshold is by Powell + if (abs(zmat(knew, j)) > 1.0E-20_RP * maxval(abs(zmat))) then ! This threshold is by Powell grot = planerot(zmat(knew, [1_IK, j])) zmat(:, [1_IK, j]) = matprod(zmat(:, [1_IK, j]), transpose(grot)) end if diff --git a/fortran/cobyla/cobylb.f90 b/fortran/cobyla/cobylb.f90 index f64abb7db9..9507f914fa 100644 --- a/fortran/cobyla/cobylb.f90 +++ b/fortran/cobyla/cobylb.f90 @@ -368,7 +368,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et ! Evaluate PREREM, which is the predicted reduction in the merit function. ! In theory, PREREM >= 0 and it is 0 iff CPEN = 0 = PREREF. This may not be true numerically. prerem = preref + cpen * prerec - trfail = (.not. prerem > 1.0E-6 * min(cpen, ONE) * rho) ! PREREM is tiny/negative or NaN. + trfail = (.not. prerem > 1.0E-6_RP * min(cpen, ONE) * rho) ! PREREM is tiny/negative or NaN. if (shortd .or. trfail) then ! Reduce DELTA if D is short or D fails to render PREREM > 0. The latter can happen due to @@ -387,7 +387,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et distsq(1:n) = [(sum((x - (sim(:, n + 1) + sim(:, j)))**2), j=1, n)] ! Implied do-loop !!MATLAB: distsq(1:n) = sum((x - (sim(:,1:n) + sim(:, n+1)))**2, 1) % Implicit expansion j = int(minloc(distsq, dim=1), kind(j)) - if (distsq(j) <= (1.0E-4 * rhoend)**2) then + if (distsq(j) <= (1.0E-4_RP * rhoend)**2) then f = fval(j) constr = conmat(:, j) cstrv = cval(j) @@ -582,7 +582,7 @@ subroutine cobylb(calcfc, iprint, maxfilt, maxfun, amat, bvec, ctol, cweight, et distsq(1:n) = [(sum((x - (sim(:, n + 1) + sim(:, j)))**2), j=1, n)] ! Implied do-loop !!MATLAB: distsq(1:n) = sum((x - (sim(:,1:n) + sim(:, n+1)))**2, 1) % Implicit expansion j = int(minloc(distsq, dim=1), kind(j)) - if (distsq(j) <= (1.0E-4 * rhoend)**2) then + if (distsq(j) <= (1.0E-4_RP * rhoend)**2) then f = fval(j) constr = conmat(:, j) cstrv = cval(j) diff --git a/fortran/common/powalg.f90 b/fortran/common/powalg.f90 index 646b2ffbd7..2bed8d1b22 100644 --- a/fortran/common/powalg.f90 +++ b/fortran/common/powalg.f90 @@ -1272,7 +1272,7 @@ subroutine updateh(knew, kref, d, xpt, idz, bmat, zmat, info) ! Powell's condition in NEWUOA/LINCOA for the IF ... THEN below: IF (ZMAT(KNEW, J) /= 0) THEN ! A possible alternative: IF (ABS(ZMAT(KNEW, J)) > 1.0E-20 * ABS(ZMAT(KNEW, JL))) THEN - if (abs(zmat(knew, j)) > 1.0E-20 * maxval(abs(zmat))) then ! Threshold comes from Powell's BOBYQA + if (abs(zmat(knew, j)) > 1.0E-20_RP * maxval(abs(zmat))) then ! Threshold comes from Powell's BOBYQA ! Multiply a Givens rotation to ZMAT from the right so that ZMAT(KNEW, [JL,J]) becomes [*,0]. grot = planerot(zmat(knew, [jl, j])) !!MATLAB: grot = planerot(zmat(knew, [jl, j])') zmat(:, [jl, j]) = matprod(zmat(:, [jl, j]), transpose(grot)) diff --git a/fortran/common/string.f90 b/fortran/common/string.f90 index 1ee89be358..3fb2afd66c 100644 --- a/fortran/common/string.f90 +++ b/fortran/common/string.f90 @@ -151,7 +151,7 @@ function real2str_scalar(x, ndgt, nexp) result(s) if (present(nexp)) then nexp_loc = nexp else - nexp_loc = ceiling(log10(real(range(x) + 0.1))) ! Use + 0.1 in case RANGE(X) = 10^k. + nexp_loc = ceiling(log10(real(range(x) + 0.1_RP))) ! Use + 0.1 in case RANGE(X) = 10^k. end if nexp_loc = min(nexp_loc, floor(real(MAX_NUM_STR_LEN - 5) / 2.0)) @@ -254,7 +254,7 @@ function real2str_vector(x, ndgt, nexp, nx) result(s) if (present(nexp)) then nexp_loc = nexp else - nexp_loc = ceiling(log10(real(range(x)) + 0.1)) ! Use + 0.1 in case RANGE(X) = 10^k. + nexp_loc = ceiling(log10(real(range(x)) + 0.1_RP)) ! Use + 0.1 in case RANGE(X) = 10^k. end if nexp_loc = min(nexp_loc, floor(real(MAX_NUM_STR_LEN - 5) / 2.0)) diff --git a/fortran/lincoa/lincob.f90 b/fortran/lincoa/lincob.f90 index a148ce94f3..c2d2630e4f 100644 --- a/fortran/lincoa/lincob.f90 +++ b/fortran/lincoa/lincob.f90 @@ -429,7 +429,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b ! Set QRED to the reduction of the quadratic model when the move D is made from XOPT. QRED ! should be positive. If it is nonpositive due to rounding errors, we will not take this step. qred = -quadinc(d, xpt, gopt, pq, hq) ! QRED = Q(XOPT) - Q(XOPT + D) - trfail = (.not. qred > 1.0E-6 * rho**2) ! QRED is tiny/negative or NaN. + trfail = (.not. qred > 1.0E-6_RP * rho**2) ! QRED is tiny/negative or NaN. if (shortd .or. trfail) then ! In this case, do nothing but reducing DELTA. Afterward, DELTA < DNORM may occur. @@ -534,7 +534,7 @@ subroutine lincob(calfun, iprint, maxfilt, maxfun, npt, Aeq, Aineq, amat, beq, b ! verify a curvature condition that really indicates that recent models are sufficiently ! accurate. Here, however, we are not really sure whether they are accurate or not. Therefore, ! ACCURATE_MOD is not the best name, but we keep it to align with the other solvers. - accurate_mod = all(dnorm_rec <= rho) .or. all(dnorm_rec(2:size(dnorm_rec)) <= 0.2 * rho) + accurate_mod = all(dnorm_rec <= rho) .or. all(dnorm_rec(2:size(dnorm_rec)) <= 0.2_RP * rho) ! Powell's version (note that size(dnorm_rec) = 5 in his implementation): !accurate_mod = all(dnorm_rec <= HALF * rho) .or. all(dnorm_rec(3:size(dnorm_rec)) <= TENTH * rho) ! CLOSE_ITPSET: Are the interpolation points close to XOPT? diff --git a/fortran/newuoa/newuob.f90 b/fortran/newuoa/newuob.f90 index 1f96ac12d1..f744453fea 100644 --- a/fortran/newuoa/newuob.f90 +++ b/fortran/newuoa/newuob.f90 @@ -287,7 +287,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm ! Set QRED to the reduction of the quadratic model when the move D is made from XOPT. QRED ! should be positive. If it is nonpositive due to rounding errors, we will not take this step. qred = -quadinc(d, xpt, gopt, pq, hq) - trfail = (.not. qred > 1.0E-6 * rho**2) ! QRED is tiny/negative, or NaN. + trfail = (.not. qred > 1.0E-6_RP * rho**2) ! QRED is tiny/negative, or NaN. if (shortd .or. trfail) then ! In this case, do nothing but reducing DELTA. Afterward, DELTA < DNORM may occur. @@ -306,7 +306,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm distsq = [(sum((x - (xbase + xpt(:, k)))**2, dim=1), k=1, npt)] ! Implied do-loop !!MATLAB: distsq = sum((x - (xbase + xpt))**2, 1) % Implicit expansion k = int(minloc(distsq, dim=1), kind(k)) - if (distsq(k) <= (1.0E-3 * rhoend)**2) then + if (distsq(k) <= (1.0E-3_RP * rhoend)**2) then f = fval(k) else ! Evaluate the objective function at X, taking care of possible Inf/NaN values. @@ -553,7 +553,7 @@ subroutine newuob(calfun, iprint, maxfun, npt, eta1, eta2, ftarget, gamma1, gamm distsq = [(sum((x - (xbase + xpt(:, k)))**2, dim=1), k=1, npt)] ! Implied do-loop !!MATLAB: distsq = sum((x - (xbase + xpt))**2, 1) % Implicit expansion k = int(minloc(distsq, dim=1), kind(k)) - if (distsq(k) <= (1.0E-3 * rhoend)**2) then + if (distsq(k) <= (1.0E-3_RP * rhoend)**2) then f = fval(k) else ! Evaluate the objective function at X, taking care of possible Inf/NaN values. diff --git a/fortran/tests/test_bobyqa.f90 b/fortran/tests/test_bobyqa.f90 index 4bcf756441..baae61f114 100644 --- a/fortran/tests/test_bobyqa.f90 +++ b/fortran/tests/test_bobyqa.f90 @@ -238,10 +238,10 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) iprint = int(randn(), kind(iprint)) maxfun = int(1.0E2_RP * rand() * real(n, RP), kind(maxfun)) maxhist = int(TWO * rand() * real(max(10_IK * n, maxfun), RP), kind(maxhist)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then maxhist = -maxhist end if - if (rand() <= 0.8) then + if (rand() <= 0.8_RP) then ftarget = -TEN**abs(real(min(range(ftarget), 12), RP) * rand()) elseif (rand() <= 0.5) then ! Note that the value of rand() changes. ftarget = REALMAX @@ -251,9 +251,9 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) rhobeg = noisy(prob % Delta0) rhoend = max(1.0E-5_RP, rhobeg * 10.0_RP**(6.0_RP * rand() - 5.0_RP)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then rhoend = rhobeg - elseif (rand() <= 0.1) then ! Note that the value of rand() changes. + elseif (rand() <= 0.1_RP) then ! Note that the value of rand() changes. rhobeg = ZERO end if diff --git a/fortran/tests/test_cobyla.f90 b/fortran/tests/test_cobyla.f90 index 2d79ea6250..dc996e254a 100644 --- a/fortran/tests/test_cobyla.f90 +++ b/fortran/tests/test_cobyla.f90 @@ -262,21 +262,21 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) iprint = int(randn(), kind(iprint)) maxfun = int(1.0E2_RP * rand() * real(n, RP), kind(maxfun)) maxhist = int(TWO * rand() * real(max(10_IK * n, maxfun), RP), kind(maxhist)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then maxhist = -maxhist end if maxfilt = int(TWO * rand() * real(maxfun, RP), kind(maxfilt)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then maxfilt = 0 end if - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then ctol = randn() * TEN**(-abs(TWO * randn())) - elseif (rand() <= 0.1) then ! Note that the value of rand() changes. + elseif (rand() <= 0.1_RP) then ! Note that the value of rand() changes. ctol = REALMAX else ctol = ZERO end if - if (rand() <= 0.8) then + if (rand() <= 0.8_RP) then ftarget = -TEN**abs(real(min(range(ftarget), 12), RP) * rand()) elseif (rand() <= 0.5) then ! Note that the value of rand() changes. ftarget = REALMAX @@ -286,9 +286,9 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) rhobeg = noisy(prob % Delta0) rhoend = max(1.0E-5_RP, rhobeg * 10.0_RP**(6.0_RP * rand() - 5.0_RP)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then rhoend = rhobeg - elseif (rand() <= 0.1) then ! Note that the value of rand() changes. + elseif (rand() <= 0.1_RP) then ! Note that the value of rand() changes. rhobeg = ZERO end if diff --git a/fortran/tests/test_lincoa.f90 b/fortran/tests/test_lincoa.f90 index 4f46790953..402852e9be 100644 --- a/fortran/tests/test_lincoa.f90 +++ b/fortran/tests/test_lincoa.f90 @@ -303,21 +303,21 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) iprint = int(randn(), kind(iprint)) maxfun = int(1.0E2_RP * rand() * real(n, RP), kind(maxfun)) maxhist = int(TWO * rand() * real(max(10_IK * n, maxfun), RP), kind(maxhist)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then maxhist = -maxhist end if maxfilt = int(TWO * rand() * real(maxfun, RP), kind(maxfilt)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then maxfilt = 0 end if - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then ctol = randn() * TEN**(-abs(TWO * randn())) - elseif (rand() <= 0.1) then ! Note that the value of rand() changes. + elseif (rand() <= 0.1_RP) then ! Note that the value of rand() changes. ctol = REALMAX else ctol = ZERO end if - if (rand() <= 0.8) then + if (rand() <= 0.8_RP) then ftarget = -TEN**abs(real(min(range(ftarget), 12), RP) * rand()) elseif (rand() <= 0.5) then ! Note that the value of rand() changes. ftarget = REALMAX @@ -327,9 +327,9 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) rhobeg = noisy(prob % Delta0) rhoend = max(1.0E-5_RP, rhobeg * 10.0_RP**(6.0_RP * rand() - 5.0_RP)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then rhoend = rhobeg - elseif (rand() <= 0.1) then ! Note that the value of rand() changes. + elseif (rand() <= 0.1_RP) then ! Note that the value of rand() changes. rhobeg = ZERO end if diff --git a/fortran/tests/test_newuoa.f90 b/fortran/tests/test_newuoa.f90 index 8f90dd845d..9451ea0cde 100644 --- a/fortran/tests/test_newuoa.f90 +++ b/fortran/tests/test_newuoa.f90 @@ -215,10 +215,10 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) iprint = int(randn(), kind(iprint)) maxfun = int(1.0E2_RP * rand() * real(n, RP), kind(maxfun)) maxhist = int(TWO * rand() * real(max(10_IK * n, maxfun), RP), kind(maxhist)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then maxhist = -maxhist end if - if (rand() <= 0.8) then + if (rand() <= 0.8_RP) then ftarget = -TEN**abs(real(min(range(ftarget), 12), RP) * rand()) elseif (rand() <= 0.5) then ! Note that the value of rand() changes. ftarget = REALMAX @@ -228,9 +228,9 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) rhobeg = noisy(prob % Delta0) rhoend = max(1.0E-5_RP, rhobeg * 10.0_RP**(6.0_RP * rand() - 5.0_RP)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then rhoend = rhobeg - elseif (rand() <= 0.1) then ! Note that the value of rand() changes. + elseif (rand() <= 0.1_RP) then ! Note that the value of rand() changes. rhobeg = ZERO end if diff --git a/fortran/tests/test_uobyqa.f90 b/fortran/tests/test_uobyqa.f90 index 1291b356f6..a1dc92cba8 100644 --- a/fortran/tests/test_uobyqa.f90 +++ b/fortran/tests/test_uobyqa.f90 @@ -199,10 +199,10 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) iprint = int(randn(), kind(iprint)) maxfun = int(1.0E2_RP * rand() * real(n, RP), kind(maxfun)) maxhist = int(TWO * rand() * real(max(10_IK * n, maxfun), RP), kind(maxhist)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then maxhist = -maxhist end if - if (rand() <= 0.8) then + if (rand() <= 0.8_RP) then ftarget = -TEN**abs(real(min(range(ftarget), 12), RP) * rand()) elseif (rand() <= 0.5) then ! Note that the value of rand() changes. ftarget = REALMAX @@ -212,9 +212,9 @@ subroutine test_solver(probs, mindim, maxdim, nrand, randseed, testdim) rhobeg = noisy(prob % Delta0) rhoend = max(1.0E-5_RP, rhobeg * 10.0_RP**(6.0_RP * rand() - 5.0_RP)) - if (rand() <= 0.1) then + if (rand() <= 0.1_RP) then rhoend = rhobeg - elseif (rand() <= 0.1) then ! Note that the value of rand() changes. + elseif (rand() <= 0.1_RP) then ! Note that the value of rand() changes. rhobeg = ZERO end if diff --git a/fortran/uobyqa/uobyqb.f90 b/fortran/uobyqa/uobyqb.f90 index 2d8eeed04f..6ce0ed58e6 100644 --- a/fortran/uobyqa/uobyqb.f90 +++ b/fortran/uobyqa/uobyqb.f90 @@ -294,7 +294,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r ! Set QRED to the reduction of the quadratic model when the move D is made from XOPT. QRED ! should be positive. If it is nonpositive due to rounding errors, we will not take this step. qred = -quadinc(pq, d, xpt(:, kopt)) ! QRED = Q(XOPT) - Q(XOPT + D) - trfail = (.not. qred > 1.0E-6 * rho**2) ! QRED is tiny/negative or NaN. + trfail = (.not. qred > 1.0E-6_RP * rho**2) ! QRED is tiny/negative or NaN. if (shortd .or. trfail) then ! Powell's code does not reduce DELTA as follows. This comes from NEWUOA and works well. @@ -310,7 +310,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r distsq = [(sum((x - (xbase + xpt(:, k)))**2, dim=1), k=1, npt)] ! Implied do-loop !!MATLAB: distsq = sum((x - (xbase + xpt))**2, 1) % Implicit expansion k = int(minloc(distsq, dim=1), kind(k)) - if (distsq(k) <= (1.0E-4 * rhoend)**2) then + if (distsq(k) <= (1.0E-4_RP * rhoend)**2) then f = fval(k) else ! Evaluate the objective function at X, taking care of possible Inf/NaN values. @@ -474,7 +474,7 @@ subroutine uobyqb(calfun, iprint, maxfun, eta1, eta2, ftarget, gamma1, gamma2, r distsq = [(sum((x - (xbase + xpt(:, k)))**2, dim=1), k=1, npt)] ! Implied do-loop !!MATLAB: distsq = sum((x - (xbase + xpt))**2, 1) % Implicit expansion k = int(minloc(distsq, dim=1), kind(k)) - if (distsq(k) <= (1.0E-4 * rhoend)**2) then + if (distsq(k) <= (1.0E-4_RP * rhoend)**2) then f = fval(k) else ! Evaluate the objective function at X, taking care of possible Inf/NaN values.