Simulating large, regular rectangular arrays with equidistant interspacing is challenging as the computational complexity grows quickly with array size. More so when considering multiple excitations, for example, to find embedded element patterns (EEPs). However, the array, appropriately meshed, leads to a multilevel Toeplitz structure in the method-of-moment (MoM) impedance matrix that can be used to mitigate the increased complexity. The problem is further complicated when elements are electrically connected, requiring specialized methods to handle currents flowing between elements. These connections also lead to strong coupling, putting high demands on the solvers used. This article presents two different accelerated solvers that both utilize the matrix structure via a novel mesh-partitioning algorithm, reducing storage and computational costs. The first is an iterative method based on multilevel fast Fourier transform (MLFFT), with an option to solve for all excitations at once, reducing solution time at the cost of memory. The MLFFT acceleration is described in detail. The second is based on a fast direct Toeplitz solver, adapted to a block-matrix structure, and shows excellent stability in scaling. Both methods are evaluated on two different array element types, for arrays with up to 900 elements. The results are compared with conventional direct and iterative matrix solvers. Improvements are seen in both time and required storage. The choice of the most efficient method depends on the problem parameters. Two different preconditioners are evaluated for connected elements. The accelerated methods handle connected elements while greatly outperforming regular matrix inversion methods with negligible reduction in accuracy.