0

I have this Fortran code :

module example

    use iso_c_binding
    implicit none
    
contains

    subroutine array_by_ref_modifying(array, nbrows, nbcols, coeff) BIND(C, NAME='array_by_ref_modifying')
        !DEC$ ATTRIBUTES DLLEXPORT :: array_by_ref_modifying
        integer, intent(in) :: nbrows, nbcols, coeff
        integer, intent(inout) :: array(nbrows, nbcols)
        integer :: i, j
        do j = 1, nbcols
            do i = 1, nbrows
                array(i, j) = coeff * array(i, j)
            enddo
        enddo
    end subroutine array_by_ref_modifying

end module example

than I compiled to a TestLib.dll that I call from Python as follows :

redist_path = r"C:\Program Files (x86)\Intel\oneAPI\compiler\2021.3.0\windows\redist\intel64_win\compiler"
dll_full_name = r"C:\DLLS\TestLib.dll"

import os
os.add_dll_directory(redist_path)

import ctypes as ct
import numpy as np

fortlib = ct.CDLL(dll_full_name)

nbrows = 2
nbcols = 6
pnbrows = ct.pointer( ct.c_int(nbrows))  
pnbcols = ct.pointer( ct.c_int(nbcols))
myarr = np.array([[1, 2, 3, 4, 5, 6], [7, 8, 9, 10, 11, 12]], dtype=ct.c_int)

print(myarr)

coeff = 2
pcoeff = ct.pointer( ct.c_int(coeff))

pyfunc_array_by_ref_modifying = getattr(fortlib, "array_by_ref_modifying")
pyfunc_array_by_ref_modifying(np.ctypeslib.as_ctypes(myarr), pnbrows, pnbcols, pcoeff)

print(myarr)

The python code outputs :

[[ 1  2  3  4  5  6]
 [ 7  8  9 10 11 12]]
 [[ 2  4  6  8 10 12]
 [14 16 18 20 22 24]]

as expected. Now, what I expect less is that, if I replace the python script "natural" bit

nbrows = 2
nbcols = 6

(leading in the Fortran to a (1:2,1:6) array) with the "less natural" bit

nbrows = 6
nbcols = 2

(leading in the Fortran to a (1:6,1:2) array, thereby showing that a correct array is passed to Fortran) the python script still outputs

[[ 1  2  3  4  5  6]
 [ 7  8  9 10 11 12]]
 [[ 2  4  6  8 10 12]
 [14 16 18 20 22 24]]

As far as I understand, in the line

pyfunc_array_by_ref_modifying(np.ctypeslib.as_ctypes(myarr), pnbrows, pnbcols, pcoeff)

the np.ctypeslib.as_ctypes(myarr) does some pretreatment/transposition. It this really the case ?

I don't want this pretreatment to be done as it has performance overhead that can be significative in real code with big array dimensions.

(I guess that the "correct" way is with the "natural bit", which does not imply any pretreatment.)

CristiFati
  • 38,250
  • 9
  • 50
  • 87
Olórin
  • 3,367
  • 2
  • 22
  • 42
  • Not sure if relevant: what happens if you set function's *argtypes* and *restype*? Anyway my guess is you have an *Undefined Behavior* here, because the array is *2 X 6* (**established at its creation time**), and by switching dimensions you could run into trouble (considering inner *1D* arrays alignments), thing that didn't happen in your case, *Fortran* did the correct multiplication (on the "transposed" array), and you get the multiplied array. – CristiFati Jul 15 '21 at 16:23
  • I'm sorry, seems that I'm missing something here. What exactly is the problem? What doesn't work as it should? – CristiFati Jul 16 '21 at 14:11
  • The problem is that everything works up to a switch of 2 and 6, while only the case with 2 rows and 6 columns should work. – Olórin Jul 16 '21 at 19:42
  • Could you edit the question adding what you think it should output when switching the dimension order? And explain (in a couple of words) why? – CristiFati Jul 16 '21 at 20:27
  • Everything is literally already explained in the question. Isn't it ? – Olórin Jul 17 '21 at 07:19
  • 1
    Hmm, it isn't clear for me what's wrong. The only thing in the question that indicates that smth would be wrong is "*the python script still outputs*", then combined with one of the previous comment, that it shouldn't work like it does when switching the rows/columns. That's why I was asking what do you think it should output in the 2nd case. – CristiFati Jul 17 '21 at 10:38

