Приведение матрицы к диагональному виду методом Гаусса - C (СИ)

Формулировка задачи:

Здравствуйте. Написал программу,которая расскладывает любой файл по битам и составляет из них матрицу заданных размеров. Задача - привести эту матрицу к диагональному виду методом Гаусса. Помогите пожалуйста с алгоритмом. Вот код программы:

#include <stdio.h>
#include <stdlib.h>

#define E_ALLOC "\nCan't allocate memory\n"
#define BIT_ERROR 10

typedef struct T_BYTE{		//Байт, как последовательность бит
	unsigned bit7: 1;
	unsigned bit6: 1;
	unsigned bit5: 1;
	unsigned bit4: 1;
	unsigned bit3: 1;
	unsigned bit2: 1;
	unsigned bit1: 1;
	unsigned bit0: 1;
}T_BYTE;

T_BYTE readbyte(FILE*);
unsigned get_bit(T_BYTE byte,unsigned n); //Получить n-ый бит из байта byte

unsigned **matrix_malloc(size_t,size_t);
void matrix_free(unsigned**,size_t);
void print_matrix(unsigned**,size_t,size_t);

int main(int argc,char *argv[]){
	if(argc==1){
		printf("No input files\n");
		exit(0);
	}
	
	unsigned **matrix;
	size_t m,n; //размер матрицы
	FILE *f;

	f=fopen(argv[1],"r");

        if(f==NULL){
                printf("File \"%s\" doesn't exist",argv[1]);
                matrix_free(matrix,m);
                exit(0);
        }

	printf("Input size of matrix: ");
	scanf("%u%u",&m,&n);
	matrix=matrix_malloc(m,n);

	if(matrix==NULL){
		perror(E_ALLOC);
		exit(0);
	}

	size_t i,j,pos=0;
	T_BYTE byte=readbyte(f);
	
	/*Заполнение матрицы данными из файла.
	  pos хранит номер бита из байта byte
	  и с помощью get_bit мы получаем значение этого байта  
	  если pos==8, то считываем из файла следующий байт
	*/	

	for(i=0;i<m;++i){
		for(j=0;j<n;++j){
			if(pos>7){
				byte=readbyte(f); 
				pos=0;
			}
			matrix[i][j]=get_bit(byte,pos++);
		}
	}

	print_matrix(matrix,m,n);
		
	matrix_free(matrix,m);
	fclose(f);
	return 0;
}

unsigned **matrix_malloc(size_t m,size_t n){
	unsigned **matrix;
	matrix=(unsigned**)malloc(sizeof(unsigned*)*m);
	
	if(matrix==NULL)
		return NULL;

	size_t i;
	for(i=0;i<m;++i){
		matrix[i]=(unsigned*)malloc(sizeof(unsigned)*n);

		if(matrix[i]==NULL){
			if(i){
				size_t j;
				for(j=0;j<i;++j)
					free(matrix[j]);
			}		
			free(matrix);
			return NULL;
		}				

	}
	return matrix;	
}

void matrix_free(unsigned **matrix,size_t m){
	size_t i;
	for(i=0;i<m;++i)
		free(matrix[i]);
	free(matrix);	
}

void print_matrix(unsigned **matrix,size_t m,size_t n){
	size_t i,j;
	for(i=0;i<m;++i){
		for(j=0;j<n;++j)
			printf("%u ",*(*(matrix+i)+j));
		printf("\n");
	}
}

T_BYTE readbyte(FILE *f){
	union{
		T_BYTE byte;
		char ch;
	}buf;
	if(feof(f))
		buf.ch=0;
	else
		buf.ch=getc(f);
	return buf.byte;
}

unsigned get_bit(T_BYTE byte,unsigned n){
	switch(n){
		case 0:
			return byte.bit0;
		case 1:
			return byte.bit1;
		case 2:
			return byte.bit2;
		case 3:
			return byte.bit3;
		case 4:
			return byte.bit4;
		case 5:
			return byte.bit5;
		case 6:
			return byte.bit6;
		case 7: 
			return byte.bit7;
		default:
			return BIT_ERROR;	
	}
}

Решение задачи: «Приведение матрицы к диагональному виду методом Гаусса»

Листинг программы

    const int n = 5;
    double A[n][n] , b[n];
 
    // здесь должна быть инициализация матрицы A и матрицы свободных коэффициентов b
 
    for (int k = 0; k < n; k ++)
    {
        for (int i = k + 1; i < n; i ++)
        {
            double mu = A[i][k] / A[k][k];
            for (int j = 0; j < n; j ++)
                A[i][j] -= A[k][j] * mu;
            b[i] -= b[k] * mu;
        }
    }

Объяснение кода листинга программы

В коде объявлены две матрицы: A[n][n] и b[n].
Матрица A представляет собой квадратную матрицу размером n на n, которая, вероятно, содержит коэффициенты некоторого квадратного уравнения.
Матрица b представляет собой вектор свободных членов этого квадратного уравнения.
Код реализует алгоритм приведения матрицы к диагональному виду методом Гаусса.
Внешний цикл for (int k = 0; k < n; k ++) инициализирует главный элемент на диагонали матрицы A[k][k] и выполняет все необходимые операции для приведения матрицы к диагональному виду.
Внутренний цикл for (int i = k + 1; i < n; i ++) перебирает строки матрицы A, начиная с позиции (k+1,k+1) и до конца строки.
В каждой итерации внутреннего цикла вычисляется значение mu = A[i][k] / A[k][k], которое представляет собой соответствующий элемент главной диагонали матрицы.
Затем, используя это значение, все элементы строки i матрицы A вычитаются на mu, умноженное на соответствующий элемент строки k.
Аналогично, все элементы вектора b, начиная с позиции i, вычитаются на mu, умноженное на соответствующий элемент вектора b, находящегося на позиции k.
Код не содержит никаких проверок на деление на ноль, поэтому предполагается, что матрица A не содержит нулевых элементов.
Код не содержит никаких комментариев или пояснений, поэтому трудно сказать, что именно он должен делать или какие результаты он должен давать.