Approximate matrix inverse of quadratic complexity ( A⁻¹ ≈ D⁻¹ - D⁻¹ E D⁻¹ )

ainv.m

function x = ainv(m)
% ainv - approximate inverse
% by Christian Himpe 2016
% released under BSD 2-Clause License ( opensource.org/licenses/BSD-2-Clause )

    d = diag(m);
    d(d~=0) = 1.0./d(d~=0);
    n = numel(d);
    x = bsxfun(@times,m,-d);
    x = bsxfun(@times,x,d');
    x(1:n+1:end) = d;
end