Home > Mobile >  Is this a correct python port of a fortran function call
Is this a correct python port of a fortran function call

Time:08-09

I'm trying to port fortran code to python (call me crazy), and am wondering whether I'm handling the input of a variable to a function in a correct way.

note that I've read:

Fortran77: what does the asterisk (*) means in subroutine argument?

How can a scalar be passed to a vector (1D array) to a Fortran subroutine?

What is * in FORTRAN 77

, but my question stands.

Focusing on the variable wrk(ia) used in curfit.f to input to subroutine fpcurf.f. wrk is a 1-d array.

curfit.f source

 c  ..array arguments..
      real*8 ... wrk(lwrk) ...
 c  ... no values assigned to wrk here, it is still empty initialized ...
 
 call fpcurf(...,
 * wrk(ia),...)
 ...

(from the links I understand the * is just a line continuation.

curfit.py

...
n, c, fp, ier = call_fpcurf(..., wrk[ia], ...)
...

fpcurf.f takes the variable wrk(ia) as variable a, which is defined as an array. Now I understand that this somehow becomes a "dummy" array (whatever that means).

fpcurf.f source

  recursive subroutine fpcurf(...,
 *   ...,a,...)
  implicit none
c  ..
c  ..scalar arguments..
      ...
c  ..array arguments..
      real*8 ...,
     * ...,a(nest,k1),...

Since wrk is empty initialized in curfit.f, and the element wrk(ia) is inputted to fpcurf.f as a, which has array dimensions (nest, k1), I basically ignore the inputted value of a in fpcurf.py and simply initialize a new empty / zeros array.

fpcurf.py

def fpcurfpy(... a, ...):
    ...
    a = np.zeros(shape=(nest, k1), dtype=np.float64)
    # alt.
    # a = np.empty(shape=(nest, k1), dtype=np.float64)
    ...

Is this python code a correct representation of the fortran code? Assume a is not returned from fpcurf.f , i.e. not used in curfit.f

More specifically: is a = np.zeros(shape=(nest, k1), dtype=np.float64) in my python file the best representation of what's going on in the fortran code, especially how a is treated in fpcurf.f

CodePudding user response:

In the previous century (except its last decade) it was common to re-use one large array for various local working spaces that various subroutines needed. One just created one array and passed it around so that the working space was shared and memory was saved.

In the current century you can use local working spaces and allocate it dynamically or have some object that contains various data structures that the algorithms needs.

Therefore, there is no need to pass around a. It should be created locally or in some class that represents the library or algorithm.


If you do pass a global working array around, do not redefine it locally, it will just redefine the reference and ignore what it received from above:

 def fun(a): 
     a = 5 
     print(a)

 a = []

 fun(a)                                                                                                                                                                                                                                                                                                          
5

In the original code the element wrk(ia) is passed. That is a legal way how to pass a part of the array that begins at element wrk(ia). In modern Fortran we would more likely pass wrk(ia:) instead. It is not clear why this is done. Perhaps to separate the section that is used in fpcurf from other sections used elsewhere.

The function interprets the passed array as a two-dimensional one a(nest, k1). This is again legal. Both of these actions are allowed by the sequence association rules which allow re-interpreting sequences of array elements. The receiving functions just see a part of the wrk array as a 2D array that is used for some internal work.

This way of passing around one large working array is archaic and should not be transferred to new code.

  • Related