Consider this code.
A statement function:
diff(u,v)= u-v
and an if statement:
IF(diff(hmax+factor*h(lmax),hmax).GT.ZERO)
So, it seems what is happening is a test for:
hmax + factor*h(lmax) - hmax > zero
which reduces to:
factor*h(lmax) > zero
I suspect there is some deep numerical reason for doing it the old way? But what is it?
hmax + factor*h(lmax) - hmax >= zero tests whether factor*h(lmax) is less than eps(hmax)
Yeah, so the function was probably a way to prevent the compiler from optimizing this away I guess?
So in modern times, we could probably use:
if (factor*h(lmax) < epsilon(hmax)) ...
It could be just sloppy math and/or programming. But it could be looking for when factor*h(lmax) is small relative to hmax. If this is from the f77 era or before, then it was difficult to get information such as epsilon() and nearest(), so programmers resorted to these kind of tricks to test for convergence.
Timely, we are almost done implementing statement functions in LFortran: Added support for ```statement functions``` by Pranavchiku · Pull Request #911 · lfortran/lfortran · GitHub.
The line with the IF statement also occurs in the TOMS 587 code and the book “Solving Least Squares Problems” by the same authors, Lawson and Hanson. In the TOMS 587 code, the function DIFF is an external function, presumably to hide its details from the compiler.
More precisely, I think, whether factor*h(lmax) is >= hmax * epsilon(hmax), assuming that none of the quantities are negative.
I wonder if any current compiler actually can detect this y+x-y trick and optimize it away (defeating the original intent).
hopefully it doesn’t without a fastmath flag.
Not yet, I think. That would require that Fortran compilers become capable of algebraic reasoning. I tried the following contraption to detect if they did:
logical :: lxpr1, lxpr2
...
lxpr1 = diff(hmax + factor*h(lmax), hmax) .gt. zero
lxpr2 = factor*h(lmax) > hmax*epsilon(0.0d0)
if(lxpr1.neqv.lxpr2)print *,'DIAG 3: ',factor,h(lmax),hmax,lxpr1,lxpr2
and likewise for the other two locations in SLSQP where such tests are performed. Intel Fortran and Gfortran, with maximum optimization levels specified, did not fall into the trap.
I then replaced the DIFF statement function with a preprocessor replacement for the ASF:
#define diff(u,v) ((u)-(v))
Again, no error.
I think the outside set of parentheses is important in that macro, or in the expression written out in fortran. Without the outside parentheses, the compiler can indeed rearrange the expression and eliminate the hmax-hmax term from the expression. K&R C (i.e. before C89) could do the same rearrangement even with the parentheses in the expression.
Fortran compilers have been doing these kinds of optimizations since at least the 1970s. These kinds of rearrangements are discussed explicitly in the fortran standard. I remember reading IBM fortran documentation in the late 1970s that discussed this, and these optimizations became even more common with the supercomputers in the 1980s. For the current standard, see section 10.1.5.2.4.
2 Once the interpretation of a numeric intrinsic operation is established, the processor may evaluate any mathematically equivalent expression, provided that the integrity of parentheses is not violated.
There is a short table after this paragraph that shows some examples of allowed rearrangements.
So, in the refactored SLSQP I propose to change:
if (diff(one+fac,one)>zero) --> if (fac>=epsilon(zero))
if ( diff(unorm+abs(a(npp1,j))*factor,unorm)>zero ) --> if (abs(a(npp1,j))*factor >= unorm*epsilon(zero))
if ( diff(hmax+factor*h(lmax),hmax)>zero ) --> if (factor*h(lmax) >= hmax*epsilon(zero))
It seems more explicit and less hacky.
Do you need abs(fac) and abs(factor*h(lmax))? If you know they are positive, then it is alright as is. Also, it might be safer to use epsilon(fac) and epsilon(hmax) than epsilon(zero). Some people write code where the zero is some other KIND, or sometimes even a different type, using fortran’s default conversion to fix things up. So if you use the variable as the argument, you are sure to use the right epsilon value.
Hi @jacobwilliams, as a researcher working on optimization with a particular interest in SQP (see, e.g., the very recent thesis of my Ph.D. student T. M. Ragonneau and our COBYQA solver), I am extremely delighted to see your refactorization of SLSQP, an important solver based on SQP.
Out of curiosity, I have a question regarding your refactorization related to the topic under discussion.
I fully agree with the spirit of these modifications. They should be done and they must be done. However, after such changes, the code will behave slightly differently compared with the original F77 implementation due to finite-precision arithmetic, which is usually unharmful.
The question is, with such a difference, how do you verify the faithfulness of your refactorization? Due to the nonlinearity of the iterative procedure, a tiny difference in the middle may lead to a dramatically different sequence of iterates, especially if the optimization problem being solved is nonconvex, where the iterates may converge to different points. It is thus difficult (if even possible) to tell whether the difference in the computed result is caused by the unharmful changes or by bugs hidden somewhere.
Note that I am not concerned by the difference in the computed result, but by the fact that it is hard to tell whether the difference corresponds to bugs or not.
Let me also stress I am not doubting the faithfulness of your refactorization, but asking about your methodology for verifying the faithfulness, which I would like to learn. IMHO, this is a critical question when we refactor old code. Without a systematic, reliable, and automated way of verification, it may be inevitable to introduce bugs provided that the original code is complicated enough.
The same question lies in the very center of PRIMA, my project of modernizing Powell’s optimization solvers. Thus I hope to hear your options and those of everyone reading this.
Thank you.
As I indicated in this thread earlier, it may be helpful to exercise the code for, say, a few months, with some redundancy programmed in, by calculating the iteration termination condition both ways (new and old), and putting in a check for consistency.
Yeah, I thought about that. Probably I could add the abs just to be safe. Probably a couple extra abs aren’t going to slow anything down that much. But since the original code didn’t have that maybe it isn’t necessary? I need to study the code some more…
I think the kind is OK, since in this code, all the reals are the same kind (which can be changed by a preprocessor directive if desired).
I can’t speak to how @jacobwilliams is doing it, but here’s how I’d do it.
Write a series of unit tests of the procedure (or procedures that are intended to work together), to make sure you understand the existing behaviour. Make various breaking changes to the code to see that you have sufficiently covered the behaviour with your test cases (this is sometimes referred to as fuzz testing). Now you can refactor with confidence that if your tests pass, you haven’t broken anything.
There are a few things that can go wrong in this process and must be evaluated on a case by case basis though.
Testing can only prompt these questions, it can’t answer them. A sufficient background knowledge in the domain is required to be able to answer them. Luckily, once answered, they are now documented in the test suite (or at least they are if you’re writing well structured tests 
).
The three most important aspects with any change management of software: test, test, and test.
Totally agree. The question here is how to test.