arraybase 1
dim a(2, 2)
a[1,1] = 1 : a[1,2] = 2 : a[2,1] = 3 : a[2,2] = 4
dim b(2, 2)
b[1,1] = 0 : b[1,2] = 5 : b[2,1] = 6 : b[2,2] = 7
call kronecker_product(a, b)

print
dim x(3, 3)
x[1,1] = 0 : x[1,2] = 1 : x[1,3] = 0
x[2,1] = 1 : x[2,2] = 1 : x[2,3] = 1
x[3,1] = 0 : x[3,2] = 1 : x[3,3] = 0
dim y(3, 4)
y[1,1] = 1 : y[1,2] = 1 : y[1,3] = 1 : y[1,4] = 1
y[2,1] = 1 : y[2,2] = 0 : y[2,3] = 0 : y[2,4] = 1
y[3,1] = 1 : y[3,2] = 1 : y[3,3] = 1 : y[3,4] = 1
call kronecker_product(x, y)
end

subroutine kronecker_product(a, b)
	ua1 = a[?][]
	ua2 = a[][?]

	ub1 = b[?][]
	ub2 = b[][?]

	for i = 1 to ua1
		for k = 1 to ub1
			print "[";
			for j = 1 to ua2
				for l = 1 to ub2
					print rjust(a[i, j] * b[k, l], 2);
					if j = ua1 and l = ub2 then
						print "]"
					else
						print " ";
					endif
				next
			next
		next
	next
end subroutine