Greetings,
I am a member of the SciPy maintainer team and recently I have been translating PROPACK and ARPACK libraries. In the meantime, since it is causing a bit too much of trouble, I started touching slsqp
fortran code located in https://github.com/scipy/scipy/blob/main/scipy/optimize/slsqp/slsqp_optmz.f
Here is the very preliminary code WIP:ENH:optimize:Replace SQP solver fortran code by ilayn · Pull Request #19121 · scipy/scipy · GitHub
I tend to think of myself as knowlegable on the topic due to my PhD but I am quite confused about this piece of code. Mostly because of the flatbuffer practice that is, w(i+1)
is the start of one array but w(j+mq)
is another and so on, that makes everything almost impossible to read. Still, I can find my way out in most cases however I stumbled upon something that is very surprising.
I am sure there is a reason that some machine in year 1964 AD ran out of punch cards or whatever, so they allowed this as valid code but I just want to make sure that this is not a secret gotcha
IF = n*n+1 <----- What is this booby trap here? Is it allowed?
DO 10 i=1,n3
i1 = n1-i
diag = SQRT (l(i2))
w(i3) = ZERO
CALL dcopy_ (i1 , w(i3), 0, w(i3), 1)
CALL dcopy_ (i1-n2, l(i2), 1, w(i3), n)
CALL dscal_sl (i1-n2, diag, w(i3), n)
w(i3) = diag
w(IF-1+i) = (g(i) - ddot_sl (i-1, w(i4), 1, w(IF), 1))/diag <---- IF is a variable now?
Excuse my astonishment but does this mean that
IF = 10
THEN = 20
ENDIF = 30
GOTO = 40
IF IF IF, THEN, ENDIF
C OK the following is just plain stupid but why not
10 GOTO GOTO
20 GOTO GOTO
30 GOTO GOTO
40 ...
should compile?
I checked the modernized version of this code in https://github.com/jacobwilliams/slsqp/blob/master/src/slsqp_core.f90#707 and this assignment is still there so I guess it is still “modern” practice?
While I’m here there is also this type of declarations in some code we already removed so has less of an effect but still equally perplexing for me.
These statements are
What are they trying to achieve here with such code?
Excuse my ignorance if these are common practices but I have no idea how to google these as I tried “cast function with equality”, “Fortran type agnostic variable”, “why Fortran why?” but no avail.