function M=SPAI(A,Z)
% SPAI computes a sparse approximate inverse
%   M=SPAI(A,Z); computes a sparse approximate inverse M of A with the 
%   sparsity pattern in the matrix Z using the SPAI technique.

n=size(A,2);
I=speye(n);
M=Z;
for i=1:n
  id=find(Z(:,i)~=0);                   % find non-zero entries in Z
  M(id,i)=A(:,id)\I(:,i);               % fill them with least squares 
end