Help me understand this Fortran code

Help
1

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
2

@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
3

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.

4

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
5

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?

6

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
7

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.

8

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 xs and ys 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
9

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
10

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.

11

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
12

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.

13

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++. :slight_smile:

14

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

2 Likes
15

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).

16

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

  • Newer features (like assumed shape) have made it mostly unnecessary
  • 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
17

RoyBattyunnamed

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

4 Likes
18

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