This paper proposes a design of matrix machines for solving large-scale band-limited matrix equations. The design provides a guide-line for specifying functionality of each processing element PE, array arrangements of PEs, communication structures and parallel algorithms. The basic feature of the design is (1) a rotary array arrangement of data-flow type PEs with the slant width 2B-1 and circumference Z, on which the diagonal strip of a band-limited matrix data is wound around, and (2) wavefront-sequential executions of tasks, each of which is consisted of a group of pipelined wavefronts and circulates around the rotary array. The wavefront algorithm first produces the Gaussian solution via LU decomposition and then its accuracy is improved by applying the incomplete LU (conjugate) residual iteration.