Инверсия двоичной матрицы в C

Question

У меня есть двоичная матрица (нули и единицы). D[][] размера nxn, где n большое (примерно 1500 - 2000). Я хочу найти обратную матрицу в C.

Поскольку я новичок в C, я начал с матрицы 3 x 3 и начал работать над обобщением ее на N x N. Это работает для значений int, однако, поскольку я работаю с двоичными 1 и 0 s. В этой реализации мне нужны unsigned int значения.

Я мог бы найти много решений для int значений, но я не нашел ни одного решения для unsigned int. Я хотел бы найти инверсию двоичной матрицы N x N без использования каких-либо внешних библиотек, таких как blas / lapack. Было бы замечательно, если бы кто-нибудь мог обеспечить вывод на матрицу M x N.

Обратите внимание, что мне нужна обратная матрица, а не псевдообратная.

/* 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;
        }
    }

}

Это просто пример, который я попробовал, и у меня есть -1 значения в обратной матрице, что является моей главной задачей. У меня есть матрица двоичных значений 1000 x 1000, и обратное значение также должно быть в двоичном виде.

Матрица:

 1  1  0 
 0  1  0 
 1  1  1

Матрица LU разложена:

 1   1   0  
 0   1   0  
 1   0   1  

Обратная матрица:

 1  -1   0  
 0   1   0  
-1   0   1

Продукт:

 1   0   0  
 0   1   0  
 0   0   1