1 Answers1

1

As I still don't understand what the problem is, I'll just assume, and try to answer "blindly".

I'll start by answering the punctual question (somewhere at the end):
numpy.ctypeslib.as_ctypes converts a NumPy array (or object, or ...) into a CTypes one. But conversion only happens on the metadata (as it is different for the 2 modules). The array contents (or the pointer, or the memory address where the actual array data is (or begins)) is left unchanged (not copied / modified / altered, ...).

References:

  1. [NumPy.Docs]: C-Types Foreign Function Interface (numpy.ctypeslib)

    1.1. The source code (somewhere at the end): [GitHub]: numpy/numpy - numpy/numpy/ctypeslib.py

  2. [Python.Docs]: ctypes - A foreign function library for Python

So, no transposition is done.
I'm not sure what you mean by "pretreatment" (do all the checks and operations performed by as_ctypes fit there?). Same about "natural bit".

Also note that as_ctypes (1st argument of pyfunc_array_by_ref_modifying) is completely unaware of the rest of them (pnbrows and pnbcols in particular).

Even if not directly impacting (I think it's a matter of "luck"), here's something you might want to check out: [SO]: C function called from Python via ctypes returns incorrect value (@CristiFati's answer).

Going to (what I think is) the real problem (that lead to the question):

  • The array is passed to the Fortran subroutine by reference (its buffer start memory address). I don't have an official source to state this (it was an assumption since the beginning), but it's like in C passing a pointer (I guess it has something to do with C bind from the subroutine declaration)
  • No metadata (rows, columns) is passed (otherwise the next 2 arguments would be useless)

The (2D) array is stored in memory as an 1D one: 1st row, followed by the 2nd one, 3rd, ..., and the last one.

In order to reach element array[i][j] (i = [0, row_count), j = [0, column_count)), the following formula (pointer logic) is used:
array_start_address + array_element_size_in_bytes * (i * column_count + j).

Hoping to clear some of the confusion, here's a small C demo. I replicated (what I thought it was) the Fortran subroutine behavior. I also increased the row count to make things clearer.

main00.c:

#include <stdio.h>
#include <string.h>

#define ROWS 3
#define COLS 6
//#pragma align 4


typedef int Array2D[ROWS][COLS];

/*
void modifyArray(Array2D array, int rows, int cols, int coef)
{
    for (int i = 0; i < rows; ++i)
        for (int j = 0; j < cols; ++j)
            array[i][j] *= coef;
}
//*/

void modifyPtr(int *pArray, int rows, int cols, int coef)
{  // This is the equivalent.
    for (int i = 0; i < rows; ++i)
        for (int j = 0; j < cols; ++j)
            (*(pArray + (i * cols + j))) *= coef;
}

void printArray(Array2D array, int rows, int cols, char *head)
{
    printf("\n%s:\n", head);
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j)
            printf("% 3d  ", array[i][j]);
        printf("\n");
    }
}

void modify(Array2D arr, int rows, int cols, int coef, char *head)
{
    printf("\nRows: %d, Cols: %d, Coef: %d", rows, cols, coef);
    Array2D arr_dup;
    memcpy(arr_dup, arr, sizeof(Array2D));
    modifyPtr(arr_dup, rows, cols, coef);
    printArray(arr_dup, ROWS, COLS, head);
}

int main()
{
    Array2D arr = {  // ROWS X COLS
        { 1, 2, 3, 4, 5, 6 },
        { 7, 8, 9, 10, 11, 12 },
        { 13, 14, 15, 16, 17, 18 },
    };
    char *modifiedText = "Modified";

    //printf("Array size: %d %d %d\n", sizeof(arr), sizeof(arr[0]), sizeof(arr[0][0]));

    printArray(arr, ROWS, COLS, "Original");

    modify(arr, 3, 6, 2, modifiedText);
    printArray(arr, ROWS, COLS, "Original");
    modify(arr, 2, 6, 2, modifiedText);
    modify(arr, 1, 6, 2, modifiedText);
    modify(arr, 3, 5, 2, modifiedText);
    modify(arr, 3, 4, 2, modifiedText);
    modify(arr, 5, 2, 2, modifiedText);
    modify(arr, 7, 1, 2, modifiedText);

    printf("\nDone.\n");
    return 0;
}

