使用LU分解求解Ax=b时输出inf值，排查代码或算法问题

你的代码同时存在多个语法和逻辑错误，和LUP分解算法本身无关，修复后即可正常运行：

• 浮点数绝对值计算错误：你使用abs()函数计算双精度浮点数的绝对值，abs()是针对整型的库函数，传入的浮点数会被强制截断为整数，所有小于1的正数都会被转为0，直接导致选主元逻辑失效，误判矩阵为奇异矩阵。后续计算中出现除0操作，生成NaN（即你输出中的-1.#IND00），这是导致输出异常的核心原因，需要替换为处理浮点数的fabs()函数。
• 数组越界访问：你定义的置换数组P长度为N，合法索引范围是0~N-1，但在LUPDecompose函数中你用for (i = 0; i <= N; i)初始化P，还写入了P[N]++统计置换次数，越界修改了其他变量的内存空间，会破坏矩阵、向量的原始值，导致计算错误。如果需要统计置换次数，可以单独新增参数传递，不要越界访问P数组。
• 缺失错误处理逻辑：LUPDecompose函数在矩阵奇异时会返回0表示分解失败，你在主函数调用时没有接收判断返回值，分解失败后依然执行求解逻辑，必然得到错误结果。
• 常量N未定义：你的代码中直接使用常量N，既没有宏定义也没有作为函数参数传递，本身无法通过编译，你测试用的是5阶矩阵，需要在代码开头新增#define N 5的宏定义。

1. 在代码开头添加N的宏定义：

#define N 5
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

2. 将LUPDecompose中选主元的abs(A[k][i])修改为fabs(A[k][i])
3. 修改P数组的初始化逻辑，删除越界访问代码：

// 原初始化代码
// for (i = 0; i <= N; i++)
//    P[i] = i;
// 修改为
for (i = 0; i < N; i++)
    P[i] = i;
// 删除P[N]++;这一行

4. 主函数新增分解结果判断：

if (!LUPDecompose( Am, 0.0001, Pm)) {
    printf("矩阵奇异，分解失败\n");
    return 1;
}
LUPSolve(Am, Pm,  bm, X);

修复后运行即可得到正确的求解结果。

内容的提问来源于stack exchange，提问作者nechi