4 ms·

True for small matrices but above certain size you wouldn’t want to use these methods or loops anyway, and they’re not optimised if your matrix is sparse either

True for small matrices but above certain size you wouldn’t want to use these methods or loops anyway, and they’re not optimised if your matrix is sparse either. So practically you end up using a library and you’re back to basically being equivalent with C.

leereeves 3y ago
The Fortran matmul intrinsic doesn't use an optimal algorithm?
- duped 3y ago
  Even theoretically this is only possible if the dimensions and sparseness are known at compile time, which can be pretty heavy assumptions.
  - jcranmer 3y ago
    The implementation of every BLAS implementation I'm aware of basically looks like this: if (dims in range A) gemm_implA else if (dims in range B) gemm_implB else if (dims in range C) gemm_implC // several more implementations... The immediate call function is invariably a switch over several implementations, and by the time you reach that call, you have actual knowledge of dimensions.
    - duped 3y ago
      The branch there can be expensive for small matrices.
- kergonath 3y ago
  The compiler is responsible for the implementation of the matmul intrinsic. Most of them have a “decent for small matrices and easy to inline” version by default and an option to convert it to call to BLAS under certain conditions (like larger matrices). Then, the performance is that of your BLAS library. The assumption (quite good in my experience) is that people use the intrinsic naturally because it is easy and natural, but turn to BLAS calls when performance becomes important. Which is not the case for the vast majority of matrix multiplication, even in high performance computing. There is no universal optimal algorithm.
- deleted 3y ago
  [deleted]
andy99 3y ago
It's possible to link in an external BLAS library that will be used by the intrinsics instead. See https://fortran-lang.discourse.group/t/compiler-option-matmul-speedup/6518/8 https://fortran-lang.discourse.group/t/compiler-option-matmu... But yes I agree, at some point it's always going to be back to something more customized.
- microtherion 3y ago
  And this is where the performance starts converging again, because BLAS is also available from C. So if you can express your computation in terms of BLAS primitives, it does not matter much from what language you're calling them.
  - kergonath 3y ago
    Yes, but you retain Fortran’s expressiveness. You can write a matrix product A = x * B, regardless of whether the matrix multiplication is handled by a compiler’s multiplication function, any flavour of BLAS, or ScaLAPACK over MPI. I don’t think Fortran has that much of an edge in performance compared to decently optimised C, but the Fortran code to get there is much more readable and easier to grasp. And with much fewer footguns (think allocatables versus pointer arithmetic).
    - lmm 3y ago
      If your goal is an expressive, safe, low-footgun language for expressing BLAS computations, you'd probably use NumPy rather than Fortran.
      - kergonath 3y ago
        NumPy is very Fortran-like, it's probably not the best example. Sure, it avoids some Fortran issues, but it also inherits some baggage from Python. Personally, working with both regularly, I don't really find it more expressive or safer than the common Fortran equivalents. And in Fortran there's less mental overhead about whether a variable is a list or a NumPy array, which does not matter until it does. Ditto for row/column-major order. Overall, in terms of ergonomics and ease of use, both are fairly similar.

leereeves 3y ago

The Fortran matmul intrinsic doesn't use an optimal algorithm?

duped 3y ago
Even theoretically this is only possible if the dimensions and sparseness are known at compile time, which can be pretty heavy assumptions.
- jcranmer 3y ago
  The implementation of every BLAS implementation I'm aware of basically looks like this: if (dims in range A) gemm_implA else if (dims in range B) gemm_implB else if (dims in range C) gemm_implC // several more implementations... The immediate call function is invariably a switch over several implementations, and by the time you reach that call, you have actual knowledge of dimensions.
  - duped 3y ago
    The branch there can be expensive for small matrices.
kergonath 3y ago
The compiler is responsible for the implementation of the matmul intrinsic. Most of them have a “decent for small matrices and easy to inline” version by default and an option to convert it to call to BLAS under certain conditions (like larger matrices). Then, the performance is that of your BLAS library. The assumption (quite good in my experience) is that people use the intrinsic naturally because it is easy and natural, but turn to BLAS calls when performance becomes important. Which is not the case for the vast majority of matrix multiplication, even in high performance computing. There is no universal optimal algorithm.
deleted 3y ago
[deleted]

andy99 3y ago

It's possible to link in an external BLAS library that will be used by the intrinsics instead. See https://fortran-lang.discourse.group/t/compiler-option-matmul-speedup/6518/8 https://fortran-lang.discourse.group/t/compiler-option-matmu... But yes I agree, at some point it's always going to be back to something more customized.

microtherion 3y ago
And this is where the performance starts converging again, because BLAS is also available from C. So if you can express your computation in terms of BLAS primitives, it does not matter much from what language you're calling them.
- kergonath 3y ago
  Yes, but you retain Fortran’s expressiveness. You can write a matrix product A = x * B, regardless of whether the matrix multiplication is handled by a compiler’s multiplication function, any flavour of BLAS, or ScaLAPACK over MPI. I don’t think Fortran has that much of an edge in performance compared to decently optimised C, but the Fortran code to get there is much more readable and easier to grasp. And with much fewer footguns (think allocatables versus pointer arithmetic).
  - lmm 3y ago
    If your goal is an expressive, safe, low-footgun language for expressing BLAS computations, you'd probably use NumPy rather than Fortran.
    - kergonath 3y ago
      NumPy is very Fortran-like, it's probably not the best example. Sure, it avoids some Fortran issues, but it also inherits some baggage from Python. Personally, working with both regularly, I don't really find it more expressive or safer than the common Fortran equivalents. And in Fortran there's less mental overhead about whether a variable is a list or a NumPy array, which does not matter until it does. Ditto for row/column-major order. Overall, in terms of ergonomics and ease of use, both are fairly similar.