stdlib_linalg_kronecker.fypp Source File

#:include "common.fypp"
#:set RCI_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES + INT_KINDS_TYPES
submodule (stdlib_linalg) stdlib_linalg_kronecker

  implicit none

contains

  #:for k1, t1 in RCI_KINDS_TYPES
    pure module function kronecker_product_${t1[0]}$${k1}$(A, B) result(C)
      ${t1}$, intent(in) :: A(:,:), B(:,:)
      ${t1}$ :: C(size(A,dim=1)*size(B,dim=1),size(A,dim=2)*size(B,dim=2))
      integer :: m1, n1, maxM1, maxN1, maxM2, maxN2
      
      maxM1 = size(A, dim=1)
      maxN1 = size(A, dim=2)
      maxM2 = size(B, dim=1)
      maxN2 = size(B, dim=2)
      

      do n1 = 1, maxN1
         do m1 = 1, maxM1
            ! 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
    end function kronecker_product_${t1[0]}$${k1}$
  #:endfor
end submodule stdlib_linalg_kronecker