How to fix the errors in compiling ODEPACK in debug mode?

[CRquantum]

Dear all,

I a little curious about the famous odepack Fortran code which contain several robust ODE solvers, which I believe is still heavily used today (it must be wrapped by other Python, Julia packages in some way. ) Code can be downloaded below,
https://people.sc.fsu.edu/~jburkardt/f77_src/odepack/odepack.html

GitHub - jacobwilliams/odepack: Unofficial Mirror of the ODEPACK Library

https://computing.llnl.gov/projects/odepack/software

So I tried to compile it and run it. I use Intel OneAPI + VS2017 in windows 10.
In release mode, it compile and run no problem. But in debug mode, I see error

'The type of the actual argument differs from the type of the dummy argument.

This is a very common issue in some F77 code, which use assumed-size arrays.

So in this case, the program did

      CALL DPREPI (NEQ, Y, S, RWORK(LYH), ...)

Now in the subroutine DPREPI, the argument in RWORK(LYH) position is called YH which is actually an assume-size array,

      SUBROUTINE DPREPI (NEQ, Y, S, YH ...)
      ...
      INTEGER NEQ
      DOUBLE PRECISION Y, S, YH,  ...
      DIMENSION  NEQ(*), Y(*), S(*), YH(*) ...

So I am just confused, in the release mode, did the Fortran code actually load the elements from LYH in RWORK array, i.e., RWORK(LYH:) as YH in subroutine DPREPI?
Or does Fortran just load the single element RWORK(LYH) as the YH in subroutine DPREPI?
I am also confused that, NEQ is an integer when CALL DPREPI, but again in subroutine DPREPI, NEQ is an array actually. Does Fortran automatically understand that when input an interger neq, NEQ(*) will just become an integer? What is the logic?

Overall, my question is,
How to properly solve this kind of error?
It seems not correct to change the definition of NEQ,…, YH in subroutine DPREPI. Do I define some interface of subroutine DPREPI so that it can handle both NEQ, …,YH as number and as array?

Thank you very much in advance!

PS.
I do see that in @jacobwilliams odepack github page, it says it seems allow those strange ‘treating number as array’ behavior.
‘’
ODEPACK contains a few instances where ANSI Fortran 77 is violated: (a) In various places in the LSODES and LSODIS solvers, a call to a subroutine has a subscripted real array as an argument where the subroutine called expects an integer array. Calls of this form occur in Subroutine DLSODES (to DSTODE), in DIPREP (to DPREP), in Subroutine DLSODIS (to DSTODI), and in DIPREPI (to DPREPI). Another such call occurs in the DLSODES demonstration program, from the main program to Subroutine SSOUT. This is done in order to use work space in an efficient manner, as the same space is sometimes used for real work space and sometimes for integer work space. If your compiler does not accept this feature, one possible way to get the desired result is to compile the called routines and calling routines in separate jobs, and then combine the binary modules in an appropriate manner. If this procedure is still not acceptable under your system, it will be necessary to radically alter the structure of the array RWORK within the LSODES or LSODIS solver package. (See also Note 5 below.) (b) Each ODEPACK solver treats the arguments NEQ, Y, RTOL, and ATOL as arrays, even though the length may be only 1. Moreover, except for Y, the usage instructions say that these arguments may be either arrays or scalars. If your system does not allow such a mismatch, then the documentation of these arguments should be changed accordingly.
‘’

But I am just curious, how to set the compiler flag, and how does compiler treat such cases? Like argument is defined as an array, but the actual input is a number.

[sblionel]

The error has nothing to do with assumed-size arrays. It is literally telling you that the types don’t match. I don’t see the declaration of RWORK in your screenshot so I don’t know what it is - but it probably isn’t DOUBLE PRECISION.

In a debug configuration, you get Generated Interface Checking. For procedures that aren’t module procedures, the compiler generates an interface in a special .intmod file for the procedure. If you then call a procedure by that name without an explicit interface visible, it loads the generated module and compares.

In other words, this is an error in the code, and you should fix it.