Output:

[cfati@CFATI-5510-0:e:\Work\Dev\StackOverflow\q068314707]> sopr.bat
### Set shorter prompt to better fit when pasted in StackOverflow (or other) pages ###

[prompt]> "c:\Install\pc032\Microsoft\VisualStudioCommunity\2019\VC\Auxiliary\Build\vcvarsall.bat" x64 > nul

[prompt]> dir /b
code00.py
example.f90
main00.c

[prompt]>
[prompt]> cl /nologo /MD /W0 main00.c  /link /NOLOGO /OUT:main00_pc064.exe
main00.c

[prompt]> dir /b
code00.py
example.f90
main00.c
main00.obj
main00_pc064.exe

[prompt]>
[prompt]> main00_pc064.exe

Original:
  1    2    3    4    5    6
  7    8    9   10   11   12
 13   14   15   16   17   18

Rows: 3, Cols: 6, Coef: 2
Modified:
  2    4    6    8   10   12
 14   16   18   20   22   24
 26   28   30   32   34   36

Original:
  1    2    3    4    5    6
  7    8    9   10   11   12
 13   14   15   16   17   18

Rows: 2, Cols: 6, Coef: 2
Modified:
  2    4    6    8   10   12
 14   16   18   20   22   24
 13   14   15   16   17   18

Rows: 1, Cols: 6, Coef: 2
Modified:
  2    4    6    8   10   12
  7    8    9   10   11   12
 13   14   15   16   17   18

Rows: 3, Cols: 5, Coef: 2
Modified:
  2    4    6    8   10   12
 14   16   18   20   22   24
 26   28   30   16   17   18

Rows: 3, Cols: 4, Coef: 2
Modified:
  2    4    6    8   10   12
 14   16   18   20   22   24
 13   14   15   16   17   18

Rows: 5, Cols: 2, Coef: 2
Modified:
  2    4    6    8   10   12
 14   16   18   20   11   12
 13   14   15   16   17   18

Rows: 7, Cols: 1, Coef: 2
Modified:
  2    4    6    8   10   12
 14    8    9   10   11   12
 13   14   15   16   17   18
Done.

Back to Fortran (saved your script locally) + Python (created a new one).

code00.py:

#!/usr/bin/env python

import ctypes as cts
import sys

import numpy as np


IntPtr = cts.POINTER(cts.c_int)

DLL_NAME = "./example.{:s}".format("dll" if sys.platform[:3].lower() == "win" else "so")


def modify(arr, rows, cols, coef, modify_func):
    print("\nRows: {:d}, Cols {:d}, Coef: {:d}".format(rows, cols, coef))
    arr_dup = arr.copy()
    arr_ct = np.ctypeslib.as_ctypes(arr_dup)
    rows_ct = cts.c_int(rows)
    cols_ct = cts.c_int(cols)
    coef_ct = cts.c_int(coef)
    modify_func(arr_ct, cts.byref(rows_ct), cts.byref(cols_ct), cts.byref(coef_ct))
    print("Modified array:\n {:}".format(arr_dup))


def main(*argv):
    dll00 = cts.CDLL(DLL_NAME)
    func = getattr(dll00, "array_by_ref_modifying")
    #func.argtypes = (cts.c_void_p, cts.c_void_p, cts.c_void_p, cts.c_void_p)
    func.argtypes = (cts.c_void_p, IntPtr, IntPtr, IntPtr)
    func.restype = None

    arr = np.array([
        [1, 2, 3, 4, 5, 6],
        [7, 8, 9, 10, 11, 12],
        [13, 14, 15, 16, 17, 18],
    ], dtype=cts.c_int)

    print("Original array:\n {:}".format(arr))

    modify(arr, 3, 6 , 2, func)
    modify(arr, 2, 6 , 2, func)
    modify(arr, 1, 6 , 2, func)
    modify(arr, 3, 5 , 2, func)
    modify(arr, 3, 4 , 2, func)
    modify(arr, 5, 2 , 2, func)
    modify(arr, 7, 1 , 2, func)


