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