As for passing a scalar to an array, that’s not allowed in the standard and I would expect an error from that too, depending on compilation order. You ARE allowed to pass an array element to an array, in which case you get “sequence association” (as long as the array you’re passing isn’t assumed-shape or POINTER.) I discuss this in Doctor Fortran in “I’ve Come Here For An Argument” - Doctor Fortran (stevelionel.com)

I’ll mention that these sorts of errors were not caught by older compilers, and many developers depended on compilers not detecting such games.

[themos]

I would add that in addition to not detecting such games-errors, developers got used to compilers implementing this extension-to-the-standard in the same way (always pass the address of the object and always assume that an address has been passed and hope for the best). And that’s how we ended up with compiled executables that changed the value of constants, like 2.

Just write Fortran. And if the hack you want is not doable, pass it to C and try it there.

You are NOT allowed to associate (i.e., pass) arrays to dummy scalars. Cure: index the array so that you associate one element. That element will be potentially modified.

You are NOT allowed to associate scalars to dummy arrays. Cure: find out how big the array must be and construct one, and pass it. Then figure out what the developer intended to happen to the scalar, probably update it with the first element of the returned array.

You ARE allowed to associate an array element with a dummy assumed-size (or explicit shape) array. That is “sequence association”. But you better make sure that there are enough subsequent elements after the array element you passed, as they will potentially be modified.

You ARE allowed to use the newer assumed-rank or assumed-type features of the language.

[CRquantum]

Thank you very much @sblionel the Dr. Fortran and @themos !
Cool! The world is a better place because of you! LOL.

Uhm, well, some of what you said are very advanced, I guess I need some time to fully digest it

Right, as Dr. Fortran @sblionel pointed out, upon checking, indeed,
As show below, on the left hand side, line 9498, I changed the last RWORK(LWM) to int(RWORK(LWM)), because as show at line 9526 on the right hand side, this place is the argument called IWK in the subroutine DPREPI, and IWK is an integer array. That is why there is a ‘type of actual argument …’ error as shown before.
So as I said, changed the last RWORK(LWM) to int(RWORK(LWM)).
My expectation is that, it should serve as IWK. But it did not, it gives another error, as below,

I mean, it looks like, for other arguments, the compiler may take things like RWORK(LWM) as the passing the array RWORK starting from position of LWM.
However, if I do int(RWORK(LWM)), compiler simply pass an integer to IWK, which cause an error bc IWK is defined as array.

However, in release mode, things just fine

Uhm, do you think there are quick fix? Thank you very much!

Happy late New Year!

PS.
On the other hand, I tried to create an array called IWK and passing it to

call DPREPI(..., IWK, ...)

But it just looks very stupid, because it seems I have to manually assign the size of IWK to a fixed number, say IWK(1000). If I want to define

integer  IWK(size(RWORK))

it gives me an error saying that size cannot be calculated. Damn, LOL!

[ivanpribec]

You may also consider

I hit the same issues in 2018. Thankfully, the replies from @sblionel and @mecej4 are still available at the Intel Community forum.

[sblionel]

As I wrote earlier, Fortran does not allow you to pass a scalar to an array. And, again, the reason it works in release mode is that a release configuration doesn’t enable the cross-source-file generated interface check.

If the value is an input, and not written to, you can pass [int(RWORK(LWM))] , making it an array constructor. Any changes made in the called routine will be discarded on return.

[CRquantum]

Thank you very much @sblionel - Dr. Fortran, and @ivanpribec !
The array constructor

[int(RWORK(LWM))] 

or

[int(RWORK(LWM:))] 

is a cool technique!

The odepack code seems indeed reply on these old F77 games.
I also read the post @ivanpribec pointed out, also Dr. Fortran @sblionel have explanation there.
It looks like the most easy way for odepack to run on modern compiler is to select the corresponding files, and disable interface check, like below,

Then everything work fine!
But obviously some of ways odepack was written is not standard Fortran, but it works LOL.

Also thanks @ivanpribec , yes, I have the modern f90 version of dvode,

Has anyone used DVODE and is it good?