if __name__ == "__main__":
    print("Python {:s} {:03d}bit on {:s}\n".format(" ".join(elem.strip() for elem in sys.version.split("\n")),
                                                   64 if sys.maxsize > 0x100000000 else 32, sys.platform))
    rc = main(*sys.argv[1:])
    print("\nDone.\n")
    sys.exit(rc)

Output:

[prompt]> "f:\Install\pc032\Intel\OneAPI\Version\compiler\2021.3.0\windows\bin\intel64\ifort.exe" /c example.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.3.0 Build 20210609_000000
Copyright (C) 1985-2021 Intel Corporation.  All rights reserved.


[prompt]> link /NOLOGO /DLL /OUT:example.dll /LIBPATH:"f:\Install\pc032\Intel\OneAPI\Version\compiler\2021.3.0\windows\compiler\lib\intel64_win" example.obj
   Creating library example.lib and object example.exp

[prompt]> dir /b
code00.py
example.dll
example.exp
example.f90
example.lib
example.mod
example.obj
main00.c
main00.obj
main00_pc064.exe

[prompt]>
[prompt]> "e:\Work\Dev\VEnvs\py_pc064_03.08.07_test0\Scripts\python.exe" code00.py
Python 3.8.7 (tags/v3.8.7:6503f05, Dec 21 2020, 17:59:51) [MSC v.1928 64 bit (AMD64)] 064bit on win32

Original array:
 [[ 1  2  3  4  5  6]
 [ 7  8  9 10 11 12]
 [13 14 15 16 17 18]]

Rows: 3, Cols 6, Coef: 2
Modified array:
 [[ 2  4  6  8 10 12]
 [14 16 18 20 22 24]
 [26 28 30 32 34 36]]

Rows: 2, Cols 6, Coef: 2
Modified array:
 [[ 2  4  6  8 10 12]
 [14 16 18 20 22 24]
 [13 14 15 16 17 18]]

Rows: 1, Cols 6, Coef: 2
Modified array:
 [[ 2  4  6  8 10 12]
 [ 7  8  9 10 11 12]
 [13 14 15 16 17 18]]

Rows: 3, Cols 5, Coef: 2
Modified array:
 [[ 2  4  6  8 10 12]
 [14 16 18 20 22 24]
 [26 28 30 16 17 18]]

Rows: 3, Cols 4, Coef: 2
Modified array:
 [[ 2  4  6  8 10 12]
 [14 16 18 20 22 24]
 [13 14 15 16 17 18]]

Rows: 5, Cols 2, Coef: 2
Modified array:
 [[ 2  4  6  8 10 12]
 [14 16 18 20 11 12]
 [13 14 15 16 17 18]]

Rows: 7, Cols 1, Coef: 2
Modified array:
 [[ 2  4  6  8 10 12]
 [14  8  9 10 11 12]
 [13 14 15 16 17 18]]

Done.

Conclusions:

  • The array is modified in the same way in the 2 examples, confirming my assumption

  • No matter what row and column numbers are passed, the array is modified from the beginning, left to right, top to bottom. This might be a bit confusing (when picturing the 2D representation and for example passing a column number that's smaller than the actual one). That's why it works with a dimension greater than the actual one. The only thing to be careful about is not to go beyond the number of elements (row_count * column_count), because that would yield Undefined Behavior (and it might crash)

  • I'm not sure about this one, but mentioning it anyway: there might be some other Undefined Behavior cases, e.g. when each row in the array would be padded by compiler in order to be properly aligned (similar to #pragma pack). Example: a char array with 7 columns (a row would have 7 bytes) might be 8 bytes aligned. Not sure how pointer logic copes with that. But I guess (as Fortran arrays are column major) that's what numpy.ctypeslib.ndpointer's flags argument takes care of

Might also want to check:

CristiFati
  • 38,250
  • 9
  • 50
  • 87