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
