This study introduces a fragmentation-based linear-scaling method for strongly correlated systems, specifically the divide-and-conquer Hartree-–Fock-–Bogoliubov (DC-HFB) approach. Two energy gradient formulations of the DC-HFB method are derived and implemented, enabling efficient optimization of molecular geometries in large systems. This method is applied to graphene nanoribbons (GNRs) to explore their geometries and polyradical characters. Numerical results demonstrate that the present DC-HFB method has a potential to treat the static electron correlation and predict diradical character in GNRs, offering new avenues for studying large-scale strongly correlated systems.