Solving many small systems

function pixsolve(A, b) result(x)
    implicit none
    real*8    :: A(:,:,:), b(:,:), x(size(b,1),size(b,2))
    integer*4 :: i, n, m, piv(size(b,2)), err
    n = size(A,1)
    m = size(A,2)
    x = b
    do i = 1, n
        call dgesv(m, 1, A(i,:,:), m, piv, x(i,:), m, err)
    end do
end function


! f2py -c -m foo{,.f90} -llapack -lblas
! python
! >>> import numpy as np
! >>> from foo import pixsolve
! >>> n = 6
! >>> A = np.random.rand(n, 3, 3)
! >>> b = np.random.rand(n, 3)
! >>> x = pixsolve(A, b)  # slightly faster (5%) if A, b are in Fortran order
! >>> max(np.linalg.norm(np.dot(A[i], x[i]) - b[i]) for i in xrange(n))
! 2.603703785810335e-16