紧凑存储的杜利特尔分解法Doolittle(LU分解法)_解线性方程组的直接解法

标签：计算方法实验

/*紧凑存储的杜利特尔分解法Doolittle：    如果初始矩阵不要求保留的话，可以紧凑存储。    因为每个a[i][j]用来计算u[i][j]或l[i][j]之后不需要了，u[i][j]或l[i][j]算出后可存入a[i][j]所占的单元；    y[i]可存入b[i]所占的单元，最后x[i]再取代y[i]存入最初b[i]所占的单元。    注意用 '' 标记部分。*/#include 
   
    #include 
    
     const int maxn = 15;int main(){    double a[maxn][maxn], b[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[i][1] /= temp;  //L的第一列元素, U的第一行元素不变    for(k = 2; k <= n; k++){        for(j = k; j <= n; j++){  //计算U的第k行元素            for(r = 1, sum = 0; r <= k - 1; r++)  sum += (a[k][r] * a[r][j]);            a[k][j] -= sum;        }        for(i = k; i <= n; i++){  //计算L的第k列元素            for(r = 1, sum = 0; r <= k - 1; r++)  sum += (a[i][r] * a[r][k]);            if(i == k)  continue;  对角线在计算U时已被计算出来            else  a[i][k] = (a[i][k] - sum) / a[k][k];        }    }    /*打印L   U, L左下, U右上    for(i = 1; i <= n; i++){        for(j = 1; j <= n; j++)  printf("%10f", a[i][j]);        printf("\n");    }    */    for(i = 2; i <= n; i++){  //求解Ly = b, y[1] = b[1]        for(k = 1, sum = 0; k <= i - 1; k++)  sum += (a[i][k] * b[k]);        b[i] -= sum;    }    b[n] /= a[n][n];  //求解Ux = y    for(i = n - 1; i >= 1; i--){        for(k = i + 1, sum = 0; k <= n; k++)  sum += (a[i][k] * b[k]);        b[i] = (b[i] - sum) / a[i][i];    }    for(int i = 1; i <= n; i++)  printf("%10f\n", b[i]);    return 0;}

数据文件