# Help me understand this Fortran code

Hi,

First of all I must point out that I am not a fortran developer. All I need is to understand the fortran code and rewrite it in some other language. So I am here looking for help to better understand what is going on.

If you take a look at SUVM subroutine and focus on the `RPROPS` argument. According to my understanding, the `RPROPS` is a one-dimensional array

``````      DIMENSION
1    IPROPS(*)          ,RPROPS(*)          ,RSTAVA(MSTRE+1)    ,
2    STRAT(MSTRE)       ,STRES(MSTRE)
``````

and the same is suggested also in the documentation:

`RPROPS`. Array of real material properties. It contains the elastic properties (Young’s modulus, E, and Poisson+s ratio, v) and the pairs {x, y} of sampling points along the hardening curve.

RPROPS = [E, v, x_1, y_1, … , x_n, y_n]

With that out of the picture, let’s take a loook at the line 56 of SUVM subroutine

``````SIGMAY=PLFUN(EPBARN,NHARD,RPROPS(IPHARD))
``````

where PLFUN function is called. The way this function is called in the above example makes me think the arguments of PLFUN are all scalars:

• EPBARN is a scalar (double).
• NHARD is a scalar (integer).
• RPROPS(IPHARD) is according to my understanding just an element from an array, thus a scalar. (note that IPHARD is a constant integer within the SUVM subroutine).

## Now the big However

However, the definition of PLFUN is as follows:

``````      DOUBLE PRECISION FUNCTION PLFUN(X, NPOINT, XFX)
C
INTEGER NPOINT, I
DOUBLE PRECISION X, XFX(2,*)
C***********************************************************************
C PIECEWISE LINEAR FUNCTION DEFINED BY A SET OF NPOINT PAIRS
C {X,F(X)} STORED IN THE MATRIX XFX (DIM. 2*NPOINT).
C***********************************************************************
C do something ....
RETURN
END
``````

where, as you can see, the argument is a two-dimensional MATRIX.

## Confused.

So I don’t get it. In the SUVM subroutine, the PLFUN is called with scalars, while the definition of the PLFUN should have a matrix.

2 Likes

@mjancic thank for the question and welcome to the forum! My understanding is that in F77 you can pass the element of an array into a function that expects an array (or matrix) and it will take the array from that element. In modern Fortran you would instead do

``````SIGMAY=PLFUN(EPBARN,NHARD,RPROPS(IPHARD:))
``````

Notice the colon, which passes a subsection of the array starting from IPHARD.

1 Like

Hi,

Ok, that makes some sense, but not really a lot if I try to understand the meaning of the code. In the `SUVM` subroutine the `PLFUN` is called as

``````SIGMAY=PLFUN(EPBARN,NHARD,RPROPS(IPHARD)),
``````

where `IPHARD` is a constant set to 4 as seen from line

``````PARAMETER(IPHARD=4  ,MSTRE=4)
``````

A bit of background: The `PLFUN` is supposed to be a piecewise linear interpolation function. So it should take some x and return the y value for a given set of linearly interpolated points {x,y}.

My concerns:
If the documentation is correct, then

RPROPS = [ E , v , x_1, y_1, … , x_n, y_n]

and if you are correct then

``````RPROPS(IPHARD)
``````

returns an array `[y_1, x_2, y_2, … , x_n, y_n]` thus skipping the first x_1 which to me makes no sense at all.

