Weird behavior with FFT on array views

[david-hoffman]

This issue was originally posted here: numpy/numpy#11762

After some investigation it seems the problem is with the mkl accelerated libraries. I hope this is the right place to raise the issue.

Reproducing code example:

import numpy as np

np.set_printoptions(precision=2)

# make random data
a = np.random.rand(3, 3) + 1j * np.random.rand(3, 3)

a_view = a[np.newaxis]
# make a new array with one empty axis
a_new = a.reshape(1, 3, 3)

# take ffts  along the last axis of these arrays
b = np.fft.fftn(a, axes=(-1,)).squeeze()
b_view = np.fft.fftn(a_view, axes=(-1,)).squeeze()
b_new = np.fft.fftn(a_new, axes=(-1,)).squeeze()

print("b looks like:")
print(b)
print("b_view looks like:")
print(b_view)
print("b_new looks like:")
print(b_new)

# check that they all are equal
assert np.array_equal(a_new, a_view), "Arrays not equal ..."
assert np.array_equal(b, b_new), "New array failed ..."
assert np.array_equal(b, b_view), "View failed ..."

Error messages:

   ...: 
b looks like:
[[ 2.08+1.42j -0.19+0.05j -0.01+0.6j ]
 [ 2.34+1.96j -0.27-0.16j  0.4 -0.09j]
 [ 2.33+1.79j  0.44-0.62j  0.09+0.38j]]
b_view looks like:
[[ 2.08+1.42j -0.19+0.05j -0.01+0.6j ]
 [ 2.08+1.42j -0.19+0.05j -0.01+0.6j ]
 [ 2.08+1.42j -0.19+0.05j -0.01+0.6j ]]
b_new looks like:
[[ 2.08+1.42j -0.19+0.05j -0.01+0.6j ]
 [ 2.34+1.96j -0.27-0.16j  0.4 -0.09j]
 [ 2.33+1.79j  0.44-0.62j  0.09+0.38j]]
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-1-6af2fa2a5e77> in <module>()
     25 assert np.array_equal(a_new, a_view), "Arrays not equal ..."
     26 assert np.array_equal(b, b_new), "New array failed ..."
---> 27 assert np.array_equal(b, b_view), "View failed ..."

AssertionError: View failed ...

Conda environment

