#include <stdio.h>  
#include <stdlib.h>  
#include <math.h>  
#include <time.h>  
  
// 稀疏矩阵数据结构  
typedef struct _triplet TRIPLET;  
struct _triplet  
{  
    int i;  
    int j;  
    double value;  
};  
  
typedef struct _s S;  
struct _s  
{  
    int nRow;  
    int nCol;  
    int nz;  
    TRIPLET* ts;  
};  
  
S* sCreate(int nRow, int nCol, int nz)  
{  
    S* s = (S*)malloc(sizeof(S));  
    s->nRow = nRow;  
    s->nCol = nCol;  
    s->nz = nz;  
    s->ts = (TRIPLET*)malloc(nz * sizeof(TRIPLET));  
    for(int i = 0; i < nz; i++)  
    {  
        s->ts[i].i = 0;  
        s->ts[i].j = 0;  
        s->ts[i].value = 0;  
    }  
    return s;  
}  
  
void sFree(S* s)  
{  
    free(s->ts);  
    free(s);  
}  
  
// 雅克比迭代法求解稀疏线性方程组  
int jacobi_sparse(S* A, double* b, double* x, int max_iter, double tolerance)  
{  
    int n = A->nRow;  
    double* x_new = (double*)malloc(n * sizeof(double));  
    double* diag = (double*)malloc(n * sizeof(double));  
  
    // 提取对角线元素  
    for(int i = 0; i < n; i++) {  
        diag[i] = 1.0; // 初始化为1，后面会更新  
    }  
  
    for(int k = 0; k < A->nz; k++) {  
        if(A->ts[k].i == A->ts[k].j) {  
            diag[A->ts[k].i] = A->ts[k].value;  
        }  
    }  
  
    int iter;  
    for(iter = 0; iter < max_iter; iter++) {  
        // 初始化新解向量  
        for(int i = 0; i < n; i++) {  
            x_new[i] = b[i];  
        }  
  
        // 计算非对角线元素的贡献  
        for(int k = 0; k < A->nz; k++) {  
            int i = A->ts[k].i;  
            int j = A->ts[k].j;  
            if(i != j) {  
                x_new[i] -= A->ts[k].value * x[j];  
            }  
        }  
  
        // 除以对角线元素  
        for(int i = 0; i < n; i++) {  
            x_new[i] /= diag[i];  
        }  
  
        // 计算误差  
        double max_error = 0.0;  
        for(int i = 0; i < n; i++) {  
            double error = fabs(x_new[i] - x[i]);  
            if(error > max_error) {  
                max_error = error;  
            }  
        }  
  
        // 更新解向量  
        for(int i = 0; i < n; i++) {  
            x[i] = x_new[i];  
        }  
  
        // 检查收敛  
        if(max_error < tolerance) {  
            break;  
        }  
    }  
  
    free(x_new);  
    free(diag);  
    return iter;  
}  
  
// 生成严格对角占优的稀疏矩阵（保证雅克比迭代收敛）  
S* generate_diagonally_dominant_sparse_matrix(int n, double density)  
{  
    // 估计非零元素个数  
    int nz = (int)(n * n * density);  
    if(nz < 3 * n) nz = 3 * n; // 至少每行3个元素  
  
    S* A = sCreate(n, n, nz);  
  
    int count = 0;  
    srand(time(NULL));  
  
    // 先设置对角线元素（保证对角占优）  
    for(int i = 0; i < n; i++) {  
        A->ts[count].i = i;  
        A->ts[count].j = i;  
        A->ts[count].value = (n + 1) + (rand() % 100) * 0.1; // 较大的对角线元素  
        count++;  
    }  
  
    // 添加非对角线元素  
    while(count < nz) {  
        int i = rand() % n;  
        int j = rand() % n;  
  
        if(i != j) {  
            // 检查是否已经存在该位置的元素  
            int exists = 0;  
            for(int k = 0; k < count; k++) {  
                if(A->ts[k].i == i && A->ts[k].j == j) {  
                    exists = 1;  
                    break;  
                }  
            }  
  
            if(!exists) {  
                A->ts[count].i = i;  
                A->ts[count].j = j;  
                A->ts[count].value = (rand() % 100 - 50) * 0.01; // 较小的非对角线元素  
                count++;  
            }  
        }  
    }  
  
    return A;  
}  
  
