#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];
}
}
}
Doolittle分解
于 2023-04-01 21:07:54 首次发布