feat: extend intrinsic matmul

doc/specs/stdlib_intrinsics.md

@@ -155,4 +155,42 @@ The output is a scalar of the same type and kind as to that of `x` and `y`.
 
 ```fortran
 {!example/intrinsics/example_dot_product.f90!}
-```
+```
+
+### `stdlib_matmul` function
+
+#### Description
+
+The extension of the intrinsic function `matmul` to handle more than 2 and less than or equal to 5 matrices, with error handling using `linalg_state_type`.
+The optimal parenthesization to minimize the number of scalar multiplications is done using the Algorithm as outlined in Cormen, "Introduction to Algorithms", 4ed, ch-14, section-2.
+The actual matrix multiplication is performed using the `gemm` interfaces.
+It supports only `real` and `complex` matrices.
+
+#### Syntax
+
+`res = ` [[stdlib_intrinsics(module):stdlib_matmul(interface)]] ` (m1, m2, m3, m4, m5, err)`
+
+#### Status
+
+Experimental
+
+#### Class
+
+Function.
+
+#### Argument(s)
+
+`m1`, `m2`: 2D arrays of the same kind and type. `intent(in)` arguments.
+`m3`,`m4`,`m5`: 2D arrays of the same kind and type as the other matrices. `intent(in), optional` arguments.
+`err`: `type(linalg_state_type), intent(out), optional` argument. Can be used for elegant error handling. It is assigned `LINALG_VALUE_ERROR` 
+        in case the matrices are not of compatible sizes.
+
+#### Result
+
+The output is a matrix of the appropriate size.
+
+#### Example
+
+```fortran
+{!example/intrinsics/example_matmul.f90!}
+```

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,12 @@
+program example_matmul
+    use stdlib_intrinsics, only: stdlib_matmul
+    complex :: x(2, 2), y(2, 2), z(2, 2)
+    x = reshape([(0, 0), (1, 0), (1, 0), (0, 0)], [2, 2])
+    y = reshape([(0, 0), (0, 1), (0, -1), (0, 0)], [2, 2]) ! pauli y-matrix
+    z = reshape([(1, 0), (0, 0), (0, 0), (-1, 0)], [2, 2])
+
+    print *, stdlib_matmul(x, y) ! should be iota*z
+    print *, stdlib_matmul(y, z, x) ! should be iota*identity
+    print *, stdlib_matmul(x, x, z, y) ! should be -iota*x
+    print *, stdlib_matmul(x, x, z, y, y) ! should be z
+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

@@ -8,6 +8,7 @@ module stdlib_intrinsics
     !!Alternative implementations of some Fortran intrinsic functions offering either faster and/or more accurate evaluation.
     !! ([Specification](../page/specs/stdlib_intrinsics.html))
     use stdlib_kinds
+    use stdlib_linalg_state, only: linalg_state_type
     implicit none
     private
 
@@ -146,6 +147,49 @@ 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` and `complex`.
+        !!
+        !! Note: The matrices must be of compatible shapes to be multiplied
+        #:for k, t, s in R_KINDS_TYPES + C_KINDS_TYPES
+            pure module function stdlib_matmul_pure_${s}$ (m1, m2, m3, m4, m5) result(r)
+                ${t}$, intent(in) :: m1(:,:), m2(:,:)
+                ${t}$, intent(in), optional :: m3(:,:), m4(:,:), m5(:,:)
+                ${t}$, allocatable :: r(:,:)
+            end function stdlib_matmul_pure_${s}$
+
+            module function stdlib_matmul_${s}$ (m1, m2, m3, m4, m5, err) result(r)
+                ${t}$, intent(in) :: m1(:,:), m2(:,:)
+                ${t}$, intent(in), optional :: m3(:,:), m4(:,:), m5(:,:)
+                type(linalg_state_type), intent(out) :: err
+                ${t}$, allocatable :: r(:,:)
+            end function stdlib_matmul_${s}$
+        #:endfor
+    end interface stdlib_matmul
+    public :: stdlib_matmul
+
+    ! internal interface
+    interface stdlib_matmul_sub
+        #:for k, t, s in R_KINDS_TYPES + C_KINDS_TYPES
+            pure module subroutine stdlib_matmul_sub_${s}$ (res, m1, m2, m3, m4, m5, err)
+                ${t}$, intent(out), allocatable :: res(:,:)
+                ${t}$, intent(in) :: m1(:,:), m2(:,:)
+                ${t}$, intent(in), optional :: m3(:,:), m4(:,:), m5(:,:)
+                type(linalg_state_type), intent(out), optional :: err
+            end subroutine stdlib_matmul_sub_${s}$
+        #:endfor
+    end interface stdlib_matmul_sub
 
 contains