# packages in environment at /Users/david/anaconda3:
#
# Name                    Version                   Build  Channel
_license                  1.1                      py36_1  
alabaster                 0.7.11                     py_3    conda-forge
anaconda                  custom           py36ha4fed55_0  
anaconda-client           1.7.1                      py_0    conda-forge
anaconda-navigator        1.8.7                    py36_0  
anaconda-project          0.8.2                      py_1    conda-forge
appdirs                   1.4.3                      py_1    conda-forge
appnope                   0.1.0            py36hf537a9a_0  
appscript                 1.0.1            py36h1de35cc_1  
asn1crypto                0.24.0                     py_1    conda-forge
astroid                   2.0.4                    py36_0  
astropy                   3.0.4            py36h1de35cc_0  
atomicwrites              1.1.5                    py36_0  
attrs                     18.1.0                     py_1    conda-forge
automat                   0.7.0                    py36_0  
babel                     2.6.0                      py_1    conda-forge
backcall                  0.1.0                    py36_0  
backports                 1.0                      py36_1  
backports.functools_lru_cache 1.5                        py_1    conda-forge
backports.shutil_get_terminal_size 1.0.0                      py_3    conda-forge
beautifulsoup4            4.6.3                    py36_0  
bitarray                  0.8.3            py36h1de35cc_0  
bkcharts                  0.2              py36h073222e_0  
blas                      1.0                         mkl  
blaze                     0.11.3                   py36_0  
bleach                    2.1.3                    py36_0  
blinker                   1.4                        py_1    conda-forge
blosc                     1.14.4               hfc679d8_0    conda-forge
bokeh                     0.13.0                   py36_0    conda-forge
boto                      2.49.0                   py36_0  
boto3                     1.7.70                     py_0    conda-forge
botocore                  1.10.71                    py_0    conda-forge
bottleneck                1.2.1            py36h1d22016_1  
bz2file                   0.98                     py36_0  
bzip2                     1.0.6                h1de35cc_5  
ca-certificates           2018.4.16                     0    conda-forge
cairo                     1.14.12              hc4e6be7_4  
certifi                   2018.8.13                py36_0  
cffi                      1.11.5           py36h5e8e0c9_1    conda-forge
chardet                   3.0.4                    py36_3    conda-forge
clangdev                  6.0.0                 default_0    conda-forge
click                     6.7                        py_1    conda-forge
cloudpickle               0.5.3                    py36_0  
clyent                    1.2.2                      py_1    conda-forge
colorama                  0.3.9                      py_1    conda-forge
colorspacious             1.1.2                     <pip>
conda                     4.5.10                   py36_0  
conda-env                 2.6.0                         1  
constantly                15.1.0           py36h28b3542_0  
contextlib2               0.5.5                      py_2    conda-forge
cryptography              2.3              py36hdbc3d79_0  
cryptography-vectors      2.3                      py36_1    conda-forge
curl                      7.61.0               ha441bb4_0  
cycler                    0.10.0                     py_1    conda-forge
cython                    0.28.5           py36h0a44026_0  
cytoolz                   0.9.0.1          py36h1de35cc_1  
dask                      0.18.2                   py36_0  
dask-core                 0.18.2                   py36_0  
datashape                 0.5.4                    py36_1  
dbus                      1.13.2               h760590f_1  
decorator                 4.3.0                    py36_0  
distributed               1.22.1                   py36_0  
docutils                  0.14             py36hbfde631_0  
entrypoints               0.2.3                    py36_2  
et_xmlfile                1.0.1            py36h1315bdc_0  
expat                     2.2.5                hfc679d8_1    conda-forge
fastcache                 1.0.2            py36h1de35cc_2  
ffmpeg                    4.0                  h01ea3c9_0  
filelock                  3.0.4                      py_1    conda-forge
flask                     1.0.2                    py36_1  
flask-cors                3.0.6                      py_0    conda-forge
fontconfig                2.13.0               h5d5b041_1  
freetype                  2.9.1                hb4e5f40_0  
fribidi                   1.0.4                h1de35cc_0  
gensim                    3.4.0            py36h1de35cc_0  
get_terminal_size         1.0.0                h7520d66_0  
gettext                   0.19.8.1             h15daf44_3  
gevent                    1.3.5            py36h1de35cc_0  
giflib                    5.1.4                h1de35cc_1  
glib                      2.56.1               h35bc53a_0  
glob2                     0.6                      py36_0  
gmp                       6.1.2                hb37e062_1  
gmpy2                     2.0.8            py36h6ef4df4_2  
gnutls                    3.5.19               h2a4e5f8_1    conda-forge
graphite2                 1.3.11               h2098e52_2  
graphviz                  2.40.1               hefbbd9a_2  
greenlet                  0.4.14           py36h1de35cc_0  
h5py                      2.8.0            py36hb794570_1    conda-forge
harfbuzz                  1.8.4                hb8d4a28_0  
hdbscan                   0.8.15           py36h7eb728f_0    conda-forge
hdf5                      1.10.2               hc401514_1    conda-forge
heapdict                  1.0.0                    py36_2  
html5lib                  1.0.1                    py36_0  
hyperlink                 18.0.0                   py36_0  
icu                       58.2                 h4b95b61_1  
idna                      2.7                      py36_2    conda-forge
imageio                   2.3.0                      py_1    conda-forge
imagesize                 1.0.0                      py_1    conda-forge
incremental               17.5.0                   py36_0  
intel-openmp              2018.0.3                      0  
ipykernel                 4.8.2                    py36_0  
ipython                   6.5.0                    py36_0  
ipython_genutils          0.2.0                      py_1    conda-forge
ipywidgets                7.4.0                    py36_0  
isort                     4.3.4                    py36_0  
itsdangerous              0.24                       py_2    conda-forge
jasper                    2.0.14               h636a363_1  
jbig                      2.1                  h4d881f8_0  
jdcal                     1.4                        py_1    conda-forge
jedi                      0.12.1                   py36_0    conda-forge
jinja2                    2.10                       py_1    conda-forge
jmespath                  0.9.3                      py_1    conda-forge
jpeg                      9c                   h470a237_0    conda-forge
jsonschema                2.6.0                    py36_1    conda-forge
jupyter                   1.0.0                    py36_4  
jupyter_client            5.2.3                      py_1    conda-forge
jupyter_console           5.2.0                    py36_1  
jupyter_core              4.4.0                    py36_0  
jupyterlab                0.33.11                  py36_0  
jupyterlab_launcher       0.11.2                   py36_0  
keyring                   13.2.1                   py36_0  
kiwisolver                1.0.1                    py36_1    conda-forge
krb5                      1.16.1               h24a3359_6  
lazy-object-proxy         1.3.1            py36h1de35cc_2  
libcurl                   7.61.0               hf30b1f0_0  
libcxx                    6.0.0                         0    conda-forge
libcxxabi                 4.0.1                hebd6815_0  
libedit                   3.1.20170329         hb402a30_2  
libffi                    3.2.1                h475c297_4  
libgfortran               3.0.1                h93005f0_2  
libiconv                  1.15                 hdd342a3_7  
libidn11                  1.33                 hf837533_0    conda-forge
libopenblas               0.2.20               hdc02c5d_7  
libopus                   1.2.1                h169cedb_0  
libpng                    1.6.34               ha92aebf_1    conda-forge
libprotobuf               3.6.0                hd9629dc_0  
libsodium                 1.0.16               h3efe00b_0  
libssh2                   1.8.0                h322a93b_4  
libtiff                   4.0.9                he6b73bb_1    conda-forge
libvpx                    1.7.0                h378b8a2_0  
libwebp                   0.5.2                         7    conda-forge
libxml2                   2.9.8                h422b904_2    conda-forge
libxslt                   1.1.32               h88dbc4e_1    conda-forge
line_profiler             2.1.2            py36h470a237_1    conda-forge
llvm-meta                 6.0.0                         0    conda-forge
llvmdev                   6.0.0                h137f3e6_4  
llvmlite                  0.24.0           py36hc454e04_0  
locket                    0.2.0                      py_2    conda-forge
lxml                      4.2.4            py36hef8c89e_0  
lzo                       2.10                 h362108e_2  
markupsafe                1.0              py36h1de35cc_1  
matplotlib                2.2.3            py36h54f8f79_0  
mccabe                    0.6.1                      py_1    conda-forge
mistune                   0.8.3            py36h470a237_2    conda-forge
mkl                       2018.0.3                      1  
mkl-service               1.1.2            py36h6b9c3cc_4  
mkl_fft                   1.0.5                    py36_0    conda-forge
mkl_random                1.0.1            py36h5d10147_1  
more-itertools            4.3.0                    py36_0  
mpc                       1.1.0                h6ef4df4_1  
mpfr                      4.0.1                h3018a27_3  
mpmath                    1.0.0                    py36_2  
mrcfile                   1.0.4                     <pip>
msgpack-python            0.5.6            py36h2d50403_2    conda-forge
multipledispatch          0.6.0                    py36_0  
navigator-updater         0.2.1                    py36_0  
nbconvert                 5.3.1                      py_1    conda-forge
nbformat                  4.4.0                      py_1    conda-forge
ncurses                   6.1                  hfc679d8_1    conda-forge
nettle                    3.3                           0    conda-forge
networkx                  2.1                        py_1    conda-forge
nltk                      3.3.0                    py36_0  
nodejs                    9.11.1                        0    conda-forge
nose                      1.3.7                    py36_2  
notebook                  5.6.0                    py36_0  
numba                     0.39.0           py36h6440ff4_0  
numexpr                   2.6.7            py36h4f467ca_0  
numpy                     1.15.0           py36h648b28d_0  
numpy-base                1.15.0           py36h8a80b8c_0  
numpydoc                  0.8.0                      py_1    conda-forge
oauthlib                  2.1.0                      py_0    conda-forge
odo                       0.5.1                      py_1    conda-forge
olefile                   0.45.1                     py_1    conda-forge
openblas                  0.2.20                        8    conda-forge
opencv                    3.4.1            py36h6fd60c2_1  
openh264                  1.7.0                         0    conda-forge
openpyxl                  2.5.5                    py36_0  
openssl                   1.0.2p               h1de35cc_0  
packaging                 17.1                     py36_0  
pandas                    0.23.4           py36h6440ff4_0  
pandoc                    2.2.2                hde52d81_1    conda-forge
pandocfilters             1.4.2                      py_1    conda-forge
pango                     1.42.2               h060686c_0  
parso                     0.3.1                      py_0    conda-forge
partd                     0.3.8                      py_1    conda-forge
path.py                   11.0.1                     py_0    conda-forge
pathlib2                  2.3.2                    py36_0  
patsy                     0.5.0                      py_1    conda-forge
pcre                      8.42                 h378b8a2_0  
pep8                      1.7.1                    py36_0  
pexpect                   4.6.0                    py36_0  
pickleshare               0.7.4            py36hf512f8e_0  
pillow                    5.2.0            py36hb68e598_0  
pip                       18.0                     py36_1    conda-forge
pixman                    0.34.0               hca0a616_3  
pkginfo                   1.4.2                      py_1    conda-forge
plotly                    3.1.0                    py36_0  
pluggy                    0.7.1                    py36_0  
ply                       3.11                       py_1    conda-forge
prometheus_client         0.3.1                    py36_0  
prompt_toolkit            1.0.15           py36haeda067_0  
psutil                    5.4.6            py36h1de35cc_0  
ptyprocess                0.6.0                    py36_0  
py                        1.5.4                    py36_0  
pyasn1                    0.4.4                    py36_0  
pyasn1-modules            0.2.2                    py36_0  
pycodestyle               2.4.0                      py_1    conda-forge
pycosat                   0.6.3            py36h470a237_1    conda-forge
pycparser                 2.18                       py_1    conda-forge
pycpd                     0.4                       <pip>
pycrypto                  2.6.1            py36h1de35cc_9  
pycurl                    7.43.0.2         py36hdbc3d79_0  
pyflakes                  2.0.0                    py36_0  
pygments                  2.2.0                      py_1    conda-forge
pyhamcrest                1.9.0                      py_2    conda-forge
pyjwt                     1.6.4                      py_0    conda-forge
pylint                    2.1.1                    py36_0  
pyodbc                    4.0.23           py36hfc679d8_1    conda-forge
pyopenssl                 18.0.0                   py36_0  
pyparsing                 2.2.0                      py_1    conda-forge
pyqt                      5.9.2            py36h655552a_0  
pysocks                   1.6.8                    py36_1    conda-forge
pytables                  3.4.4            py36h247b57e_1    conda-forge
pytest                    3.7.1                    py36_0  
pytest-arraydiff          0.2              py36h39e3cac_0  
pytest-astropy            0.4.0                    py36_0  
pytest-doctestplus        0.1.3                      py_0    conda-forge
pytest-openfiles          0.3.0                      py_0    conda-forge
pytest-remotedata         0.3.0                    py36_0  
python                    3.6.6                h5001a0f_0    conda-forge
python-crfsuite           0.9.6            py36h470a237_0    conda-forge
python-dateutil           2.7.3                      py_0    conda-forge
python-graphviz           0.8.4                    py36_2    conda-forge
python.app                2                        py36_8  
pytz                      2018.5                   py36_0  
pywavelets                0.5.2            py36h7eb728f_2    conda-forge
pyyaml                    3.13             py36h1de35cc_0  
pyzmq                     17.1.2           py36h1de35cc_0  
qt                        5.9.6                h74ce4d9_0  
qtawesome                 0.4.4              pyh8a2030e_1    conda-forge
qtconsole                 4.3.1            py36hd96c0ff_0  
qtpy                      1.4.2              pyh8a2030e_1    conda-forge
readline                  7.0                  hc1231fa_4  
requests                  2.19.1                   py36_1    conda-forge
requests-oauthlib         1.0.0                      py_1    conda-forge
retrying                  1.3.3                    py36_2  
rope                      0.11.0                   py36_0  
ruamel_yaml               0.15.46          py36h1de35cc_0  
s3transfer                0.1.13                   py36_0  
scikit-image              0.14.0           py36h0a44026_1  
scikit-learn              0.19.1           py36hf9f1f73_0  
scipy                     1.1.0            py36hf1f7d93_0  
seaborn                   0.9.0                    py36_0  
send2trash                1.5.0                    py36_0  
service_identity          17.0.0           py36h28b3542_0  
setuptools                40.0.0                   py36_1    conda-forge
simplegeneric             0.8.1                    py36_2  
simpleitk                 1.1.0            py36h0a44026_0    simpleitk
singledispatch            3.4.0.3          py36hf20db9d_0  
sip                       4.19.8           py36h0a44026_0  
six                       1.11.0                   py36_1  
smart_open                1.6.0                      py_1    conda-forge
snowballstemmer           1.2.1                      py_1    conda-forge
sortedcollections         1.0.1                      py_1    conda-forge
sortedcontainers          2.0.4                      py_1    conda-forge
sphinx                    1.7.6                    py36_0  
sphinxcontrib             1.0                      py36_1  
sphinxcontrib-websupport  1.1.0                    py36_1  
spyder                    3.3.1                    py36_1  
spyder-kernels            0.2.6                    py36_0  
sqlalchemy                1.2.10           py36h470a237_1    conda-forge
sqlite                    3.24.0               h2f33b56_0    conda-forge
statsmodels               0.9.0            py36h1d22016_0  
sympy                     1.2                      py36_0  
tblib                     1.3.2                      py_1    conda-forge
terminado                 0.8.1                    py36_1  
testpath                  0.3.1            py36h625a49b_0  
tk                        8.6.8                         0    conda-forge
toolz                     0.9.0                    py36_0  
tornado                   5.1              py36h470a237_1    conda-forge
tqdm                      4.24.0                   py36_0  
traitlets                 4.3.2            py36h65bd3ce_0  
twisted                   17.5.0                   py36_0  
twython                   3.7.0                    py36_0  
typed-ast                 1.1.0            py36h1de35cc_0  
typing                    3.6.4                    py36_0  
unicodecsv                0.14.1                     py_1    conda-forge
unixodbc                  2.3.6                h3efe00b_0  
urllib3                   1.23                     py36_1    conda-forge
viscm                     0.7                       <pip>
wcwidth                   0.1.7                      py_1    conda-forge
webencodings              0.5.1                    py36_1  
werkzeug                  0.14.1                   py36_0  
wheel                     0.31.1                   py36_1    conda-forge
widgetsnbextension        3.4.0                    py36_0  
wrapt                     1.10.11          py36h1de35cc_2  
x264                      1!152.20180717       h470a237_0    conda-forge
xlrd                      1.1.0                      py_2    conda-forge
xlsxwriter                1.0.5                    py36_0  
xlwings                   0.11.8                   py36_0  
xlwt                      1.3.0                      py_1    conda-forge
xz                        5.2.4                h1de35cc_4  
yaml                      0.1.7                hc338f04_2  
zeromq                    4.2.5                hfc679d8_5    conda-forge
zict                      0.1.3                    py36_0  
zlib                      1.2.11               h470a237_3    conda-forge
zope.interface            4.5.0            py36h470a237_0    conda-forge

