csr格式存储的矩阵a利用选主元lu分解求解ax=b的c代码

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

void lu_decomposition(double** a, int n, int* p) {
    int i, j, k, row;
    double max, temp;
    for (i = 0; i < n; i++) {
        p[i] = i;
    }
    for (i = 0; i < n; i++) {
        max = 0.0;
        row = i;
        for (j = i; j < n; j++) {
            if (abs(a[j][i]) > max) {
                max = abs(a[j][i]);
                row = j;
            }
        }
        if (max == 0.0) {
            printf("Singular matrix\n");
            exit(EXIT_FAILURE);
        }
        if (row != i) {
            for (j = 0; j < n; j++) {
                temp = a[i][j];
                a[i][j] = a[row][j];
                a[row][j] = temp;
            }
            k = p[i];
            p[i] = p[row];
            p[row] = k;
        }
        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];
            }
        }
    }
}

void lu_solve(double** a, double* b, int* p, int n, double* x) {
    int i, j;
    double sum;
    for (i = 0; i < n; i++) {
        x[i] = b[p[i]];
    }
    for (i = 0; i < n; i++) {
        sum = x[i];
        for (j = 0; j < i; j++) {
            sum -= a[i][j] * x[j];
        }
        x[i] = sum;
    }
    for (i = n - 1; i >= 0; i--) {
        sum = x[i];
        for (j = i + 1; j < n; j++) {
            sum -= a[i][j] * x[j];
        }
        x[i] = sum / a[i][i];
    }
}

int main() {
    int n = 3;
    double** a = (double**)malloc(n * sizeof(double*));
    int i, j, k;
    for (i = 0; i < n; i++) {
        a[i] = (double*)malloc(n * sizeof(double));
    }
    int* p = (int*)malloc(n * sizeof(int));
    double* b = (double*)malloc(n * sizeof(double));
    double* x = (double*)malloc(n * sizeof(double));

    a[0][0] = 2.0;
    a[0][1] = 1.0;
    a[0][2] = 1.0;
    a[1][0] = 4.0;
    a[1][1] = 3.0;
    a[1][2] = 3.0;
    a[2][0] = 8.0;
    a[2][1] = 7.0;
    a[2][2] = 9.0;
    b[0] = 3.0;
    b[1] = 7.0;
    b[2] = 16.0;

    lu_decomposition(a, n, p);
    lu_solve(a, b, p, n, x);

    printf("Solution:\n");
    for (i = 0; i < n; i++) {
        printf("%lf\n", x[i]);
    }

    for (i = 0; i < n; i++) {
        free(a[i]);
    }
    free(a);
    free(p);
    free(b);
    free(x);

    return 0;
}