# Looking at some old code **URL:** \ **Category:** Uncategorized\ **Created:** [October 20, 2022, 8:38pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580 "2022-10-20T20:38:46Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [October 20, 2022, 8:38pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/1 "2022-10-20T20:38:46Z")
Consider [this code](https://github.com/scipy/scipy/blob/main/scipy/optimize/slsqp/slsqp_optmz.f). A statement function: ```fortran diff(u,v)= u-v ``` and an if statement: ```fortran IF(diff(hmax+factor*h(lmax),hmax).GT.ZERO) ``` So, it seems what is happening is a test for: ```fortran hmax + factor*h(lmax) - hmax > zero ``` which reduces to: ```fortran factor*h(lmax) > zero ``` I suspect there is some deep numerical reason for doing it the old way? But what is it? ---
**Author:** ![oscardssmith](https://avatars.discourse-cdn.com/v4/letter/o/b9e5f3/32.png) [@oscardssmith](https://fortran-lang.discourse.group/u/oscardssmith)\ **Post date:** [October 20, 2022, 9:04pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/2 "2022-10-20T21:04:32Z")
`hmax + factor*h(lmax) - hmax >= zero` tests whether `factor*h(lmax)` is less than `eps(hmax)` ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [October 20, 2022, 9:07pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/3 "2022-10-20T21:07:58Z")
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: ```fortran if (factor*h(lmax) < epsilon(hmax)) ... ``` ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [October 20, 2022, 9:08pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/4 "2022-10-20T21:08:19Z")
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. ---
**Author:** ![certik](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/certik/32/4_2.png) [@certik](https://fortran-lang.discourse.group/u/certik)\ **Post date:** [October 21, 2022, 3:41am UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/5 "2022-10-21T03:41:53Z")
Timely, we are almost done implementing statement functions in LFortran: [Added support for ```statement functions``` by Pranavchiku · Pull Request #911 · lfortran/lfortran · GitHub](https://github.com/lfortran/lfortran/pull/911). ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [October 21, 2022, 7:33am UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/6 "2022-10-21T07:33:17Z")
The line with the IF statement also occurs in the [TOMS 587 code](http://www.netlib.no/netlib/toms/587) 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. ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [October 21, 2022, 7:36am UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/7 "2022-10-21T07:36:27Z")
> [@oscardssmith](#): > > `hmax + factor*h(lmax) - hmax >= zero` tests whether `factor*h(lmax)` is less than `eps(hmax)` More precisely, I think, whether `factor*h(lmax) is >= hmax * epsilon(hmax)`, assuming that none of the quantities are negative. ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [October 21, 2022, 1:57pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/8 "2022-10-21T13:57:02Z")
I wonder if any current compiler actually can detect this y+x-y trick and optimize it away (defeating the original intent). ---
**Author:** ![oscardssmith](https://avatars.discourse-cdn.com/v4/letter/o/b9e5f3/32.png) [@oscardssmith](https://fortran-lang.discourse.group/u/oscardssmith)\ **Post date:** [October 21, 2022, 2:22pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/9 "2022-10-21T14:22:41Z")
hopefully it doesn’t without a fastmath flag. ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [October 21, 2022, 2:37pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/10 "2022-10-21T14:37:32Z")
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: ```auto 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: ```auto #define diff(u,v) ((u)-(v)) ``` Again, no error. ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [October 21, 2022, 3:18pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/11 "2022-10-21T15:18:03Z")
> [@mecej4](#): > > I then replaced the DIFF statement function with a preprocessor replacement for the ASF: > > ```auto > #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. ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [October 21, 2022, 3:45pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/12 "2022-10-21T15:45:57Z")
> [@mecej4](#): > > Not yet, I think. That would require that Fortran compilers become capable of algebraic reasoning. 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. ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [October 21, 2022, 11:19pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/13 "2022-10-21T23:19:24Z")
So, in the [refactored SLSQP](https://github.com/jacobwilliams/slsqp) I propose to change: ```fortran if (diff(one+fac,one)>zero) --> if (fac>=epsilon(zero)) ``` ```fortran if ( diff(unorm+abs(a(npp1,j))*factor,unorm)>zero ) --> if (abs(a(npp1,j))*factor >= unorm*epsilon(zero)) ``` ```fortran if ( diff(hmax+factor*h(lmax),hmax)>zero ) --> if (factor*h(lmax) >= hmax*epsilon(zero)) ``` It seems more explicit and less hacky. ---
**Author:** ![RonShepard](https://avatars.discourse-cdn.com/v4/letter/r/a3d4f5/32.png) [@RonShepard](https://fortran-lang.discourse.group/u/RonShepard)\ **Post date:** [October 21, 2022, 11:38pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/14 "2022-10-21T23:38:10Z")
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. ---
**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\ **Post date:** [October 22, 2022, 1:31am UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/15 "2022-10-22T01:31:38Z")
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](https://tomragonneau.com/documents/thesis.pdf) and our [COBYQA solver](https://www.cobyqa.com)), 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. > [@jacobwilliams](#): > > So, in the [refactored SLSQP](https://github.com/jacobwilliams/slsqp) I propose to change: > > ```auto > if (diff(one+fac,one)>zero) --> if (fac>=epsilon(zero)) > > ``` > > ```auto > if ( diff(unorm+abs(a(npp1,j))*factor,unorm)>zero ) --> if (abs(a(npp1,j))*factor >= unorm*epsilon(zero)) > > ``` > > ```auto > if ( diff(hmax+factor*h(lmax),hmax)>zero ) --> if (factor*h(lmax) >= hmax*epsilon(zero)) > > ``` > > It seems more explicit and less hacky. 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](https://github.com/equipez/PRIMA). Thus I hope to hear your options and those of everyone reading this. Thank you. ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [October 22, 2022, 1:55am UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/16 "2022-10-22T01:55:49Z")
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. > [@mecej4](#): > > ```auto > logical :: lxpr1, lxpr2 > ... > lxpr1 = diff(hmax + factor*h(lmax), hmax) .gt. zero ! old way > lxpr2 = factor*h(lmax) > hmax*epsilon(0.0d0) ! new way > if(lxpr1.neqv.lxpr2)print *,'DIAG 3: ',factor,h(lmax),hmax,lxpr1,lxpr2 ! do they disagree? > > ``` ---
**Author:** ![jacobwilliams](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/jacobwilliams/32/10_2.png) [@jacobwilliams](https://fortran-lang.discourse.group/u/jacobwilliams)\ **Post date:** [October 22, 2022, 3:01pm UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/17 "2022-10-22T15:01:35Z")
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). ---
**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\ **Post date:** [October 23, 2022, 12:41am UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/18 "2022-10-23T00:41:21Z")
> [@zaikunzhang](#): > > how do you verify the faithfulness of your refactorization? 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. - The existing code fails one of your new unit tests. Did you misunderstand how the code is intended to behave/be used or did you find a bug, or maybe an undocumented assumption about the inputs? - No tests fail if you make (what you expect to be) a breaking change. Do you need more tests, or is that bit of code unnecessary? - Your refactoring changing the answers, but only by a small amount. What tolerance of change is acceptable? 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 😉). ---
**Author:** ![FortranFan](https://avatars.discourse-cdn.com/v4/letter/f/96bed5/32.png) [@FortranFan](https://fortran-lang.discourse.group/u/FortranFan)\ **Post date:** [October 23, 2022, 1:37am UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/19 "2022-10-23T01:37:00Z")
The **three** most important aspects with any change management of software: **test, test, and test.** ---
**Author:** ![zaikunzhang](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/zaikunzhang/32/785_2.png) [@zaikunzhang](https://fortran-lang.discourse.group/u/zaikunzhang)\ **Post date:** [October 23, 2022, 1:43am UTC](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580/20 "2022-10-23T01:43:37Z")
> [@FortranFan](#): > > The **three** most important aspects with any change management of software: **test, test, and test.** Totally agree. The question here is how to test. [Next page](https://fortran-lang.discourse.group/t/looking-at-some-old-code/4580.md?page=2)