线性方程组 Ax=b 求解(Chapter28)

1. /**
2. *线性方程组 Ax=b 求解函数
3. *@param     A       (in)系数矩阵
4. *@param     x       (out)根列向量，接受输出值
5. *@param     b       (in)方程右侧常数项
6. *@param     size    (in)矩阵的大小
7. *@return    求解成功返回非零值，失败返回零
8. *@LUDecomposition LU分解-详见上一篇或下一篇                            
9. */
10. int LinearEquation(double **A, double *x, double *b, int size)
11. {
12.     double **L, **U, *y;
13.     double temp;
14.     int i, j;
16.     L = (double**)malloc(sizeof(double*)*size);
17.     U = (double**)malloc(sizeof(double*)*size);
19.     for (i = 0; i < size; i++)
20.     {
21.         L[i] = (double*)malloc(sizeof(double)*size);
22.         U[i] = (double*)malloc(sizeof(double)*size);
23.     }
25.     if (!LUDecomposition(A, L, U, size))
26.     {
27.         fprintf(stderr, "Error: In LinearEquation./n");
28.         return 0;
29.     }
31.     y = (double*)malloc(sizeof(double)*size);
33.     for (i = 0; i < size; i++)
34.     {
35.         temp = 0;
36.         y[i] = 0;
38.         for (j = 0; j <= i-1; j++)
39.             temp += L[i][j] * y[j];
41.         y[i] += b[i] - temp;
42.     }
44.     for (i = size-1; i >= 0; i--)
45.     {
46.         temp = 0;
47.         x[i] = 0;
48.         for (j = i+1; j < size; j++)
49.             temp += U[i][j] * x[j];
51.         x[i] = (y[i] - temp) / U[i][i];
52.     }
54.     for (i = 0; i < size; i++)
55.     {
56.         free(L[i]);
57.         free(U[i]);
58.     }
60.     free(L);
61.     free(U);
62.     free(y);
64.     return 1;
65. }