// 生成测试向量  
double* generate_test_vector(int n)  
{  
    double* b = (double*)malloc(n * sizeof(double));  
    srand(time(NULL) + 1);  
  
    for(int i = 0; i < n; i++) {  
        b[i] = (rand() % 100) * 0.1;  
    }  
  
    return b;  
}  
  
// 计算残差 ||Ax - b||_2double compute_residual(S* A, double* x, double* b)  
{  
    int n = A->nRow;  
    double* Ax = (double*)calloc(n, sizeof(double));  
  
    // 计算 Ax    for(int k = 0; k < A->nz; k++) {  
        int i = A->ts[k].i;  
        int j = A->ts[k].j;  
        Ax[i] += A->ts[k].value * x[j];  
    }  
  
    // 计算残差  
    double residual = 0.0;  
    for(int i = 0; i < n; i++) {  
        residual += (Ax[i] - b[i]) * (Ax[i] - b[i]);  
    }  
  
    free(Ax);  
    return sqrt(residual);  
}  
  
// 高斯消去法（稠密矩阵版本）  
void gaussian_elimination_dense(double** A, double* b, double* x, int n)  
{  
    // 创建增广矩阵  
    double** Ab = (double**)malloc(n * sizeof(double*));  
    for(int i = 0; i < n; i++) {  
        Ab[i] = (double*)malloc((n + 1) * sizeof(double));  
        for(int j = 0; j < n; j++) {  
            Ab[i][j] = A[i][j];  
        }  
        Ab[i][n] = b[i];  
    }  
  
    // 前向消元  
    for(int k = 0; k < n; k++) {  
        // 选主元（部分选主元）  
        int max_row = k;  
        for(int i = k + 1; i < n; i++) {  
            if(fabs(Ab[i][k]) > fabs(Ab[max_row][k])) {  
                max_row = i;  
            }  
        }  
  
        // 交换行  
        if(max_row != k) {  
            double* temp = Ab[k];  
            Ab[k] = Ab[max_row];  
            Ab[max_row] = temp;  
        }  
  
        // 消元  
        for(int i = k + 1; i < n; i++) {  
            double factor = Ab[i][k] / Ab[k][k];  
            for(int j = k; j <= n; j++) {  
                Ab[i][j] -= factor * Ab[k][j];  
            }  
        }  
    }  
  
    // 回代  
    for(int i = n - 1; i >= 0; i--) {  
        x[i] = Ab[i][n];  
        for(int j = i + 1; j < n; j++) {  
            x[i] -= Ab[i][j] * x[j];  
        }  
        x[i] /= Ab[i][i];  
    }  
  
    // 释放内存  
    for(int i = 0; i < n; i++) {  
        free(Ab[i]);  
    }  
    free(Ab);  
}  
  
// 将稀疏矩阵转换为稠密矩阵  
double** sparse_to_dense(S* A)  
{  
    int n = A->nRow;  
    double** dense = (double**)malloc(n * sizeof(double*));  
  
    for(int i = 0; i < n; i++) {  
        dense[i] = (double*)calloc(n, sizeof(double));  
    }  
  
    for(int k = 0; k < A->nz; k++) {  
        int i = A->ts[k].i;  
        int j = A->ts[k].j;  
        dense[i][j] = A->ts[k].value;  
    }  
  
    return dense;  
}  
  
// 释放稠密矩阵内存  
void free_dense_matrix(double** A, int n)  
{  
    for(int i = 0; i < n; i++) {  
        free(A[i]);  
    }  
    free(A);  
}  
  
