Least-squares (LS) problems are encountered in many geophysical estimation and data analysis problems where a large number of observations (data) are combined to determine a model (some aspect of the earth structure). Examples of least squares in seismic exploration include several data processing algorithms, theoretically accurate LS migration, inversion for reservoir parameters, and background velocity estimation. A frequently encountered problem is that the volume of data in 3D is so large that the matrices required for the LS solution cannot be stored within the memory of a single computer. A new technique is described for parallel computation of the LS operator that is based on a partitioned-matrix algorithm. The classical LS method for solution of block-Toeplitz systems of normal equation (NE) to the general case of block-Hermitian and non-Toeplitz systems of NE. is generalized. Specifically, a solution of a block-Hermitian system of NE is shown that may be obtained recursively by linearly combining the solutions of lesser order that are related to the forward and backward subsystems of equations. This results in an efficient parallel algorithm in which each partitioned system can be evaluated independently. The application of the algorithm to the problem of 3D plane wave transformation is demonstrated.