Refactoring (Part 2)
I'm working on a modern Fortran version of the simulated annealing optimization algorithm, using this code as a starting point. There is a Fortran 90 version here, but this seems to be mostly just a straightforward conversion of the original code to free-form source formatting. I have in mind a more thorough modernization and the addition of some new features that I need.
Refactoring old Fortran code is always fun. Consider the following snippet:
C Check termination criteria.
QUIT = .FALSE.
FSTAR(1) = F
IF ((FOPT - FSTAR(1)) .LE. EPS) QUIT = .TRUE.
DO 410, I = 1, NEPS
IF (ABS(F - FSTAR(I)) .GT. EPS) QUIT = .FALSE.
410 CONTINUE
This can be greatly simplified to:
! check termination criteria.
fstar(1) = f
quit = ((fopt-f)<=eps) .and. (all(abs(f-fstar)<=eps)
Note that we are using some of the vector features of modern Fortran to remove the loop. Consider also this function in the original code:
FUNCTION EXPREP(RDUM)
C This function replaces exp to avoid under- and overflows and is
C designed for IBM 370 type machines. It may be necessary to modify
C it for other machines. Note that the maximum and minimum values of
C EXPREP are such that they has no effect on the algorithm.
DOUBLE PRECISION RDUM, EXPREP
IF (RDUM .GT. 174.) THEN
EXPREP = 3.69D+75
ELSE IF (RDUM .LT. -180.) THEN
EXPREP = 0.0
ELSE
EXPREP = EXP(RDUM)
END IF
RETURN
END
It is somewhat disturbing that the comments mention IBM 370 machines. This routine is unchanged in the f90 version. However, these numeric values are no longer correct for modern hardware. Using modern Fortran, we can write a totally portable version of this routine like so:
pure function exprep(x) result(f)
use, intrinsic :: ieee_exceptions
implicit none
real(dp), intent(in) :: x
real(dp) :: f
logical,dimension(2) :: flags
type(ieee_flag_type),parameter,dimension(2) :: out_of_range = [ieee_overflow,ieee_underflow]
call ieee_set_halting_mode(out_of_range,.false.)
f = exp(x)
call ieee_get_flag(out_of_range,flags)
if (any(flags)) then
call ieee_set_flag(out_of_range,.false.)
if (flags(1)) then
f = huge(1.0_dp)
else
f = 0.0_dp
end if
end if
end function exprep
This new version is entirely portable and works for any real kind. It contains no magic numbers and uses the intrinsic IEEE exceptions module of modern Fortran to protect against overflow and underflow, and the huge() intrinsic to return the largest real(dp)
number.
So, stay tuned...
See also
- S. Kirkpatrick, C. D. Gelatt Jr., M. P. Vecchi, "Optimization by Simulated Annealing", Science 13 May 1983, Vol. 220, Issue 4598, pp. 671-680
- W. L. Goffe, SIMANN: A Global Optimization Algorithm using Simulated Annealing, Studies in Nonlinear Dynamics & Econometrics, De Gruyter, vol. 1(3), pages 1-9, October 1996.
- Original simann.f source code [Netlib]