3 ms·

What is the advantage of doing this compared to writing Fortran (90+) kernels, use F2Py bindings and combining them with NumPy? For HPC you generally need to be

What is the advantage of doing this compared to writing Fortran (90+) kernels, use F2Py bindings and combining them with NumPy? For HPC you generally need to be able to have full control over compiler flags. Plus, Fortran already has way better multi-dim array handling built-in - NumPy is basically just a wrapper around that.

greglindahl 10y ago
I'm a huge Fortran fan, but even I can see that reducing the number of languages required to get high performance can be a good thing. Not all of my science colleagues are interested in learning a ton of languages.
- m_mueller 10y ago
  Well, most of the scientific computing in Python is happening on NumPy. Which uses Fortran compatible data structures. I doubt that a C++ backed python extension will be compatible with that. For that reason, if anything, I'd advocate to transpile to Fortran rather than C++ if you want more comfort for scientific users.
  - shoyer 10y ago
    NumPy natively supports both row-major (C/C++) and column-major (Fortran) arrays, so from the perspective of NumPy inter-operability there would be no advantage either way.
    - m_mueller 10y ago
      I didn't know that. How does NumPy use its Fortran kernels with row-major arrays? Reordering would have quite significant performance implications.
      - shoyer 10y ago
        NumPy has an abstract iterator interface so you don't need to write a kernel with a specific memory layout in mind: http://docs.scipy.org/doc/numpy-1.11.0/reference/c-api.iterator.html#simple-iteration-example http://docs.scipy.org/doc/numpy-1.11.0/reference/c-api.itera... For operations like adding two arrays or summing an array along an axis, this means that NumPy automatically iterates in the fastest order.
        m_mueller 10y ago
        Interesting. How does it decide the storage order when initializing an array? I'd assume due to python's dynamic nature it can't optimize for later operations there, right?
        infinite8s 10y ago
        It defaults to row major but you can specify column major on creation.
      - dagss 10y ago
        NumPy certainly does not have anything like "Fortran kernels". Elementwise arithmetic is coded in C (and not very fast...). Linear algebra uses the BLAS and LAPACK interfaces, of these only LAPACK is usually written in Fortran these days... but anyway, a Fortran matrix is a C transposed matrix and vice versa and these interfaces all support transpose as an argument, so there is no need for reordering. And if you do np.dot(A.T, A) there is no reordering either, just a change of flags to DGEMM.
        m_mueller 10y ago
        There may not be an absolute need in terms of ability to perform operations, but performance is heavily dependant on accessing in the correct order - otherwise your accessible memory bandwidth goes down the drain, as do cache hitrates. I'm guessing LAPACK might not be too dependant on order though, but I haven't studied its performance tbh.
        dagss 10y ago
        Algorithms to handle what you talk about efficiently is one of the primary reason these libraries even exist. OF COURSE they access memory in the right order; that is their job. But it is not really about right vs wrong order, but using an algorithm to tile the data efficiently. Google "Anatomy of High Performance Matrix Multiplication" for examples from OpenBLAS. This is where NumPy falls through -- "a.T + a" will be very slow no matter what you do, but it can be done efficiently with tiling.
shoyer 10y ago
> Plus, Fortran already has way better multi-dim array handling built-in - NumPy is basically just a wrapper around that. I don't think either of these statements is true. Fortran is hard to beat for performance, but from a usability perspective NumPy (and other Python libraries) has some nice features for array-oriented computation that you can't find even in modern Fortran. The most obvious examples are broadcasting [1] and advanced indexing [2]. NumPy itself contains only C, though some SciPy routines wrap Fortran code. Yes, NumPy drew heavy inspiration from Fortran, but it's hardly a wrapper. [1] http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html http://docs.scipy.org/doc/numpy/user/basics.broadcasting.htm... [2] http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html http://docs.scipy.org/doc/numpy/reference/arrays.indexing.ht...
- m_mueller 10y ago
  I guess by "indexing" you mean what NumPy calls "advanced indexing". All the basic stuff seems lifted straight from Fortran. But yes, I should probably have been more exact with my statement - NumPy sure has some more helpful syntax that saves some LOC.
  - shoyer 10y ago
    Yes, I did mean "advanced indexing" -- I edited my post to clarify.

