A subaperture stitching algorithm was developed for testing aspheric surfaces. The full aperture was divided into one central circular region plus several partially-overlapping annuli. Each annulus was composed of partially-overlapping circular subapertures. The phase map in each subaperture was obtained through the phase-shifting interferometry and retrieved by an iterative tilt-immune phase-shifting algorithm and a Zernike-polynomial-based phase-unwrapping process. All subapertures in one annulus were stitched simultaneously in least-squares sense. By eliminating the relative piston and tilt between adjacent subapertures, the sum of squared errors in the overlapped regions was minimized. The phase stitching between annuli also utilized the least-squares method in the overlapped region. Simulation results on a test wavefront with 30-wave spherical aberrations demonstrated the effectiveness of the proposed algorithm. The rms phase residue after the phase-shifting, phase-unwrapping and phase-stitching processes was 0.006 waves, which met the precision requirement of common interferometers. This algorithm should be applicable to general surfaces in subaperture stitching interferometry.