我在其中一个循环中使用了错误的索引，但不知道是哪一个。我已经进行了建议的更改

Question

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




const int N = 3;

    
void LUBKSB(double b[], double a[N][N], int N, int *indx)
{
    int i, ii, ip, j;
    double sum;
    ii = 0;
    
    for(i=0;i<N;i++)
    {
        ip = indx[i];
        sum = b[ip];
        b[ip] = b[i];
        
        if (ii)
        {
            for(j = ii;j<i-1;j++)
            {
                sum = sum - a[i][j] * b[j];
            }
        }
            else if(sum)
            {
                ii = i;
            }
            b[i] = sum;
    }
    for(i=N-1;i>=0;i--)
    {
        sum = b[i];
        for (j = i; j<N;j++)
        {
            sum = sum - a[i][j] * b[j];
        }
        b[i] = sum/a[i][i];
    }
    
    for (i=0;i<N;i++)
    {
        printf("b[%d]: %lf \n",i,b[i]);
        
    }
}

void ludecmp(double a[][3], int N)
{
    int i, imax, j, k;
    double big, dum, sum, temp, d;
    double *vv = (double *) malloc(N * sizeof(double));
    int *indx = (int *) malloc(N * sizeof(double));
    double TINY = 0.000000001;
    
    double b[3] = {2*M_PI,5*M_PI,-8*M_PI};
    d = 1.0;
    
    for(i=0;i<N;i++)
    {
        big = 0.0;
        for(j=0;j<N;j++)
        {
            temp = fabs(a[i][j]);
            if (temp > big)
            {
                big = temp;
            }
        }
            if (big == 0.0)
            {
                printf("Singular matrix\n");
                exit(1);
            }
        vv[i] = 1.0/big;
    }
    for(j=0;j<N;j++)
    {
        for(i=0;i<j-1;i++)
        {
            sum = a[i][j];
            for(int k=0;k<i-1;k++)
            {
                sum = sum - (a[i][k] * a[k][j]);
            }
            a[i][j] = sum;
        }
            big = 0.0;
            for(i=j;i<N;i++)
            {
                sum = a[i][j];
                for(k=0;k<j-1;k++)
                {
                    sum = sum - a[i][k] * a[k][j];
                }
                a[i][j] =sum;
                dum = vv[i] * fabs(a[i][j]);
                if(dum >= big)
                {
                    big = dum;
                    imax = i;
                }
            }
            if(j != imax)
            {
                for(k=0;k<N;k++)
                {
                    dum = a[imax][k];
                    a[imax][k] = a[j][k];
                    a[j][k] = dum;
                }
                d = -d;
                vv[imax] = vv[j];
            }
            indx[j] = imax;
            if (a[j][j] == 0)
            {
                a[j][j] = TINY;
            }
            if (j != N)
            {
                dum = 1.0/a[j][j];
                for(i = j; i<N; i++)
                {
                    a[i][j] = a[i][j] * dum;
                }
            }
    }
    
    LUBKSB(b,a,N,indx);
    
    free(vv);
    free(indx);
   
}
    
int main()
{
    
    int N, i, j;

    
    N = 3;
    double a[3][3] = { 1, 2, -1, 6, -5, 4, -9, 8, -7};
    
    ludecmp(a,N);
    
    
}

    
    

我正在使用这些算法来查找矩阵的 LU 分解并试图找到解 A.x = b

给定一个 N ×N 矩阵 A，表示为 {a}N,Ni,j=1，例程将其替换为 LU 分解自身的行排列。输入“a”和“N”。 “a”也输出， 修改为应用 LU 分解； {索引}N i=1 是一个输出向量，记录了 由部分旋转实现的行排列； “d”是输出，采用±1取决于 行交换的数量是偶数还是奇数。此例程结合使用 使用算法 2 求解线性方程或求矩阵求逆。

求解 N 个线性方程组 A 。 x = b。矩阵{a} N,N i,j=1 实际上是 从算法 1 获得的原始矩阵 A 的 LU 分解。向量 {indxi} ñ i=1 是 输入作为算法 1 返回的置换向量。向量 {bi} ñ i=1 作为右手边向量 B 输入，但返回解向量 X。输入 {a} N,N i,j=1, N 和 {indxi} ñ 我=1 在这个算法中没有被修改。

Answer 1

你的代码有很多问题：

1. 在你的 for 循环中，i <= N 应该是 i < N，i = N 应该是 i = N - 1。
2. double 的绝对值由 fabs 返回，而不是 abs。
3. 声明exit 应该是exit(1) 或exit(EXIT_FALILURE)。
4. 您的两个函数缺少 return 语句。

您还应该释放使用函数free 分配的内存。编译 C 程序时，您还应该启用所有警告。