General vs segmented contraction
Unless one is using a fully decontracted a.k.a "primitive" basis set, one uses contracted Gaussian-type orbitals (cGTOs), which are given as a linear combination of the primitive Gaussian-type orbitals (pGTOs):
$\chi_i^\text{cGTO}({\bf r}) = \sum_{\alpha} d_{\alpha i} \chi_\alpha^\text{pGTO}({\bf r}) $.
Alternatively, you can write the cGTOs in terms of the pGTOs as a matrix equation: ${\bf X}^\text{cGTO} = {\bf X}^\text{pGTO} {\bf D}$.
Let's assume our primitive exponents are given as a vector $(\alpha_1, \alpha_2, \dots, \alpha_N)^{\rm T}$ with $\alpha_1 > \alpha_2 > \dots > \alpha_N$. In a generally contracted basis set, ${\bf D}$ typically looks like
$\boldsymbol{D}=\left(\begin{array}{cccccc}
x & \cdots & x & x & 0 & 0\\
\vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\
x & x & x & x & 0 & 0\\
x & x & x & x & 0 & 0\\
x & x & x & x & 1 & 0\\
x & x & x & x & 0 & 1
\end{array}\right)$
where $x$ marks a non-zero entry. You have $K$ primitives (rows in D), which are contracted to $N$ atomic orbitals (columns filled with the $x$s), and then you have freed up the few most diffuse functions as primitive functions to be able to describe changes in the valence region.
Often the $N$ atomic orbitals are orthonormal to each other (${\bf D}^{\rm T} {\bf SD}={\bf 1}$ where ${\bf S}$ is the overlap), but this is not always the case; IIRC the correlation consistent sets for the transition metals have significant overlap between the contracted functions, since the first functions come from atomic Hartree-Fock calculations, while the last come from configuration interaction calculations and are not orthogonal to the Hartree-Fock orbitals.
The thing about generally contracted basis sets is that there is some flexibility in the way you can pick the functions. All that matters in the end is the span of the basis, so you are free to make transformations to your basis as long as the span stays the same. This means you're free to substract columns from each other: for example, if $x$ and $y$ form a linearly independent basis, $x$ and $x-y$ also span the same basis.
Now, note that I have two free primitives in the basis, which means these functions can be fully represented by the basis. I can just drop these functions from the general contraction, since I can eliminate the coefficients by substracting the free primitives with the correct weight from the contracted functions. (This is called a Davidson rotation and is what "optimize general contractions" mean on the Basis Set Exchange; see Chem. Phys. Lett. 1996, 260, 514−518.)
I now get
$\boldsymbol{D}=\left(\begin{array}{cccccc}
x & \cdots & x & x & 0 & 0\\
\vdots & \ddots & \vdots & \vdots & \vdots & \vdots\\
x & x & x & x & 0 & 0\\
x & x & x & x & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0\\
0 & 0 & 0 & 0 & 0 & 1
\end{array}\right)$
But, it does not stop here: I can eliminate entries in the fully contracted block, too, by adding and substracting vectors. Let's say the contracted block looks like this:
$\boldsymbol{D}=\left(\begin{array}{cccc}
x & x & x & x\\
x & x & x & x\\
x & x & x & x\\
x & x & x & x\\
x & x & x & x\\
x & x & x & x\\
x & x & x & x\\
x & x & x & x
\end{array}\right)$
I can take the first column from the left and substract it from the other ones to make the coefficient of the tightest exponent vanish. Next, I take the second column and make the coefficient of the second-tightest function go away.. and keep going until I've made it to the right-most column.
I can also do the same thing in reverse: I take the last column and make the most diffuse coefficient vanish in the other columns, and repeat until I am back on the right. Then I get
$\boldsymbol{D}=\left(\begin{array}{cccc}
x & 0 & 0 & 0\\
x & x & 0 & 0\\
x & x & x & 0\\
x & x & x & x\\
x & x & x & x\\
0 & x & x & x\\
0 & 0 & x & x\\
0 & 0 & 0 & x
\end{array}\right)$
Suddenly, you see that the contraction matrix has developed quite a few zeros. If the contraction length is not very long, $N\approx K$, the columns are mostly zeros (the extreme case is $N=K$ where you can choose ${\bf D}={\bf 1}$ i.e. to just use uncontracted primitives).
Moreover, one can get even more sparsity by discarding small elements $x$ in the columns, and (possibly) reoptimizing the contraction coefficients and exponents for a pre-specified sparsity pattern, e.g.
$\boldsymbol{D}=\left(\begin{array}{cccc}
x & 0 & 0 & 0\\
x & 0 & 0 & 0\\
x & x & 0 & 0\\
0 & x & 0 & 0\\
0 & x & x & 0\\
0 & 0 & x & x\\
0 & 0 & x & x\\
0 & 0 & 0 & x
\end{array}\right)$
Frank Jensen has written a very nice article on this a few years ago, see J. Chem. Theory Comput. 2014, 10, 3, 1074–1085.
Integral evaluation
What does this have to do with integral evaluation? The two-electron integrals $(pq|rs)$ can be evaluated analytically for Gaussian basis functions. However, the analytical integrals are known only for the primitives, which means that to get the integrals for the contracted functions one has to calculate
$(pq|rs) = \sum_{\alpha \beta \gamma \delta} d_{\alpha p} d_{\beta q} d_{\gamma r} d_{\delta s} (\alpha \beta|\gamma \delta)$, instead, where $\alpha \beta \gamma \delta$ stand for the primitive functions.
If you have a segmented basis set, where there are a lot of zeros in the contraction matrix, it is more efficient to evaluate the integrals one shell quartet at a time; that is, you recompute all the primitive integrals for each $(pq|rs)$ you're aiming for, because most primitives only contribute to a single contracted basis function. This means that you have simpler code, but also has the additional benefit that you can easily perform the contractions within the two-electron integral evaluation algorithm, which makes some steps cheaper in the computation, depending on the used integrals evaluation algorithm. Indeed, most codes assume this approach; e.g. this is how Gaussian, Psi4, Q-Chem, and Orca operate.
However, if you have a generally contracted basis set (especially when $K \gg N$ and $N \gg 1$), the programs that rely internally on segmented basis sets run into troble, as instead of computing the $K^4$ primitive shell quartets, you basically need to compute the $K^4$ quartets $N^4$ times. Screening for combinations of primitives that have small contributions helps a bit, but this is still extremely slow.
As far as I understand, the generally contracted programs compute all primitive shell quartets once, and then form the $N^4$ contractions with the ${\bf D}$'s per the equation above to get the target integrals $(pq|rs)$. (Disclaimer: I have never written a general contraction code.)
Programs that handle generally contracted basis sets well are far fewer, but they do exist. For instance, OpenMolcas handles general contractions. PySCF is also very efficient with generally contracted basis set. Both programs are open source and freely available. NWChem is also open source and freely available, and claims to support general contractions (see the opening paragraph of their documentation on basis sets), but it's so slow that I'm not sure whether this is really the case; see below.
On the commercial side, IIRC Molpro has excellent support for general contractions, and is faster than OpenMolcas.
Benchmarks
To see whether your program supports general contractions is quite simple: just try running e.g. a transition metal system with a correlation consistent basis set. For instance, in cc-pVDZ the first-row transition metals have a [20s16p8d2f|6s5p3d1f] contraction pattern, so $K=20$ and $N=6$ for s orbitals, $K=16$ and $N=5$ for p orbitals, and $K=8$ and $N=3$ for d orbitals, which make them quite heavy.
Although I'm usually against mindless benchmarks, here one makes sense. I'll take the Zn atom and Zn2 dimer (R=2.3 Å) as examples, since they're nice closed-shell systems.
I performed all benchmarks on my laptop with 4 cores, using Fedora 32 packages: psi4-1.3.2-2.fc32.x86_64, python3-pyscf-1.7.5-1.fc32.x86_64 employing the optimized integrals library qcint-4.0.5-1.fc32.x86_64, nwchem-openmpi-7.0.2-1.fc32.x86_64, and OpenMolcas-20.10-1.fc32.x86_64.
Psi4
As I said above, Psi4 is an example of a program that assumes segmented basis sets.
Setting scf_type = pk to enforce exact integrals instead of the default use of density fitting (which would make the calculation a LOT faster), the HF calculation on the zinc atom took about 19 seconds wall time. Final energy -1777.84665520791896
Running the dimer took considerably longer: 216 seconds of wall time. Final energy -3555.64477390092725
PySCF
The PySCF calculation on the Zn atom runs in about 0.7 seconds of wall time. Final energy -1777.846655207936010
The Zn2 dimer took about 1.4 seconds wall time. Final energy -3555.644773900906785
OpenMolcas
Running the Zn atom took about 2 seconds. Final energy -1777.846655208
The Zn2 dimer took about 5 seconds. Final energy -3555.644773901
NWChem
After running for 5 minutes on four cores for the Zn atom, without seeing even one iteration of SCF, I cancelled the execution...
