Linear Combination C++ in modulus

I want to compute an LU decomposition of a matrix and extract the linear combinations out of it.

I first asked a question using the library Armadillo here but as pointed out by one comment, Armadillo is not able to deal with modulus computation.

Thus I started to develop an LU using prime modulus from scratch, here is what I obtain but there is still a bug that I am not able to see.

Here is the code that I have for now. (Do not think too much about the class Matrix, it is just a way to encapsulate for now a vector<vector<int>>.

Matrix* Matrix::triangulation(Matrix & ident)
{
    unsigned int n = getNbLines();
    unsigned int m = getNbColumns();

    vector<vector<int>> mat = getMat();

    vector<vector<int>> identity = ident.getMat();
    vector<vector<int>> lower;
    vector<vector<int>> upper;

    /* ------------------------------------------------------------ */  

    /**
     * @brief
     * This code initialize a 'lower' matrix of size 'n x n'.
     * The matrix is fill with only '0'.
     */
    for(unsigned int i = 0; i < n; i++) {   
        vector<int> v(m);
        lower.push_back(v);

        for(unsigned int j = 0; j < m; j++) lower[i][j] = 0;            
    }   

    /**
     * @brief
     * This code initialize an 'upper' matrix of size 'n x m'.
     * The matrix is fill with only '0'.
     */
    for(unsigned int i = 0; i < n; i++) {   
        vector<int> v(m);
        upper.push_back(v);
        for(unsigned int j = 0; j < m; j++) upper[i][j] = 0;
    }

    /**
     * @brief
     * This code initialize an 'identity' matrix of size 'm x m'.
     * The matrix is fill with only '0'.
     */
    for(unsigned int i = 0; i < m; i++) {
        vector<int> v2(m);
        identity.push_back(v2);
        for(unsigned int j = 0; j < m; j++) identity[i][j] = 0;

        identity[i][i] = 1;
    }

    /* ------------------------------------------------------------ */

    // Decomposing matrix into Upper and Lower triangular matrix 
    for (unsigned int i = 0; i < n; i++) { 

        // Upper Triangular 
        for (unsigned int k = 0; k < m; k++) { 

            // Summation of L(i, j) * U(j, k) 
            int sum = 0; 
            for (unsigned int j = 0; j < n; j++) 
                sum = sum + ((lower[i][j] * upper[j][k])); 

            // Evaluating U(i, k)
            upper[i][k]    = (mat[i][k] - sum) % prime;
            identity[i][k] = (mat[i][k] - sum) % prime;
        } 

        // Lower Triangular 
        for (unsigned int k = 0; k < n; k++) { 
            if (i == k) {

                lower[i][i] = 1; // Diagonal as 1 
            }
            else { 

                // Summation of L(k, j) * U(j, i) 
                int sum = 0; 
                for (unsigned int j = 0; j < n; j++) 
                    sum = sum + ((lower[k][j] * upper[j][i])); 

                // Evaluating L(k, i) 
                lower[k][i] =    (((mat[k][i] - sum)) / upper[i][i]) % prime; 
                identity[k][i] = (((mat[k][i] - sum)) / upper[i][i]) % prime;                   
            } 
        } 
    }  

    ident.setMat(identity);

    return new Matrix(lower,prime);
}

I call it with the object : Matrix mat({ { 2, 1, 3, 2, 0}, { 4, 3, 0, 1, 1 }},5); So basically, I want the LU decomposition (especially the lower-triangle matrix) with all my computation done in modulus 5.

It works to extract the lower-matrix, however, the linear combinations (which are just all the operations done on an identity matrix) are not correct. Here is the trace that I have with an explanation of what I want to obtain:

c |-------------------------------------------------------------------------------------------------------|
c | Prime Number: 5

c |-------------------------------------------------------------------------------------------------------|
c | Input Matrix: 

2 1 3 2 0 
4 3 0 1 1 

c |-------------------------------------------------------------------------------------------------------|
c | Lower Matrix: 

1 0 0 0 0 
2 1 0 0 0 

c |-------------------------------------------------------------------------------------------------------|
c | Linear Combination Matrix:

2 0 3 2 0 
0 1 0 3 1 
0 0 1 0 0 
0 0 0 1 0 
0 0 0 0 1 

c |-------------------------------------------------------------------------------------------------------|
c | Expected Solution: 

3 2 3 0 3
0 1 1 3 4
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1

c |-------------------------------------------------------------------------------------------------------|
c | Explanations: 

c | 3 * c1 + 0 * c2 + 0 * c3 + 0 * c4 + 0 * c5 = c1 of Lower-Matrix
c | 2 * c1 + 1 * c2 + 0 * c3 + 0 * c4 + 0 * c5 = c2 of Lower-Matrix
c | 3 * c1 + 1 * c2 + 1 * c3 + 0 * c4 + 0 * c5 = c3 of Lower-Matrix
c | 0 * c1 + 3 * c2 + 0 * c3 + 1 * c4 + 0 * c5 = c4 of Lower-Matrix
c | 3 * c1 + 4 * c2 + 0 * c3 + 0 * c4 + 1 * c5 = c5 of Lower-Matrix

c +=======================================================================================================+

So as a small sum-up:

The lower-matrix is OK, the result is the expected one.
The linear combinations outputted are not the one expected.
An explanation of what is expected is given in the last subsection.

Question: Where is the bug in my way to apply the modification on the identity matrix, and why I do not get in output the correct linear combinations?

EDIT

Clear view of what should happened normally. But the algorithm I did (the LU decomposition) is not exactly the one I do by hand, even though it should lead to the same results. That is the real trouble here...

Let's put my comments into an actual answer: while addition and multiplication modulo prime will do what you expect (note below), subtraction has the pitfall where modulo will return negative results for negative input (e.g. (-3)%5 == -3) and for division you cannot just use integer division, you have to actually implement the inverse of multiplication (for hints, see the answer by Demosthenes in the linked previous question).

Note: unless you overflow, if prime*prime > INT_MAX , you are in trouble for multiplication as well