nmod_mat_solve_tril_right_classical, nmod_mat_solve_triu_right_classical

Substitution, for an $r \times r$ triangular matrix, with three methods. If $B$ has at least $\max(6, r/4)$ rows ($\max(6, r/2)$ for $n > 2^{31}$, or $n > 2^{15}$ with 32-bit words), as nmod_mat_solve_tril_classical (nmod_mat_solve_triu_classical) with delayed reductions on the transposed system, $B$ being transposed by blocks of 256 rows. Otherwise, if $r \le 128$ ($r \le 64$ for these larger moduli) or $B$ has at most 4 rows, each row $x$ of $X$ is computed from the last entry to the first for $L$ (from the first to the last for $U$): once $x_j$ is known, $x_j$ times the row $j$ of the triangular matrix is added to the unreduced sums of the remaining entries, each reduced once. Otherwise (a large triangular matrix and a few rows), each entry is computed with one dot product and one reduction, the columns of the triangular matrix being copied to rows 64 at a time so that the dot products read contiguous memory. The inverses of the diagonal entries are computed once.