[oleksandr-pavlyk]

Thank you @david-hoffman. A simpler reproducer is

import numpy as np
import mkl_fft

# make random data
a = np.random.rand(3, 3) + 1j * np.random.rand(3, 3)

# make a view with one empty axis
a_view = a[np.newaxis]
a_view_copy = a_view.copy()

# take ffts  along the last axis of these arrays
b_view2 = mkl_fft.fft(a_view_copy, axis=-1)
b_view = mkl_fft.fft(a_view, axis=-1)

# check that they all are equal
assert np.array_equal(a_view, a_view_copy), "Input arrays not equal ..."
assert np.array_equal(b_view2, b_view), "Result arrays not equal ..."

The trips the code is the zero stride in a_view:

In [1]: import numpy as np

In [2]: a = np.random.rand(3, 3) + 1j * np.random.rand(3, 3)

In [3]: a.reshape((1,3,3)).strides
Out[3]: (144, 48, 16)

In [4]: a[np.newaxis].strides
Out[4]: (0, 48, 16)

A fix shouldn't be too hard.

[seberg]

Just as a side note. That 0 stride could be any arbitrary number. At the same time the array will still be considered C-Contiguous by numpy (in case that has anything to do with the code here).

[oleksandr-pavlyk]

@seberg Thank you for that information. Is there a function in NumPy that makes the strides array canonical, by populating ambiguous values (zero or otherwise) with values that would apply for a generic array evaluated at the shape of incoming array (even if it has unit strides)?

