2

I came across a rather peculiar behavior when was passing arrays to LSODA procedure of ODEPACK library.

The scene: Intel Fortran compiler, one of the recent ones, OS Windows. I have a project that uses features of Fortran 2003, but dirty computational work is passed to LSODA, which is Fortran 77.

LSODA is declared as

  SUBROUTINE DLSODA (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
 1            ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, JT)
  EXTERNAL F, JAC
  INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW, JT
  DOUBLE PRECISION Y, T, TOUT, RTOL, ATOL, RWORK
  DIMENSION NEQ(*), Y(*), RTOL(*), ATOL(*), RWORK(LRW), IWORK(LIW)

note the arrays. To call all that I do:

CALL DLSODA(FDARCY, NEQ, SATCCUR, TCUR, TNEXT, ITOL, RTOL, ATOL, ITASK, &
            ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JDUM, JT)

Let's focus on SATCCUR (in Y position of LSODA). It is declared:

REAL(8), DIMENSION(:), ALLOCATABLE :: SATCCUR

Just in case, in the debugger I checked that DOUBLE PRECISION is REAL(8). At some point SATCCUR is also allocated. The problem is, that LSODA with this call sees complete garbage. The pointer Y sometimes points to something that I cannot either read or write, leading to the program crash. The only way to fix the problem was to change the call to:

CALL DLSODA(FDARCY, NEQ, SATCCUR(1), TCUR, TNEXT, ITOL, RTOL, ATOL, ITASK, &
            ISTATE, IOPT, RWORK(1), LRW, IWORK(1), LIW, JDUM, JT)

Notice that I am not passing the array, but its first element --- the pass by reference of Fortran helps here, as LSODA would effectively get the address of the first element. And I had to do it with all arrays in the call.

If I use the example that is supplied by the LSODA's documentation, everything works fine with passing whole arrays. The differences that I could spot:

  • The arrays in LSODA's example are statically allocated, whereas I do dynamic allocation.
  • In place of function F I pass the function internal to the subroutine. It makes my life much easier, as it captures parameters necessary to run it. This, I believe, is Fortran 2003 feature. But I can't see how this can create the problem as the problem occurs much earlier than F call.

So, the question is why does this happen? Is it an expected behavior?

PS. Also, for LSODA to compile I had to switch off the function interface check as it's internally passing DOUBLE PRECISION array where INTEGER array is expected --- shameless hack on the authors side!

UPDATE

I've managed to trace down the error. The problem was that in my interface declaration for DLSODA I have put

INTERFACE
    SUBROUTINE DLSODA(F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK, &
        ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, JT)
    USE GENERAL
    EXTERNAL :: F, JAC
    INTEGER :: NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, LIW, JT, IWORK
    REAL(DP) :: T, TOUT, RTOL, ATOL, RWORK
    DIMENSION RWORK(LRW), IWORK(LIW)
    REAL(DP), DIMENSION(:) :: Y
    END SUBROUTINE DLSODA
END INTERFACE

Note the declaration for Y and how it is different from RWORK. Apparently, there is a difference between

REAL(8), DIMENSION(:) :: X

and

REAL(8) X
DIMENSION X(*)

I thought this was just a matter of style. Can someone explain why there is a difference?

mobiuseng
  • 2,326
  • 1
  • 16
  • 30
  • @francescalus was writing such a long post, forgot the question! So, is there a reason why this happens? Is it expected that garbage would be passed if I am passing an allocated array? Is it a part of the standard? – mobiuseng Jan 29 '18 at 18:38
  • @francescalus I've edited the question. – mobiuseng Jan 29 '18 at 18:39
  • 1
    There is nothing inherently wrong having an (allocated) allocatable array as an actual argument to an assumed-size dummy argument; there is nothing inherently wrong (in modern Fortran) with having an internal procedure as an actual argument. That is, your problem could well be elsewhere, and memory issues can manifest in strange ways. We probably can'thelp much more unless you can reduce your program to a complete (and minimal) test case (such as shown in [mcve]). – francescalus Jan 29 '18 at 21:35
  • @francescalus OK, I thought I might have been mssing something in Fortran specs regarding passing arrays. I will try to get the problem smaller, so I can post it here. Yet, it's so strange: even with a debugger I can see how garbage appears on LSODA sight, yet on a caller side it's all fine. – mobiuseng Jan 29 '18 at 22:05
  • @francescalus I've found the source of the error, any ideas why there is a difference between `DIMENSION(:)` and `DIMENSION X(*)` declarations? – mobiuseng Jan 30 '18 at 07:42
  • 1
    There's a massive difference between `a(:)` and `a(*)`. Does [this question](https://stackoverflow.com/q/24543329/3157076) cover it for You? – francescalus Jan 30 '18 at 07:50
  • @francescalus Yes, thanks! My understanding is that with `DIMENSION(:)` declaration compiler passes extra information regarding array size, perhaps, appended in the front. So, what reaches the procedure is not the pointer to the first array element, and thus I was getting a garbage. – mobiuseng Jan 30 '18 at 08:00

0 Answers0