There are two simple approaches.

**1. If you prefer to do all of your algebra modulo ***k*,
the best way to do this is to compute the Smith Normal form of the matrix *A* (or more generally, a related matrix). You can do this by multiplication of unimodular matrices on the left and right of *A* (*i.e.* elementary row and column operations respectively) to compute greatest common divisors.

- First, compute the greatest common divisor of the first column, together with
*k*, in the (1,1) coefficient (which will change the value of the first row in general).
- Then, compute the greatest common divisor of the first row, together with
*k*, again producing the result in the (1,1) coefficient (which will change the value of the first column in general).

... And so forth, until the (1,1) coefficient divides all coefficients in the first row and also the first column; then clear the first row and first column. Repeat this for the second row/column, etc. until you have only coefficients on the diagonal; and all the while, keep track of the products of unimodular matrices which you accumulate on the left and right. Doing so yields unimodular matrices *U,V* such that *S = VAU* has only elements on the diagonal; and the system *A***x** ≡ **b** (mod *k*) is then equivalent to the system *S***z** ≡ *V***b** (mod *k*), where we perform the substitution **z** = *U*^{−1}**x**. This system is solvable if and only if each row individually is solvable; and for any solution **z**, we may obtain a solution **x** = *U***z** of the original system.

**2. An equivalent approach which is not restricted to artihmetic mod ***k**,* which is likely to be faster in practise, is to reduce to Diophantine equations. To wit: a linear congruence of the form
$$ \mathbf a_j \cdot \mathbf x \equiv b_j \pmod{k}$$
should instead be interpreted as the equivalent *Diophantine* equation
$$ \mathbf a_j \cdot \mathbf x + k d_j = b_j \\;,$$
introducing a new auxiliary variable *d*_{j} for the multiple of *k* which is neglected in the modular arithmetic. One can express an entire system *A***x** ≡ **b** (mod *k*) of such equations — involving the same modulus; but different moduli for each equation can be easily accomodated — by the system of integer equations
$$ \Bigl[\\;\\; A \\;\\;\Big|\\;\\; k I \;\;\Bigr] (\mathbf x \oplus \mathbf d) = \mathbf b \;.$$
In reference to Felix's comment above, you then compute a Hermite Normal Form (or any upper-triangular matrix, really) using unimodular transformations. If *C* = [ *A* | *kI* ], then one looks for the unimodular transformation *U* such that *C' = CU* is in Hermite normal form (in fact, any upper-triangular form will do); this may be done incrementally by column reduction to clear out rows, starting with the last one. One may then solve the system of equations
$$ C' \mathbf z = \mathbf b $$
over the reals, performing the substitution **z** = *U*^{ -1 }(**x** ⊕ **d**). The final several coefficients of **z** will be fixed, and will be integer-valued if there are any solutions at all to the original system of equations; and the first several will be free parameters, albeit ones which should range only over ℤ. Then, we may solve for **x** by computing (**x** ⊕ **d**) = *U***z**.

[This answer is completely revised, replacing an answer which I deleted after realizing that it was wrong; I subsequently added a reference to the Smith Normal Form above.]