How to define dynamic arrays for intent out in Fortran subroutines and call them in Python using f2py?

Question

I am new to using modern Fortran constructs. I have successfully applied dynamic allocation array as input and have a single output as result as verified and confirmed in this jupyter notebook cell:

import numpy as np
import locateindex

tstarr=np.array([1,3,2,4,10])
chk=3
indxval=0

indxval=locateindex.spallocindx(tstarr,chk)
print('indxval: ',indxval)

The cell returns the correct Fortran index location (starting from 1) for an input array. The Fortran subroutine spallocindx.f90 is provided below:

SUBROUTINE spallocindx(indx, xarray, chkval) !RESULT(indx)

  IMPLICIT NONE
  INTEGER, PARAMETER :: INT64=8
  INTEGER(KIND=INT64), INTENT(IN) :: xarray(:)
  INTEGER(KIND=INT64), INTENT(IN) :: chkval
  INTEGER(KIND=INT64), INTENT(OUT) :: indx
  
  indx = findloc(xarray,chkval,dim=1)
 
  END SUBROUTINE spallocindx

This subroutine is successfully compiled for use in python with the numpy.f2py routine as follows:

(venv) MyPythonDirUbuntu$ python -m numpy.f2py -c spallocindx.f90 -m locateindex

So, now I am trying to extend this concept to solving a Lower triangular matrix (essentially a functionality of python's scipy.sparse.linalg.spsolve_triangular).

My key Fortran subroutine is

SUBROUTINE fwdsubv2(xsol,ElIdArr,Lval,rhsvec,nsol) !RESULT(xsol)

  IMPLICIT NONE
  INTEGER, PARAMETER :: REAL_KIND=8
  INTEGER, PARAMETER :: INT64=8

  INTEGER(KIND=INT64), INTENT(IN) :: nsol
  INTEGER(KIND=INT64), INTENT(IN) :: ElIdArr(:)

  REAL(KIND=REAL_KIND), INTENT(IN) :: Lval(:), rhsvec(:)
  
  REAL(KIND=REAL_KIND), ALLOCATABLE, INTENT(OUT) :: xsol(:)

  REAL(KIND=REAL_KIND) :: rowsum, dgval 
  INTEGER(KIND=INT64) :: elidstart, elidend, dgelid, iz, offdgindx, dgindx
  INTEGER(KIND=INT64) :: ir, remraw, colid

  ALLOCATE(xsol(nsol), source=0.0_REAL_REAL_KIND)

  xsol(1)=rhsvec(1)/Lval(1)

  DO ir=2,nsol
    rowsum = 0.0
    dgelid=(ir-1)*nsol+ir
    CALL spallocindx(dgindx, ElIdArr, dgelid)

    IF (dgindx == 0) THEN
        print '(2g0)', 'WARNING: Diagonal value missing for Lower Triangle row ', ir
        dgval=0.000000001
     ELSE
        dgval=Lval(dgindx)
     END IF

    elidstart=(ir-1)*nsol+1
    elidend=(ir-1)*(nsol+1)

    DO iz = elidstart,elidend
        CALL spallocindx(offdgindx, ElIdArr, iz)
        IF (offdgindx > 0) THEN
           remraw=MOD(iz,nsol)
           IF (remraw > 0) THEN
              colid=remraw
           ELSE
              colid=nsol
           END IF
           rowsum = rowsum - Lval(iz) * xsol(colid)/dgval
        END IF
     END DO

     xsol(ir) = rhsvec(ir)/dgval + rowsum

  END DO
  
END SUBROUTINE fwdsubv2

The internal subroutine spalloc.f90 works and has been verified by jupyter cell above. When trying to compile with f2py, using as before:

(venv) MyPythonDirUbuntu$ python -m numpy.f2py -c fwdsubv2.f90 spallocindx.f90 -m tstfwdsub

I get the following error:

fwdsubv2.f90:19:53:

   19 |       ALLOCATE(xsol(nsol), source=0.0_REAL_REAL_KIND)
      |                                                     1
Error: Missing kind-parameter at (1)
error: Command "/usr/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops -I/tmp/tmpp07jlue4/src.linux-x86_64-3.10 -I/home/exouser/venv/lib/python3.10/site-packages/numpy/core/include -I/home/exouser/venv/include -I/usr/include/python3.10 -c -c fwdsubv2.f90 -o /tmp/tmpp07jlue4/fwdsubv2.o" failed with exit status 1

I have tried different variations of the source initialization like source=0.0 or source=REAL_KIND(0.0) or source=REAL_REAL_KIND(0.0) but showing the above example since it seems closest to what is recommended by the Fortran Best Practices linked here.

Looking forward to learning how to do this correctly.

Answer 1

This is purely a problem of incorrect Fortran and it has nothing to do with f2py or Python.

Your kind constant is named REAL_KIND and therefore the literal must look like

0.0_REAL_KIND

There is no reason to duplicate the word REAL_. The constant is called REAL_KIND, not REAL_REAL_KIND.

To be precise, you should also use the kind suffix in

dgval=0.000000001_REAL_KIND

In the other occurence where you set = 0.0 it does not matter that much and you can just use = 0.