I'm looking for a library that performs matrix operations on large sparse matrices w/o sacrificing numerical stability. Matrices will be 1000+ by 1000+ and values of the matrix will be between 0 and 1000. I will be performing the index calculus algorithm so I will be generating (sparse) row vectors of the matrix serially. As I develop each row, I will need to test for linear independence. Once I fill my matrix with the desired number of linearly independent vectors, I will then need to transform the matrix into reduced row echelon form.
The problem now is that my implementation uses Gaussian elimination to determine linear independence (ensuring row echelon form once all my row vectors have been found). However, given the density and size of the matrix, this means the entries in each new row become exponentially larger over time, as the lcm of the leading entries must be found in order to perform cancellation. Finding the reduced form of the matrix further exacerbates the problem.
So my question is, is there an algorithm, or better yet an implementation, that can test linear independence and solve the reduced row echelon form while keeping the entries as small as possible? An efficient test for linear independence is especially important since in the index calculus algorithm it is performed by far the most.