满秩矩阵的LU分解-非递归法(Chapter 28)

1. /**
2. *LU分解函数(V2非递归法)
3. *@param     A       (in)输入矩阵
4. *@param     L       (out)输出L阵
5. *@param     U       (out)输出U阵
6. *@param     size    (in)矩阵的大小
7. *@return    分解成功返回非零值，失败返回零
8. */
9. int LUDecomposition(double **A, double **L, double **U, int size)
10. {
11.     int i, j, k;
13.     double **Ax;
15.     Ax = (double**)malloc(sizeof(double*)*size);
17.     for (i = 0; i < size; i++)
18.         Ax[i] = (double*)malloc(sizeof(double)*size);
20.     for (i = 0; i < size; i++)
21.         for (j = 0; j < size; j++)
22.             Ax[i][j] = A[i][j];
24.     for (i = 1; i < size; i++)
25.         for (j = 0; j < i; j++)
26.             U[i][j] = 0;
28.     for (i = 0; i < size; i++)
29.     {
30.         L[i][i] = 1;
31.         for (j = i+1; j < size; j++)
32.             L[i][j] = 0;
33.     }
35.     for (k = 0; k < size; k++)
36.     {
37.         U[k][k] = Ax[k][k];
38.         for (i = k+1; i < size; i++)
39.         {
40.             L[i][k] = Ax[i][k] / U[k][k];
41.             U[k][i] = Ax[k][i];
42.         }
43.         for (i = k+1; i < size; i++)
44.             for (j = k+1; j < size; j++)
45.                 Ax[i][j] = Ax[i][j] - L[i][k] * U[k][j];
46.     }
48.     for (i = 0; i < size; i++)
49.         free(Ax[i]);
50.     free(Ax);
52.     return 1;
53. }