This study presents a computational framework for the simultaneous optimization of the positions and lengths of several interacting geothermal borehole heat exchangers (BHEs). Thermal consistency is ensured by coupling a Finite Line Source (FLS) model for soil temperature evolution with a local ODE-based model governing fluid inlet and outlet temperatures. We minimize the uniform BHE length under operational, environmental, and geometric inequality constraints, effectively handling partially non-convex site geometries. Our implementation uses first and second derivatives for all constraints. To further enhance performance, two key acceleration strategies are presented. First, an efficient interpolation scheme is introduced for the radial dependence of the FLS solution, enabling rapid spatial evaluations over a wide range of borehole spacings. Second, Fast Fourier Transforms (FFT) are employed to accelerate the long-term temporal convolutions inherent to transient thermal simulations. The proposed framework is demonstrated through a case study, highlighting its flexibility and high computational efficiency.