subroutine krprod (l1,l2,a,m1,m2,b,n1,n2,c)
      double precision a(1), b(1), c(1)
c     kronecker product c= a x b is formed.
c     a is l1 x l2, b is m1 x m2, c is n1 x n2, all matrices are
c     as vectors in column representation.
      n1=l1*m1
      n2=l2*m2
      n=1
      do 10 j=1,l2
      do 10 l=1,m2
      do 10 i=1,l1
      do 10 k=1,m1
      ii=l1*(j-1)+i
      jj=m1*(l-1)+k
      c(n)=a(ii)*b(jj)
      n=n+1
   10 continue
      return
      end