Is this yet another example of losing a potential Fortran programmer because they see this terrible obsolete FORTRAN 77 code that looks like it was written 30 years ago? (see NEK for computational fluid dynamics moving to C++? - #3 by jacobwilliams)

2 Likes

If RPROPS is indeed as given, then you would expect IPHARD to have the value 3. But are you sure RPROPS is filled in this way? Can you see that in the code?

In Fortran 77, it was not uncommon to pass a 1D array to a subroutine expecting a 2D array (remembering that Fortran is column major), or to pass a single element, which is treated as a pointer to the array starting with that element. I am probably not expressing things precisely. The following code demonstrates how arrays can be passed in the old way:

``````subroutine sub(mat,ncol)
integer, intent(in) :: ncol,mat(2,*)
integer             :: i
print*,"in sub, ncol =",ncol
do i=1,ncol
print*,i,mat(:,i)
end do
print*
end subroutine sub

program main
integer, parameter :: nrows = 2, ncol = 3, n = nrows*ncol
integer            :: i,j,k,mat(n)
forall (i=1:n) mat(i) = i
call sub(mat,ncol)
do j=0,ncol-1
k = 2*j + 1
print*,"calling sub with k =",k
call sub(mat(k),ncol-j)
end do
end program main
``````

It gives output using gfortran or Intel Fortran of

`````` in sub, ncol =           3
1           1           2
2           3           4
3           5           6

calling sub with k =           1
in sub, ncol =           3
1           1           2
2           3           4
3           5           6

calling sub with k =           3
in sub, ncol =           2
1           3           4
2           5           6

calling sub with k =           5
in sub, ncol =           1
1           5           6
``````
2 Likes

Maybe, but there are lots of modern Fortran codes in the Finite Elements and Computational Fluid Dynamics sections of Fortran code on GitHub, so I don’t think the state of the language or tooling is to blame. Some groups will modernize their code and some will not.

This kind of thing is why I recommend never to have explicit shape arrays as dummy arguments.
It either:

• Doesn’t do what you’d expect (i.e. check that the shape of the actual argument matches that of the dummy argument)
• Is trying to do something “clever” that other people won’t understand

Basically, in days of old, it was specified “how” memory was laid out for arrays. They are contiguous, with the lowest index increasing first for subsequent elements in memory.

So a matrix of the form:

``````integer :: A(2,2)
A(1,1) = 1
A(1,2) = 2
A(2,1) = 3
A(2,2) = 4
``````

is laid out in memory like

``````[  1  ] [  3  ] [  2  ] [  4  ]
A(1, 1) A(2, 1) A(1, 2) A(2, 2)
``````

This allowed efficient passing of arrays by using a single pointer. When passing an array with an explicit shape, the actual can be one of either

• The whole array, in which case the dummy argument receives it starting at the beginning
• An element of the array, in which case the dummy argument receives it starting at that point

In this case, it’s using the second option, and so on the calling side you have any array like

``````[ E ] [ v ] [x_1] [y_1] [x_2] [y_2] ...
A(1)  A(2)  A(3)  A(4)  A(5)  A(6) ...
``````

which would then be received as

``````[ x_1 ] [ y_1 ] [ x_2 ] [ y_2 ] ...
A(1, 1) A(2, 1) A(1, 2) A(2, 2) ...
``````

or at least that’s what I’m guessing is intended, but if your starting index is off, then you’ve probably got a bug.

This is the kind of bug I actually find fairly frequently in old codes, because the `x`s and `y`s always happened to be similar enough in their particular use case, meaning nothing would be sufficiently off to notice.

If you change the dummy argument to be assumed shape, this kind of trick doesn’t work anymore, because it has to pass along the shape information as well (i.e. how big is the first dimension, vs the second, etc.) So if you change the dummy argument to `double precision xfx(:,:)` you’ll get a compile time error of the form `array rank mismatch, passed rank(0) to rank(2)`.

3 Likes

Actually you’d need a reshape (assuming you make the dummy argument assumed shape like i suggest above). I.e.

``````SIGMAY=PLFUN(EPBARN,NHARD,reshape(RPROPS(IPHARD:), [2, NHARD]))
``````
1 Like

But I think reshape can only be passed to `intent(in)` argument? If it is `intent(inout)` which I believe was the case above, then that probably would not work.

That’s true, but if you’re using assumed shape, then you have to pass an array of the same rank, so you’ll need an intermediate variable if you want `intent` `out`, or `inout`.

1 Like

Hi @mjancic,
Please note that this is a update procedure for elasto-plastic piecewise-linear material. So, we have to track our EPBARN (strain) - and see whether it crossed elastic limit or in which segment of the defined piece wise linear, it is located during the analysis - after that move along that line.
Note that RPROPS(2) - is youngs modulus
RPROPS(3) - is Poisson ratio
from - RPROP(4) - your function coordinates start - that’s why that is passed. If you notice in PLFUN - XFX is 2 x NPOINTS - that means it’s first column - contains all x values and second column contains all y values.
So, XFX(1,i) - has your x coordinates and XFX(2,i) has your y coordinates.
As said earlier - there is automatic reshaping going on.
[xi,y1,x2,y2,x3,y3,x4,y4] → is being converted to [x1 y1
x2 y2
x3 y3
x4 y4]
As you enter PLFUN it first checks whether your strain(EPBARN) crossed all the x points - then it returns the last y coordinate (XFX(2,NPOINTS)) - NPOINTS is the number of points - so when used as index it gives the last y - value (Reached plastic limit). Then it checks if strain didn’t cross even first value (X1) - if yes - then it returns XFX(2,1) (First y value). If the strain is in between these two limits - it will interpolate - y = y1 + (slope)*(x-x1).
Hope this helps.
By the way in which language are you rewriting ?
The code is good and can be easily modified to Modern Fortran. It references some section from a book (I guess) what is the reference.

I must be missing something really basic here.

Assuming an array

``````RPROPS = [ *E* , *v* , x_1, y_1, … , x_n, y_n]
C index= [  1  ,  2  ,  3  , 4 ,..., n-1 , n ]
``````

So `RPROPS(1)` is Young’s modulus and `RPROPS(2)` is Poisson ratio. That’s according to my understanding. Can you elaborate why `RPROPS(2)` is not the Young’s modulus and not Poisson ratio?

FreeFem++.

The code doesn’t use RPROPS(1) at all - see line 32 and 33

2 Likes

Maybe I am missing something but sequence association is perfectly fine in modern Fortran.

``````
Module interfaces
Interface
Subroutine do(rprops)
Double Precision :: rprops(*)
End Subroutine
Function fun(xfx)
Double Precision :: fun, xfx(2, *)
End Function
End Interface
End Module

program test
use interfaces, only:do
double precision,dimension(9) :: rprops
data rprops/-1.0d0, 3.4d-13, 2.7d-6, 2.0d0, 42.0d0, 4*-42.0d0/
call do(rprops)
end program

subroutine do(rprops)
use interfaces, only:fun
double precision, dimension(*) :: rprops
double precision :: sigmay
integer, parameter :: iphard = 4

sigmay = fun(rprops(iphard))
print *,sigmay
end subroutine do

double precision function fun(xfx)
double precision:: xfx(2,*)
fun = xfx(2,1)
end function fun
``````

The standard says (on Sequence association):

`Sequence association only applies when the dummy argument is an explicit-shape or assumed-size array.`

In this case, it is an assumed-size array, XFX(2,*).

`An actual argument represents an element sequence if it is ... an array element designator ...`

This is our case, RPROPS(IPHARD)

`the element sequence consists of that array element and each element that follows it in array element order.`

The sequence is of length 6: RPROPS(4),RPROPS(5)…,RPROPS(9)

`An actual argument that represents an element sequence and corresponds to a dummy argument that is an array is sequence associated with the dummy argument. The rank and shape of the actual argument need not agree with the rank and shape of the dummy argument,`

they don’t and that’s ok.

`but the number of elements in the dummy argument shall not exceed the number of elements in the element sequence of the actual argument.`

number of elements is not defined for assumed-size dummy so keep reading.

`If the dummy argument is assumed-size, the number of elements in the dummy argument is exactly the number of elements in the element sequence.`

So XFX(2,*) is associated with RPROPS(4:9) and you must not reference XFX with second index larger than 3 (the compiler cannot check for you).

Correct. I’m not suggesting that you can’t do it in modern Fortran, but that I recommend avoiding it in modern Fortran because:

• It’s not clear to future readers of the code (especially one’s new to Fortran) what is going on
• In cases where the ranks don’t match it’s not very intuitive (you have to think real hard about how the indices line up)

Not to mention

the compiler cannot check for you

allowing out-of-bounds access errors to possibly go unnoticed, if not likely to occur in the first place. I’ve seen plenty of code that

• reads off the end of an array because the actual argument was smaller
• looks like it reads off the end of an array (reads past the dummy’s declared size), but didn’t actually because the actual argument was big enough

Finding things like the second one in old code is a bit of a roller coaster ride. “Oh that’s an array argument, even though it looked like they were passing a single element. Ok, it’s declared to have `n` values. Wait, but it actually loops over `n*2` values. Oh, but the actual arguments always have `n*2` values.”

1 Like

People have been abusing early compilers’ inability to check for conforming code and the result can be a total disaster.

4 Likes

Hello @mjancic, nice to find you here.

As @Ashok has written before if you check closely in routine `SUVM` on lines 32 and 33 you find,

``````C Set some material properties
YOUNG=RPROPS(2)
POISS=RPROPS(3)
``````

meaning, that from the fourth element and further will be the (x_i,y_i) samplings points of the hardening curve.

1 Like