We describe elements of a fast integral equation solver for large periodic and partly periodic finite array systems. A key element of the algorithm is utilization (in a rigorous way) of a block-Toeplitz structure of the impedance matrix in conjunction with either conventional Method of Moments (MoM), Fast Multipole Method (FMM), or Fast Fourier Transform (FFT)-based Adaptive Integral Method (AIM) compression techniques. We refer to the resulting algorithms as the (block-)Toeplitz-MoM, (block-)Toeplitz-AIM, or (block-)Toeplitz-FMM algorithms. While the computational complexity of the Toeplitz-AIM and Toeplitz-FMM algorithms is comparable to that of their non-Toeplitz counterparts, they offer a very significant (about two orders of magnitude for problems of the order of five million unknowns) storage reduction. In particular, our comparisons demonstrate, that the Toeplitz-AIM algorithm offers significant advantages in problems of practical interest involving arrays with complex antenna elements. This result follows from the more favorable scaling of the Toeplitz-AIM algorithm for arrays characterized by large number of unknowns in a single array element and applicability of the AIM algorithm to problems requiring strongly sub-wavelength resolution.