Doolittle分解

#include <stdio.h>
#include <math.h>

void Doolittle(int n, double A[][100], double L[][100], double U[][100]);
double compute_det(int n, double A[][100]);

int main() {
 	int n;
  	double A[100][100], L[100][100], U[100][100];
  	
  	scanf("%d", &n);
	for(int i = 0; i < n; i++){
		for(int j = 0; j < n; j++){
			scanf("%lf", &A[i][j]);
		}
	}
	Doolittle(n, A, L, U);
//	printf("\n");
//	printf("%lf", compute_det(n-1, A));
	for(int i = 0; i < n; i++){
		for(int j = 0;j < n; j++){
			printf("%lf ", L[i][j]);	
		}	
		printf("\n");
	}
	for(int i = 0; i < n; i++){
		for(int j = 0;j < n; j++){
			printf("%lf ", U[i][j]);	
		}
		printf("\n");
	}
}


double compute_det(int n, double A[][100]) {
	double ret = 0;
	double B[100][100];
	if (n == 1) {
		ret = A[0][0];
	} else {
		for(int i=0; i<n; i++) {
			for (int j=0; j<i; j ++) {
				for (int k=1; k<n; k ++) {
					B[j][k-1] = A[j][k]; 
				}
			}
			for (int j=i; j<n; j ++) {
				for (int k=1; k<n; k ++) {
					B[j][k-1] = A[j+1][k]; 
				}
			}
			ret += A[i][0] * pow(-1, i+0) * compute_det(n-1, B);
		}
	}
	return ret;
}

void Doolittle(int n, double A[][100], double L[][100], double U[][100]){
  	if (compute_det(n-1, A) == 0) {
		return;
	}
	for (int k=0; k<n; k++) {
		L[k][k] = 1;
		for (int i=k; i<n; i++) {
			double sum=0;
			for (int j=0; j<k; j++) {
				sum += L[k][j] * U[j][i];
			}
			U[k][i] = A[k][i] - sum;
		}
		for (int i_=k+1; i_<n; i_++) {
			double sum_=0;
			for (int j_=0; j_<k; j_++) {
				sum_ += L[i_][j_] * U[j_][k];
			}
			L[i_][k] = (A[i_][k] - sum_) / U[k][k];
		}
	}
}