[seberg]

Converting to memoryview will do that. Otherwise, the argument I think is:

• You use the fact that it is contiguous: no need to use strides, just use itemsize instead. Since itemsize + shape give you all information you need.
• You use the strides (support uncontiguous arrays): In that case you will always end up using them together with a shape, so that a funny value does not matter.

We used to make this relaxation with C-contiguous definition because it was causing other weird things and also fast paths being missed, and the only problems I seriously remember were funny mixtures of the two approaches (Such as getting the itemsize from the last dimension after checking contiguouity).

EDIT: I am very sure about the memoryview thing, but I am frankly not sure we have tests for it.

[seberg]

@oleksandr-pavlyk the code here says /* we must have C- or F- contiguous layout */ (not sure where it is checked), so you can definetely use the itemsize!

[oleksandr-pavlyk]

@seberg I am not certain what you mean by memoryview thing. Is that the memoryview in Cython?
I found that ar.reshape(ar.shape).strides gives the canonical version of ar.strides, but I was hoping there is a C function that would do that based on array of strides, shapes and itemsize.

The code checks C- or F- contiguous layout with int any_contig = (PyArray_ISONESEGMENT(x)) ? 1 : 0;.

You are right, itemsize and shape is sufficient for contiguous arrays, but I am hoping to not restrict myself to contiguous arrays one.

