  We present an efficient algorithm for calculating the properties of Ising
models in two dimensions, directly in the spin basis, without the need for
mapping to fermion or dimer models. The algorithm gives numerically exact
results for the partition function and correlation functions at a single
temperature on any planar network of N Ising spins in O(N^{3/2}) time or less.
The method can handle continuous or discrete bond disorder and is especially
efficient in the case of bond or site dilution, where it executes in O(L^2 ln
L) time near the percolation threshold. We demonstrate its feasibility on the
ferromagnetic Ising model and the +/- J random-bond Ising model (RBIM) and
discuss the regime of applicability in cases of full frustration such as the
Ising antiferromagnet on a triangular lattice.