greglindahl 10y ago

I'm a huge Fortran fan, but even I can see that reducing the number of languages required to get high performance can be a good thing. Not all of my science colleagues are interested in learning a ton of languages.

m_mueller 10y ago
Well, most of the scientific computing in Python is happening on NumPy. Which uses Fortran compatible data structures. I doubt that a C++ backed python extension will be compatible with that. For that reason, if anything, I'd advocate to transpile to Fortran rather than C++ if you want more comfort for scientific users.
- shoyer 10y ago
  NumPy natively supports both row-major (C/C++) and column-major (Fortran) arrays, so from the perspective of NumPy inter-operability there would be no advantage either way.
  - m_mueller 10y ago
    I didn't know that. How does NumPy use its Fortran kernels with row-major arrays? Reordering would have quite significant performance implications.
    - shoyer 10y ago
      NumPy has an abstract iterator interface so you don't need to write a kernel with a specific memory layout in mind: http://docs.scipy.org/doc/numpy-1.11.0/reference/c-api.iterator.html#simple-iteration-example http://docs.scipy.org/doc/numpy-1.11.0/reference/c-api.itera... For operations like adding two arrays or summing an array along an axis, this means that NumPy automatically iterates in the fastest order.
      - m_mueller 10y ago
        Interesting. How does it decide the storage order when initializing an array? I'd assume due to python's dynamic nature it can't optimize for later operations there, right?
        infinite8s 10y ago
        It defaults to row major but you can specify column major on creation.
    - dagss 10y ago
      NumPy certainly does not have anything like "Fortran kernels". Elementwise arithmetic is coded in C (and not very fast...). Linear algebra uses the BLAS and LAPACK interfaces, of these only LAPACK is usually written in Fortran these days... but anyway, a Fortran matrix is a C transposed matrix and vice versa and these interfaces all support transpose as an argument, so there is no need for reordering. And if you do np.dot(A.T, A) there is no reordering either, just a change of flags to DGEMM.
      - m_mueller 10y ago
        There may not be an absolute need in terms of ability to perform operations, but performance is heavily dependant on accessing in the correct order - otherwise your accessible memory bandwidth goes down the drain, as do cache hitrates. I'm guessing LAPACK might not be too dependant on order though, but I haven't studied its performance tbh.
        dagss 10y ago
        Algorithms to handle what you talk about efficiently is one of the primary reason these libraries even exist. OF COURSE they access memory in the right order; that is their job. But it is not really about right vs wrong order, but using an algorithm to tile the data efficiently. Google "Anatomy of High Performance Matrix Multiplication" for examples from OpenBLAS. This is where NumPy falls through -- "a.T + a" will be very slow no matter what you do, but it can be done efficiently with tiling.

shoyer 10y ago

> Plus, Fortran already has way better multi-dim array handling built-in - NumPy is basically just a wrapper around that. I don't think either of these statements is true. Fortran is hard to beat for performance, but from a usability perspective NumPy (and other Python libraries) has some nice features for array-oriented computation that you can't find even in modern Fortran. The most obvious examples are broadcasting [1] and advanced indexing [2]. NumPy itself contains only C, though some SciPy routines wrap Fortran code. Yes, NumPy drew heavy inspiration from Fortran, but it's hardly a wrapper. [1] http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html http://docs.scipy.org/doc/numpy/user/basics.broadcasting.htm... [2] http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html http://docs.scipy.org/doc/numpy/reference/arrays.indexing.ht...

m_mueller 10y ago
I guess by "indexing" you mean what NumPy calls "advanced indexing". All the basic stuff seems lifted straight from Fortran. But yes, I should probably have been more exact with my statement - NumPy sure has some more helpful syntax that saves some LOC.
- shoyer 10y ago
  Yes, I did mean "advanced indexing" -- I edited my post to clarify.