Added spec and example program.

Author: adenchfi

doc/specs/stdlib_linalg.md

@@ -160,6 +160,37 @@ Returns a rank-2 array equal to `u v^T` (where `u, v` are considered column vect
 {!example/linalg/example_outer_product.f90!}
 ```
 
+## `kronecker_product` - Computes the kronecker product of two rank-2 arrays
+
+### Status
+
+Experimental
+
+### Description
+
+Computes the Kronecker product of two rank-2 arrays
+
+### Syntax
+
+`C = [[stdlib_linalg(module):kronecker_product(interface)]](A, B)`
+
+### Arguments
+
+`A`: Shall be a rank-2 array with dimensions M1, N1
+
+`v`: Shall be a rank-2 array with dimensions M2, N2
+
+### Return value
+
+Returns a rank-2 array equal to `A \otimes B`. The shape of the returned array is `[M1*M2, N1*N2]`.
+
+### Example
+
+```fortran
+{!example/linalg/example_kronecker_product.f90!}
+```
+
+
 ## `cross_product` - Computes the cross product of two vectors
 
 ### Status

example/linalg/example_kronecker_product.f90

@@ -0,0 +1,25 @@
+program example_kronecker_product
+  use stdlib_linalg, only: kronecker_product
+  implicit none
+  integer, parameter :: m1=1, n1=2, m2=2, n2=3
+  real :: A(m1, n1), B(m2,n2)
+  real, allocatable :: C(:,:)
+
+  do j=1, n1
+     do i=1, m1
+        A(i,j) = i*j ! A = [1, 2]
+     end do
+  end do
+  
+  do j=1, n2
+     do i=1, m2      ! B = [ 1, 2, 3 ]
+        B(i,j) = i*j !     [ 2, 4, 6 ]
+     end do
+  end do
+  
+  C = kronecker_product(A, B)
+  ! C =     [ a(1,1) * B(:,:) | a(1,2) * B(:,:) ]
+  ! or in other words, 
+  ! C =     [  1.00      2.00      3.00      2.00      4.00      6.00  ]
+  !         [  2.00      4.00      6.00      4.00      8.00     12.00  ]
+end program example_kronecker_product

src/stdlib_linalg_kronecker.fypp

@@ -19,7 +19,8 @@ contains
 
       do n1=1, maxN1
          do m1=1, maxM1
-            ! We use the numpy.kron convention for ordering of the matrix elements
+            ! We use the Wikipedia convention for ordering of the matrix elements
+	    ! https://en.wikipedia.org/wiki/Kronecker_product
             C((m1-1)*maxM2+1:m1*maxM2, (n1-1)*maxN2+1:n1*maxN2) = A(m1, n1) * B(:,:)
          end do
       end do