Solving very large, very sparse linear systems over GF(2)?

Hello! I am a graduate student at UC Davis, and my area of research is
modeling graph-theoretic NP-Complete problems as systems of polynomial
equations. These problems reduce to solving a large-scale VERY sparse
linear system over GF(2).

I am very new to this area, and am trying to understand what is the
best algorithm for my matrices.

We are trying to solve Ax = b over GF(2) (which may or may not be
feasible) where A is m x n. Sometimes m > n, and sometimes m < n, but
they are almost never square. I don't have very much information about
the rank, but I know that columns either have two 1s and the rest
zero, or three 1s and the rest zero. The matrices get very, very
large... easily up into the millions, but because they are so sparse,
we can easily solve them using straightforward gaussian elimination.

However, I am looking for a more efficient method then straight-
forward gaussian elimination. I have been looking into the Montgomery
block Lanczos algorithm, the Coppersmith block Lanczos algorithm, and
the Wiedemann block algorithm. I have the following papers "Solving
Large Sparse Linear Systems over Finite Fields" (LaMacchia, Odlyzko),
"Solving Linear Equations over GF(2): Block Lanczos
Algorithm" (Coppersmith), "Solving Sparse Linear Equations over Finite
Fields" (Wiedemann), and "A Block Lanczos Algorithm for Finding
Dependencies over GF(2)" (Montgomery).

I don't know if this is a silly question, but I wonder what is the
best algorithm for me? Obviously, I would only like to implement once
and get it right the first time... I wonder if someone has a strong
opinion that they would mind quickly sharing? Stuctured elimination
with Lanczos? Block Wiedemann? Conjugate gradient? So many new
algorithms to learn and explore!

I would love a little bit of advice and direction as I jump into this
new area, if anyone has the time.

Very best regards,
Susan

.