Utilizing a modified embedded atom method potential, we performed large-scale classical molecular dynamics simulations (with up to 50 million atoms) to investigate the response of Ga single crystals to shock compression along the three major orientations of the orthorhombic A11 ground state, i.e., [001], [010], and [100]. For weak shocks with particle velocity u p < 300 m/s, these defect-free single crystals respond elastically, but for stronger shocks, they undergo a structural phase transformation and then (for even stronger shocks) melt. For intermediate shock strengths (300 m/s < u p < 1.2 km/s) a split shock wave is formed, with an elastic precursor (uniaxial compression wave) leading the slower transformation wave. The transformed region consists of a mixed phase, with stripes of the product phase embedded into the uniaxially compressed parent phase, with the ratio of product to parent phase increasing with increasing shock strength, much like in martensitic phase transformations. Upon shock release from the free surface at the end of the sample, the transformation is reversed, leaving only some defects behind, which makes it difficult to experimentally investigate the structural transformation in shock-recovered samples. We investigated the structure produced by the shock and found it to be similar to the β phase of Ga which is obtained by supercooling from the liquid state. However, the product phase shows only half the period along [100] in the calculated diffraction pattern. Further investigation showed that this is due to a different stacking sequence in this direction, namely ABCD instead of AB for the β phase.