The problem I am trying to solve is to detect if a lower-dimensional FFT along a set of axes can be viewed as a sequence of FFT on multiple equidistant datasets. If so, one can make a single call to MKL to perform them, and MKL would be able to do it more efficiently due to threading and vectorization.

This equidistant property is always true for 1D fft on rows or columns of a matrix. For higher array rank tensors, it is true for C- or F- contiguous arrays, if transformation is along the first or the last axis.

But it may be true in a more general settings, for example for `fft( ar[::4, ::4, ::4], axis=-1)``, where the input array is not contiguous.

So in hope to generalize the check in the future, I'd rather avoid making the contiguity assumption.

[seberg]

memoryview is not a cython thing as such (they just have something they call that as well), what I meant is the Python buffer protocol more generally, but it all adds up the same, since cython will use that Python protocol to define its memoryviews from numpy arrays probably.

Ok, so if you want to have arr[::4, ::4] or so also work (in some sense a "reshapable" array), things are more complicating of course, you are right. Otherwise I think the code should be much cleaner here by just checking the contiguity (in some sense the whole point of allowing funny strides is to allow you/us this optimization).
It could be nice to add a function reshape_if_possible in numpy or so, but basically it does not exist. However, I do not understand your min_nonzero here. Unless you know the shape, this is not something you can really do (or that helps you much)? And that function is not using the shape at all.

Did not look at it too long though, so maybe I am missing things.

[seberg]

@oleksandr-pavlyk, I think the general case can probably be achieved using NpyIter and removing the dimension the FFT is created over. NpyIter has all the logic to try to give you the largest possible contiguous block for each iteration. But it would also add some overhead to construct it (meaning that fast paths for the contiguous cases may still be useful).

[seberg]

As another point about the general optimization (frankly, one question is whether it is really worth it for anything but extremely short FFTs), is that you can also imagine something like arr.transpose(1, 0, 2) making things even harder.

So overall, I am getting more convinced that you should just trust NumPy here, and limit (this part of the code) yourself to asking numpy whether the array is C or F-contiguous:

/* for the FFT stride, it is OK if it is 0 or weird, since in that case the dimension must be 1 */
fft_stride = PyArray_STRIDES(x)[axis];
if ((axis == 0) != PyArray_C_CONTIGUOUS(x)) {
     /* FFT-axis is the fastest one */
     outer_stride = PyArray_DIM(x, axis) * PyArray_ITEMSIZE(x); 
}
else {
     /* FFT-axis is the slowest one */
     outer_stride = PyArray_ITEMISZE(x);
}

EDIT: Just to point out, there are probably lots of silly wrong functions names in there :). Also removed a wrong sentence, since I did that simplification already.

EDIT2: A warning, I am not sure the flag returning will return actual "boolean" values, so the != may not be strictly right.

EDIT3: OK, one more edit much later :), the fft stride is always OK...

[oleksandr-pavlyk]

@seberg Thank you for your input and pointing out that strides associated with unit shape elements may be arbitrary.

I have fixed the issue in #22, by providing a function to compute distance between equidistant datasets, that accounts for possibility that stride may not be meaningful for unit shape directions.

[seberg]

Quite frankly, glancing at the new function, I think it is a bit misguided (or solving the wrong problem). The code I gave above (or something similar) solves the issue for contiguous cases.

I guess you attempt to solve it for a more generic case then the contiguous one, but since you do not actually do all the other checks that would be necessary in the generic, non-contiguous case, there is no point.

Now that I think about it, i think the reason for this is that you misinterpreted the ISONESEGMENT macro. That macro does not give you the information that it might be, for example, a transposed array. The result of the macro is (at least in two dimensions) guaranteed to be either C- or F-Contiguous.

However, also your code assumes this, because for the code to work correctly, the FFT axis has to be either the fastest or the slowest one. A general transpose would break this.

Long story short, I think the logic, while sound, is over engineered. If someone needs something like it at some point for a much more general try, you can do it then!

[oleksandr-pavlyk]

@seberg I think there may be a miscommunication here.

Details of implementation aside, given an array, to compute FFT along specified set of axes, one needs to iterate over complimentary axes, calling FFT for every possible combination of complementary indexes.

If distances between starting pointers of each dataset are [dist, 2*dist, ..., p*dist], one can replace p calls to MKL FFT with a single call, see the relevant section of mkl-reference.

Having a single call has many performance advantages. MKL can perform transform computations over independent datasets in parallel, better utilizing compute resources than a serial calls to FFT would.