int main()  
{  
    const int n = 1000;           // 矩阵阶数  
    const double density = 0.01;  // 稀疏矩阵密度  
    const int max_iter = 1000;    // 最大迭代次数  
    const double tolerance = 1e-6; // 收敛容差  
  
    printf("=== 1000阶稀疏线性方程组求解对比 ===\n");  
  
    // 生成测试数据  
    printf("生成严格对角占优的稀疏矩阵...\n");  
    S* A = generate_diagonally_dominant_sparse_matrix(n, density);  
    double* b = generate_test_vector(n);  
    double* x_jacobi = (double*)calloc(n, sizeof(double));  
    double* x_gauss = (double*)calloc(n, sizeof(double));  
  
    printf("矩阵规模: %d x %d\n", n, n);  
    printf("非零元素: %d (密度: %.3f%%)\n", A->nz, (double)A->nz/(n*n)*100);  
  
    // 测试雅克比迭代法  
    printf("\n--- 雅克比迭代法 ---\n");  
    clock_t start = clock();  
    int iterations = jacobi_sparse(A, b, x_jacobi, max_iter, tolerance);  
    clock_t end = clock();  
    double time_jacobi = ((double)(end - start)) / CLOCKS_PER_SEC;  
  
    double residual_jacobi = compute_residual(A, x_jacobi, b);  
  
    printf("迭代次数: %d\n", iterations);  
    printf("计算时间: %.4f 秒\n", time_jacobi);  
    printf("残差范数: %.6e\n", residual_jacobi);  
  
    // 测试高斯消去法  
    printf("\n--- 高斯消去法 ---\n");  
    start = clock();  
    double** A_dense = sparse_to_dense(A);  
    gaussian_elimination_dense(A_dense, b, x_gauss, n);  
    end = clock();  
    double time_gauss = ((double)(end - start)) / CLOCKS_PER_SEC;  
  
    double residual_gauss = compute_residual(A, x_gauss, b);  
  
    printf("计算时间: %.4f 秒\n", time_gauss);  
    printf("残差范数: %.6e\n", residual_gauss);  
  
    // 对比分析  
    printf("\n=== 对比分析 ===\n");  
    printf("时间比 (高斯/雅克比): %.2f\n", time_gauss / time_jacobi);  
    printf("残差比 (雅克比/高斯): %.2e\n", residual_jacobi / residual_gauss);  
  
    // 验证解的差异  
    double max_diff = 0.0;  
    for(int i = 0; i < n; i++) {  
        double diff = fabs(x_jacobi[i] - x_gauss[i]);  
        if(diff > max_diff) {  
            max_diff = diff;  
        }  
    }  
    printf("解的最大差异: %.6e\n", max_diff);  
  
    // 内存使用分析  
    size_t sparse_memory = sizeof(S) + A->nz * sizeof(TRIPLET);  
    size_t dense_memory = n * n * sizeof(double);  
    printf("\n内存使用对比:\n");  
    printf("稀疏矩阵存储: %.2f MB\n", sparse_memory / (1024.0 * 1024.0));  
    printf("稠密矩阵存储: %.2f MB\n", dense_memory / (1024.0 * 1024.0));  
    printf("存储比 (稠密/稀疏): %.1f\n", (double)dense_memory / sparse_memory);  
  
    // 释放内存  
    sFree(A);  
    free(b);  
    free(x_jacobi);  
    free(x_gauss);  
    free_dense_matrix(A_dense, n);  
  
    return 0;  
}

精度

• 高斯消去法：由于是直接法，通常能达到机器精度（~1e-15）
• 雅克比迭代法：达到预设容差1e-6，精度较低但满足大多数工程需求

效率

• 对于严格对角占优矩阵，雅克比迭代收敛较快
• 高斯消去法的时间复杂度为O(n³)，对于n=1000需要较长时间
• 雅克比迭代每次迭代时间复杂度为O(nnz)，总体效率更高