The interpolation routine DVINDY has the following requirement, which is stated in the source code: For valid results, T must lie between TCUR - HU and TCUR. In other words, if interpolation is desired, it should be performed after each successful integration step, which covers Tcur - Hu to Tcur (both evaluated with the current value of Tcur). If, instead, you perform a number of successive integration steps and only then try to interpolate over the entire range of the independent variable,…

and you taught me how to use the DVINDY , that is fantastic! In DVODE, there are comments in the code instruct how to use lapack, and I used lapack. It works very good!

Thank you very much guys! Super!

[ivanpribec]

For the record, compiling ODEPACK outside of Visual Studio is straightforward. After downloading the needed files from Netlib, you just need the command

$ ifort -c -nogen-interfaces opksmain.f opksa1.f opksa2.f 

On Linux you would normally proceed to create a static library with

$ ar rcs libopks.a opksmain.o opksa1.o opksa2.o

which can then be linked with -lopks (for the single precision version), assuming the library is somewhere on your linker path.

Instead of ODEPACK you can also consider using SUNDIALS. This library is written in C but also provides a Fortran interface. Usage may seem a bit complex at first, but in a different thread, I think it was @nicholaswogan who mentioned that it works well.

[CRquantum]

Thank you very much @ivanpribec !
I have put the whole odepack in the below repo,

Chen Rong / odepack_test

all the solver are double precision version downloaded LLNL website,
https://computing.llnl.gov/projects/odepack/software
Besides odepack, include also dvode, dvodpk, gearbi, krysi.

All the code are original version.
Only changes are deleted tduplicated lapackor linpack whatever subroutines in dvode, dvodpk, gearbi, krysi which are already included in the two odepack auxiliary files.

Using VS2017 and open the sln file should be able to compile and run without problem both in debug and release mode.

Besides setting ‘check interface’ to NO. For gearbi, also need to in the runtime, set ‘check array bounds’ to NO.

[CRquantum]

@ivanpribec

For the record, compiling ODEPACK outside of Visual Studio is straightforward. After downloading the needed files from Netlib, you just need the command

$ ifort -c -nogen-interfaces opksmain.f opksa1.f opksa2.f 

Thank you very much @ivanpribec !
Hmmm, for gfortran, does

-nogen-interfaces

basically equal to

-std=legacy

is it right?

I checked
https://manpages.ubuntu.com/manpages/impish/man1/alpha-linux-gnu-gfortran-10.1.html

it says,

[sblionel]

For ifort in Visual Studio, simply set the project property Fortran > Diagnostics > Check Routine Interfaces to No. This disables the checks for argument consistency when calling routines without explicit interfaces. It’s similar to gfortran’s -fallow-argument-mismatch - std=legacy is a much bigger hammer.

[CRquantum]

@sblionel

std=legacy is a much bigger hammer.

Thank you very much Dr. Fortran! @sblionel .
Hmmm, may I ask, so how big is this std=legacy hammer?

[sblionel]

@CRquantum

Hmmm, may I ask, so how big is this std=legacy hammer?

According to the gfortran documentation, “The default value for std is ‘gnu’, which specifies a superset of the latest Fortran standard that includes all of the extensions supported by GNU Fortran, although warnings will be given for obsolete extensions not recommended for use in new code. The ‘legacy’ value is equivalent but without the warnings for obsolete extensions, and may be useful for old non-standard programs.” This means it won’t warn you of any extensions or problematic usages.

[urbanjost]

Hopefully, we will get time to remove the issues with the work directories in the ODEPACK on the @jacobwilliams site; it is complicated by it not just being used for scratch work space but being used to return values, and being used as type doubleprecision and integer. Before actually making significant changes we were looking to create a true set of unit tests. If anyone has small programs that can be used as test cases please contribute them. The goal is to make ODEPACK into a single module, and then hopefully use it as a demonstration of co-array usage if time permits, but before making significant changes it requires a complete test suite or development will be much slower.

[mecej4]

Indeed, without having a comprehensive set of tests and verified results, our intended improvements can end up being new bugs!

The maxim, Primum non nocere , is quite apt here.