I have to solve a system of up to 10000 equations with 10000 unknowns as fast as possible (preferably within a few seconds). I know that Gaussian elimination is too slow for that, so what algorithm is suitable for this task?All coefficients and constants are non-negative integers modulo p (where p is a prime). There is guaranteed to be only 1 solution.  I need the solution modulo p.
Fastest way to solve a system of linear equations
algorithms;linear algebra;matrices;numerical algorithms
A LU decomposition of a $n \times n$ matrix can be computed in $O(M(n))$ time, where $M(n)$ is the time to multiply two $n \times n$ matrices.  Therefore, you can find a solution to a system of $n$ linear equations in $n$ unknowns in $O(M(n))$ time.  For instance, Strassen's algorithm achieves $M(n) = O(n^{2.8})$, which is faster than Gaussian elimination.  See https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion.Rather than trying to implement this yourself, I would suggest using a library: e.g., a BLAS library.