Accepting F-contiguous arrays and passing them to C is questionable

The following code shows that F-contiguous arrays can cause to a user big pain. It is questionable whether they should be accepted by petsc4py.

Passing F-contiguous array to Vec.__getitem__ goes seemsingly fine, the returned array has the same shape. One would intuitively expect the failing assert to pass.

from petsc4py import PETSc
import numpy as np

# Create Vec = 0, 1, 2, ...
v = PETSc.Vec().create(PETSc.COMM_SELF)
v.setSizes(8)
v.setType(PETSc.Vec.Type.SEQ)
idx = np.arange(v.size, dtype=PETSc.IntType)
val = np.arange(v.size, dtype=PETSc.ScalarType)
v.setValues(idx, val)

# All arrays in C-ordering
I = np.arange(6, dtype=PETSc.IntType).reshape(2, 3)
V = v[I]
assert np.allclose(I, V)  # Ok

# Index by F-ordered array
I = I.T
V = v[I]
assert I.shape, V.shape  # Ok
assert np.allclose(I, V)  # Fails because of mixed orderings!

# This is what petsc4py client has to do now
V = v[np.ascontiguousarray(I)]
assert np.allclose(I, V)  # Ok
assert np.allclose(np.ascontiguousarray(I), V)  # Ok

Comments (10)