Inverse of a binary matrix in C

Question

I have a binary matrix (zeros and ones) D[][] of dimension nxn where n is large (approximately around 1500 - 2000). I want to find the inverse of this matrix in C.

Since I'm new to C, I started with a 3 x 3 matrix and working around to generalize it to N x N. This works for int values, however since I'm working with binary 1's and 0's. In this implementation, I need unsigned int values.

I could find many solutions for int values but I didn't come across any solution for unsigned int. I'd like to find the inverse of a N x N binary matrix without using any external libraries like blas/lapack. It'd be great if anyone could provide a lead on M x N matrix.

Please note that I need inverse of a matrix, not the pseudo-inverse.

/* To find the inverse of a matrix using LU decomposition */

/* standard Headers */
#include<math.h>
#include<stdio.h>
int main() {
    /* Variable declarations */
    int i,j;
    unsigned int n,m;
    unsigned int rows,cols;
    unsigned int D[3][3], d[3], C[3][3];
    unsigned int x, s[3][3];
    unsigned int y[3];
    void LU();    
    n = 2;
    rows=3;cols=3;

    /* the matrix to be inverted */

    D[0][0] = 1;
    D[0][1] = 1;
    D[0][2] = 0;
    D[1][0] = 0;
    D[1][1] = 1;
    D[1][2] = 0;
    D[2][0] = 1;
    D[2][1] = 1;
    D[2][2] = 1;

    /* Store the matrix value for camparison later.
     this is just to check the results, we don't need this
     array for the program to work */
    for (m = 0; m <= rows-1; m++) {
        for (j = 0; j <= cols-1; j++) {
            C[m][j] = D[m][j];
        }
    }

    /* Call a sub-function to calculate the LU decomposed matrix. Note that
     we pass the two dimensional array [D] to the function and get it back */
    LU(D, n);

    printf(" \n");
    printf("The matrix LU decomposed \n");
    for (m = 0; m <= rows-1; m++) {
        for (j = 0; j <= cols-1; j++){
        printf(" %d \t", D[m][j]);
        }
        printf("\n");
    }

    /*  TO FIND THE INVERSE */

    /* to find the inverse we solve [D][y]=[d] with only one element in
     the [d] array put equal to one at a time */

    for (m = 0; m <= rows-1; m++) {
        d[0] = 0;
        d[1] = 0;
        d[2] = 0;
        d[m] = 1;
        for (i = 0; i <= n; i++) {
            x = 0;
            for (j = 0; j <= i - 1; j++){
                x = x + D[i][j] * y[j];
            }
            y[i] = (d[i] - x);
        }

        for (i = n; i >= 0; i--) {
            x = 0;
            for (j = i + 1; j <= n; j++) {
                x = x + D[i][j] * s[j][m];
            }
            s[i][m] = (y[i] - x) / D[i][i];
        }
    }

    /* Print the inverse matrix */
    printf("The Inverse Matrix\n");
    for (m = 0; m <= rows-1; m++) {
        for (j = 0; j <= cols-1; j++){
            printf(" %d \t", s[m][j]);
        }
        printf("\n");
    }
    /* check that  the product of the matrix with its iverse results
     is  indeed a unit matrix  */
    printf("The product\n");
    for (m = 0; m <= rows-1; m++) {
        for (j = 0; j <= cols-1; j++){
            x = 0;
            for (i = 0; i <= 2; i++) {
                x = x + C[m][i] * s[i][j];
            }
            //printf(" %d    %d    %f \n", m, j, x);
            printf("%d \t",x);
        }
        printf("\n");
    }
    return 0;
}

/* The function that calcualtes the LU deomposed matrix.
 Note that it receives the matrix as a two dimensional array
 of pointers. Any change made to [D] here will also change its
 value in the main function. So there is no need of an explicit
 "return" statement and the function is of type "void". */

void LU(int (*D)[3][3], int n) {
    int i, j, k;
    int x;
    printf("The matrix \n");
    for (j = 0; j <= 2; j++) {
        printf(" %d  %d  %d \n", (*D)[j][0], (*D)[j][1], (*D)[j][2]);
    }
    for (k = 0; k <= n - 1; k++) {
        for (j = k + 1; j <= n; j++) {
            x = (*D)[j][k] / (*D)[k][k];
            for (i = k; i <= n; i++) {
                (*D)[j][i] = (*D)[j][i] - x * (*D)[k][i];
            }
            (*D)[j][k] = x;
        }
    }

}

This is just a sample example that I tried and I have -1 values in the inverse matrix which is my main concern. I have 1000 x 1000 matrix of binary values and the inverse should also be in binary.

The matrix:

 1  1  0 
 0  1  0 
 1  1  1

The matrix LU decomposed:

 1   1   0  
 0   1   0  
 1   0   1  

The Inverse Matrix:

 1  -1   0  
 0   1   0  
-1   0   1

The product:

 1   0   0  
 0   1   0  
 0   0   1