feat: extend intrinsic matmul

example/intrinsics/CMakeLists.txt

@@ -1,2 +1,3 @@
 ADD_EXAMPLE(sum)
-ADD_EXAMPLE(dot_product)
+ADD_EXAMPLE(dot_product)
+ADD_EXAMPLE(matmul)

example/intrinsics/example_matmul.f90

@@ -0,0 +1,7 @@
+program example_matmul
+ use stdlib_intrinsics, only: stdlib_matmul
+ complex :: a(2,2)
+ a = reshape([(0, 0), (0, -1), (0, 1), (0, 0)], [2, 2]) ! pauli y-matrix
+
+ print *, stdlib_matmul(a, a, a, a, a) ! should be sigma_y
+end program example_matmul

src/CMakeLists.txt

@@ -19,6 +19,7 @@ set(fppFiles
  stdlib_hash_64bit_spookyv2.fypp
  stdlib_intrinsics_dot_product.fypp
  stdlib_intrinsics_sum.fypp
+ stdlib_intrinsics_matmul.fypp
  stdlib_intrinsics.fypp
  stdlib_io.fypp
  stdlib_io_npy.fypp
@@ -32,14 +33,14 @@ set(fppFiles
  stdlib_linalg_kronecker.fypp
  stdlib_linalg_cross_product.fypp
  stdlib_linalg_eigenvalues.fypp
- stdlib_linalg_solve.fypp 
+ stdlib_linalg_solve.fypp
  stdlib_linalg_determinant.fypp
  stdlib_linalg_qr.fypp
  stdlib_linalg_inverse.fypp
  stdlib_linalg_pinv.fypp
  stdlib_linalg_norms.fypp
  stdlib_linalg_state.fypp
- stdlib_linalg_svd.fypp 
+ stdlib_linalg_svd.fypp
  stdlib_linalg_cholesky.fypp
  stdlib_linalg_schur.fypp
  stdlib_optval.fypp

src/stdlib_intrinsics.fypp

@@ -146,6 +146,39 @@ module stdlib_intrinsics
  #:endfor
  end interface
  public :: kahan_kernel
+
+ interface stdlib_matmul
+ !! version: experimental
+ !!
+ !!### Summary
+ !! compute the matrix multiplication of more than two matrices with a single function call.
+ !! ([Specification](../page/specs/stdlib_intrinsics.html#stdlib_matmul))
+ !!
+ !!### Description
+ !!
+ !! matrix multiply more than two matrices with a single function call
+ !! the multiplication with the optimal parenthesization for efficiency of computation is done automatically
+ !! Supported data types are `real`, `integer` and `complex`.
+ !!
+ !! Note: The matrices must be of compatible shapes to be multiplied
+ #:for k, t, s in I_KINDS_TYPES + R_KINDS_TYPES + C_KINDS_TYPES
+ pure module function stdlib_matmul_${s}$_3 (a, b, c) result(d)
+ ${t}$, intent(in) :: a(:,:), b(:,:), c(:,:)
+ ${t}$, allocatable :: d(:,:)
+ end function stdlib_matmul_${s}$_3
+
+ pure module function stdlib_matmul_${s}$_4 (a, b, c, d) result(e)
+ ${t}$, intent(in) :: a(:,:), b(:,:), c(:,:), d(:,:)
+ ${t}$, allocatable :: e(:,:)
+ end function stdlib_matmul_${s}$_4
+
+ pure module function stdlib_matmul_${s}$_5 (a, b, c, d, e) result(f)
+ ${t}$, intent(in) :: a(:,:), b(:,:), c(:,:), d(:,:), e(:,:)
+ ${t}$, allocatable :: f(:,:)
+ end function stdlib_matmul_${s}$_5
+ #:endfor
+ end interface stdlib_matmul
+ public :: stdlib_matmul
 
 contains
 

src/stdlib_intrinsics_matmul.fypp

@@ -0,0 +1,133 @@
+#:include "common.fypp"
+#:set I_KINDS_TYPES = list(zip(INT_KINDS, INT_TYPES, INT_KINDS))
+#:set R_KINDS_TYPES = list(zip(REAL_KINDS, REAL_TYPES, REAL_SUFFIX))
+#:set C_KINDS_TYPES = list(zip(CMPLX_KINDS, CMPLX_TYPES, CMPLX_SUFFIX))
+
+submodule (stdlib_intrinsics) stdlib_intrinsics_matmul
+ implicit none
+
+contains
+
+ ! Algorithm for the optimal parenthesization of matrices
+ ! Reference: Cormen, "Introduction to Algorithms", 4ed, ch-14, section-2
+ ! Internal use only!
+ pure function matmul_chain_order(n, p) result(s)
+ integer, intent(in) :: n, p(:)
+ integer :: s(1:n-1, 2:n), m(1:n, 1:n), l, i, j, k, q
+ m(:,:) = 0
+ s(:,:) = 0
+
+ do l = 2, n
+ do i = 1, n - l + 1
+ j = i + l - 1
+ m(i,j) = huge(1)
+
+ do k = i, j - 1
+ q = m(i,k) + m(k+1,j) + p(i)*p(k+1)*p(j+1)
+
+ if (q < m(i, j)) then
+ m(i,j) = q
+ s(i,j) = k
+ end if
+ end do
+ end do
+ end do
+ end function matmul_chain_order
+
+#:for k, t, s in I_KINDS_TYPES + R_KINDS_TYPES + C_KINDS_TYPES
+
+ pure module function stdlib_matmul_${s}$_3 (a, b, c) result(d)
+ ${t}$, intent(in) :: a(:,:), b(:,:), c(:,:)
+ ${t}$, allocatable :: d(:,:)
+ integer :: sa(2), sb(2), sc(2), cost1, cost2
+ sa = shape(a)
+ sb = shape(b)
+ sc = shape(c)
+
+ if ((sa(2) /= sb(1)) .or. (sb(2) /= sc(1))) then
+ error stop "stdlib_matmul: Incompatible array shapes"
+ end if
+
+ ! computes the cost (number of scalar multiplications required)
+ ! cost(A, B) = shape(A)(1) * shape(A)(2) * shape(B)(2)
+ cost1 = sa(1) * sa(2) * sb(2) + sa(1) * sb(2) * sc(2) ! ((AB)C)
+ cost2 = sb(1) * sb(2) * sc(2) + sa(1) * sa(2) * sc(2) ! (A(BC))
+
+ if (cost1 < cost2) then
+ d = matmul(matmul(a, b), c)
+ else
+ d = matmul(a, matmul(b, c))
+ end if
+ end function stdlib_matmul_${s}$_3
+
+ pure module function stdlib_matmul_${s}$_4 (a, b, c, d) result(e)
+ ${t}$, intent(in) :: a(:,:), b(:,:), c(:,:), d(:,:)
+ ${t}$, allocatable :: e(:,:)
+ integer :: p(5), i
+ integer :: s(3,2:4)
+
+ p(1) = size(a, 1)
+ p(2) = size(b, 1)
+ p(3) = size(c, 1)
+ p(4) = size(d, 1)
+ p(5) = size(d, 2)
+
+ s = matmul_chain_order(4, p)
+
+ select case (s(1,4))
+ case (1)
+ select case (s(2, 4))
+ case (2)
+ e = matmul(a, matmul(b, matmul(c, d)))
+ case (3)
+ e = matmul(a, matmul(matmul(b, c), d))
+ case default
+ error stop "stdlib_matmul: unexpected error unexpected s(i,j)"
+ end select
+ case (2)
+ e = matmul(matmul(a, b), matmul(c, d))
+ case (3)
+ select case (s(1, 3))
+ case (1)
+ e = matmul(matmul(a, matmul(b, c)), d)
+ case (2)
+ e = matmul(matmul(matmul(a, b), c), d)
+ case default
+ error stop "stdlib_matmul: unexpected error unexpected s(i,j)"
+ end select
+ case default
+ error stop "stdlib_matmul: unexpected error unexpected s(i,j)"
+ end select
+ end function stdlib_matmul_${s}$_4
+
+ pure module function stdlib_matmul_${s}$_5 (a, b, c, d, e) result(f)
+ ${t}$, intent(in) :: a(:,:), b(:,:), c(:,:), d(:,:), e(:,:)
+ ${t}$, allocatable :: f(:,:)
+ integer :: p(6), i
+ integer :: s(4,2:5)
+
+ p(1) = size(a, 1)
+ p(2) = size(b, 1)
+ p(3) = size(c, 1)
+ p(4) = size(d, 1)
+ p(5) = size(e, 1)
+ p(6) = size(e, 2)
+
+ s = matmul_chain_order(5, p)
+
+ select case (s(1,5))
+ case (1)
+ f = matmul(a, stdlib_matmul(b, c, d, e))
+ case (2)
+ f = matmul(matmul(a, b), stdlib_matmul(c, d, e))
+ case (3)
+ f = matmul(stdlib_matmul(a, b ,c), matmul(d, e))
+ case (4)
+ f = matmul(stdlib_matmul(a, b, c, d), e)
+ case default
+ error stop "stdlib_matmul: unexpected error unexpected s(i,j)"
+ end select
+ end function stdlib_matmul_${s}$_5
+
+#:endfor
+end submodule stdlib_intrinsics_matmul