本文共 1843 字,大约阅读时间需要 6 分钟。
/*紧凑存储的克洛特分解法Crout: 如果初始矩阵不要求保留的话,可以紧凑存储。 因为每个a[i][j]用来计算u[i][j]或l[i][j]之后不需要了,u[i][j]或l[i][j]算出后可存入a[i][j]所占的单元; x[i]可存入y[i]所占的单元,但y[i]不存入b[i]所占的单元,因为y[i]的计算需要旧的b[i]值。 注意用 '' 标记部分。*/#include#include const int maxn = 15;int main(){ double a[maxn][maxn], b[maxn], y[maxn]; int i, j, k, r, n, sum, temp; freopen("lu.txt", "r", stdin); scanf("%d", &n); for(int i = 1; i <= n; i++){ for(int j = 1; j <= n; j++) scanf("%lf", &a[i][j]); scanf("%lf", &b[i]); } /*打印数据文件 for(int i = 1; i <= n; i++){ for(int j = 1; j <= n; j++) printf("%10f", a[i][j]); printf("%10f\n", b[i]); } */ temp = a[1][1]; //////// for(i = 1; i <= n; i++) a[1][i] /= temp; //U的第一行元素, L的第一列元素不变 for(k = 2; k <= n; k++){ for(i = k; i <= n; i++){ //计算L的第k列元素 for(r = 1, sum = 0; r <= k - 1; r++) sum += (a[i][r] * a[r][k]); a[i][k] -= sum; } for(j = k; j <= n; j++){ //计算U的第k行元素 for(r = 1, sum = 0; r <= k - 1; r++) sum += (a[k][r] * a[r][j]); if(k == j) continue; ////////对角线元素在计算L时已被计算出来 else a[k][j] = (a[k][j] - sum) / a[k][k]; } } /*打印L U for(i = 1; i <= n; i++){ for(j = 1; j <= n; j++) printf("%10f", l[i][j]); printf("\t\t"); for(k = 1; k <= n; k++) printf("%10f", u[i][k]); printf("\n"); } */ y[1] = b[1] / a[1][1]; //求解Ly = b for(i = 2; i <= n; i++){ for(k = 1, sum = 0; k <= i; k++) sum += (a[i][k] * y[k]); y[i] = (b[i] - sum) / a[i][i]; } for(i = n - 1; i >= 1; i--){ //求解Ux = y for(k = i + 1, sum = 0; k <= n; k++) sum += (a[i][k] * y[k]); y[i] -= sum; } for(int i = 1; i <= n; i++) printf("%10f\n", y[i]); return 0;}
数据文件
实验结果