How to pass arrays with general rank in f2py

Help
1

Hi, (related to stackoverflow), and a bit more.

Suppose I want to use f2py to wrap something from Fortran to Python to call BLAS. In python, I can pass general rank arrays without declaring them in each def. In numpy array, only defines the rank when the first time introduces them.

If I pass python arrays to Fortran, in the Fortran subroutine, I need to declare the ranks for all arrays, including the input and output. The Python input can vary its rank. But, in one subroutine of Fortran, the rank is fixed.

Is there any solution without creating each subroutine for each rank? Assume rank(..) has limited access to the array itself. I can use assume size (*) to accept a general rank array as input. But if I want to utilize dgemm, I need to specify the row and column matrix. Even if I tried to do it with a 1D array by specifying the positions in memory, the rank of the output array, if assume size, is 1. I need to declare the rank according to what is the input in Python. I cannot customize it Fortran in one subroutine.

For example, how to generalize the Fortran code below that can match the input and output ranks from Python without specifying in the particular case (I know this case all are 2, but they can vary as tensor contraction in einsum)

subroutine matmul(A, B, C, n)
    use iso_c_binding
    use types

    implicit none
    integer(kind=c_int), intent(in) :: n
    real(kind=8), intent(in) :: A(n,n), B(n,n)
    real(kind=8), intent(inout) :: C(n,n)
    integer :: i, j, k
    
    do i = 1, n
        do j = 1, n
            C(i, j) = 0.0
            do k = 1, n
                C(i, j) = C(i, j) + A(i, k) * B(k, j)
            end do
        end do
    end do
end subroutine matmul

by
python -m numpy.f2py -c -m fortran_matmul fortran_matmul.f90
and the following python file

import numpy as np
import fortran_matmul

n = 3

A = np.random.rand(n, n)
B = np.random.rand(n, n)
C = np.zeros((n, n))

A = np.asfortranarray(A)
B = np.asfortranarray(B)
C = np.asfortranarray(C)

fortran_matmul.matmul(A, B, C, n)

print(C)