使用 LU 分解给出 inf 值求解 A.x = b

Question

我正在尝试使用 A = LU 算法求解 Ax = b 形式的线性系统，其中 A 是 (nxn) 实数矩阵，b 是 (1xn) 实数向量。这是我的实现：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int LUPDecompose(double A[N][N], double Tol, int P[N])
{

    int i, j, k, imax;
    double maxA, ptr[N], absA;

    for (i = 0; i <= N; i++)
        P[i] = i; //Unit permutation matrix, P[N] initialized with N

    for (i = 0; i < N; i++) {
        maxA = 0.0;
        imax = i;

        for (k = i; k < N; k++)
            if ((absA = abs(A[k][i])) > maxA) {
                maxA = absA;
                imax = k;
            }

        if (maxA < Tol) return 0; //failure, matrix is degenerate

        if (imax != i) {
            //pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            //pivoting rows of A
            for (int ii = 0; ii < N; ii++)
            {
                ptr[ii] = A[i][ii];
                A[i][ii] = A[imax][ii];
                A[imax][ii] = ptr[ii];
            }

            //counting pivots starting from N (for determinant)
            P[N]++;
        }

        for (j = i + 1; j < N; j++) {
            A[j][i] /= A[i][i];

            for (k = i + 1; k < N; k++)
                A[j][k] -= A[j][i] * A[i][k];
        }
    }

    return 1;  //decomposition done
}

/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
 * OUTPUT: x - solution vector of A*x=b
 */
void LUPSolve(double A[N][N], int P[N], double b[N], double x[N])
{

    for (int i = 0; i < N; i++) {
        x[i] = b[P[i]];

        for (int k = 0; k < i; k++)
            x[i] -= A[i][k] * x[k];
    }

    for (int i = N - 1; i >= 0; i--) {
        for (int k = i + 1; k < N; k++)
            x[i] -= A[i][k] * x[k];

        x[i] /= A[i][i];
    }
}

int main()
{
    double Am[N][N] = {{0.6289, 0, 0.0128, 0.3184, 0.7151},
                        {0,   1,   0,   0,   0},
                        {0.0128, 0, 0.0021, 0.0045, 0.0380},
                        {0.3184, 0, 0.0045, 0.6618, 0.3371},
                        {0.7151, 0, 0.0380, 0.3371, 1.1381}};
    double bm[N] = {1.6752, 0, 0.0574, 1.3217, 2.2283};
    int Pm[N] = {0};
    double X[N] = {0};

    LUPDecompose( Am, 0.0001, Pm);
    LUPSolve(Am, Pm,  bm, X);

    printf("%f %f %f %f %f",X[0],X[1],X[2],X[3],X[4]);
}

但是，我得到的是 inf 值。

-1.#IND00 -1.#IND00 3.166387 0.849298 0.670689

请问是代码问题还是算法问题。对解决这个问题有帮助吗？

Answer 1

"I wonder if it is a code issue or algorithm. Any help to solve this issue?"

我认为存在代码和算法问题。以下是您的代码，其中修正了 仅编译错误和警告 （请参阅行内注释）。它不调试超出 C 语法以实现干净的编译，并且运行 w/o 错误。（即运行没有 除以零 或 inf 错误。）

#define N 5 //required to be 5 by hard-coded array definitions in main()

int LUPDecompose(double A[N][N], double Tol, int P[N])
{

    int i, j, k, imax, ii;//added ii here to increase scope below
    double maxA, ptr[N], absA;

    //for (i = 0; i <= N; i++)
    for (i = 0; i < N; i++)
        P[i] = i; //Unit permutation matrix, P[N] initialized with N (actually init with i)

    for (i = 0; i < N; i++) {
        maxA = 0.0;
        imax = i;

        for (k = i; k < N; k++)
            if ((absA = fabs(A[k][i])) > maxA) {// using fabs, not abs to avoid conversion of double to int.
                maxA = absA;
                imax = k;
            }

        if (maxA < Tol) return 0; //failure, matrix is degenerate

        if (imax != i) {
            //pivoting P
            j = P[i];
            P[i] = P[imax];
            P[imax] = j;

            //pivoting rows of A
           //for (int ii = 0; ii < N; ii++)
           for ( ii = 0; ii < N; ii++)
            {
                ptr[ii] = A[i][ii];
                A[i][ii] = A[imax][ii];
                A[imax][ii] = ptr[ii];
            }

            //counting pivots starting from N (for determinant)
            //P[N]++;//N will always overflow for array with only N elements
            P[ii-1]++;//use index here instead
        }

        for (j = i + 1; j < N; j++) {
            A[j][i] /= A[i][i];

            for (k = i + 1; k < N; k++) {//extra brackets added for readability
                A[j][k] -= A[j][i] * A[i][k];
            }
        }
    }

    return 1;  //decomposition done
}

/* INPUT: A,P filled in LUPDecompose; b - rhs vector; N - dimension
 * OUTPUT: x - solution vector of A*x=b
 */
void LUPSolve(double A[N][N], int P[N], double b[N], double x[N])
{

    for (int i = 0; i < N; i++) {
        x[i] = b[P[i]];
        for (int k = 0; k < i; k++) {//extra brackets added for readability
            x[i] -= A[i][k] * x[k];
        }
    }

        for (int i = N - 1; i >= 0; i--) {
            for (int k = i + 1; k < N; k++) {//additional brackets added for readability
                x[i] -= A[i][k] * x[k];
            }

            x[i] /= A[i][i];
        }
}
    

//int main()
int main(void)//minimum signature for main includes void
{
    //Note hardcoded arrays in this code require N == 5 (#define at top)
    double Am[N][N] = {{0.6289, 0, 0.0128, 0.3184, 0.7151},
                        {0,   1,   0,   0,   0},
                        {0.0128, 0, 0.0021, 0.0045, 0.0380},
                        {0.3184, 0, 0.0045, 0.6618, 0.3371},
                        {0.7151, 0, 0.0380, 0.3371, 1.1381}};
    double bm[N] = {1.6752, 0, 0.0574, 1.3217, 2.2283};
    int Pm[N] = {0};
    double X[N] = {0};

    LUPDecompose( Am, 0.0001, Pm);
    LUPSolve(Am, Pm,  bm, X);

    printf("%f %f %f %f %f",X[0],X[1],X[2],X[3],X[4]);
    
    return 0;  //int main(void){...} requires return statement.
}

基于this calculator, with these inputs:

正确的答案是：

-0.590174531351002
0
-19.76923076923077
1.0517711171662125
2.6772727272727272

但上面代码的实际输出是：

算法相关调试留给你去执行。

使用 LU 分解给出 inf 值求解 A.x = b

Solve A.x = b using LU factorisation give inf values

c

solver