下面是一个广泛的问题。

在 NxM 的网格上给出了带有 BC 的泊松方程。

重写为 N*M 线性方程组，

转换为稀疏格式（CSR/CRS/Yale），使用Conj Grad 方法求解

至于 100x100 网格计算时间约为 0.3 秒，如果我将网格增加到 500x500，代码已经运行了 11 秒。

最初是在寻找可以用作 Navier-Stokes 求解器子程序的快速方法，即会多次运行 Poisson 求解器。

问题 - 最佳计算时间是多少以及如何正确选择网格大小。或者，也许，我应该完全改变方法，还是它有问题？

下面发布了我的 CSR 输入格式的 conj 函数：

double *cj(double *val, int *col, int *row, int valSize, int colSize, int rowSize, double *b, int ansSize){
    double *x; 
    x = (double *) malloc(ansSize * sizeof(double));
    double *r;
    r = (double *) malloc(ansSize * sizeof(double));
    double *p;
    p = (double *) malloc(ansSize * sizeof(double));
    double *Ap;
    Ap = (double *) malloc(ansSize * sizeof(double));
    double rsold, alpha, rsnew, sum;
    int colStart = 0;
    int i, j, l, k, rowStart, rowEnd, colEnd;
    for (i = 0; i < ansSize; i++){// fill x with 0
        x[i] = 0;
    }
    for (i = 0; i < ansSize; i++){
        rowStart = row[i];
        rowEnd = row[i+1];
        sum = 0;
        for (j = rowStart; j < rowEnd; j++){
            k = col[j];
            sum = sum + val[j]*x[k];
        }
        r[i] = b[i] - sum;
    }
    for (i = 0; i < ansSize; i++){
        p[i] = r[i];
    }
    rsold = 0;
    for (i = 0; i < ansSize; i++){
        rsold = rsold + r[i] * r[i];
    }
    for (i = 0; i < ansSize; i++){
        for (l = 0; l < ansSize; l++){
            rowStart = row[l];
            rowEnd = row[l+1];
            sum = 0;
            for (j = rowStart; j < rowEnd; j++){
                k = col[j];
                sum = sum + val[j]*p[k];
            }
            Ap[l] = sum;
        }
        sum = 0;
        for (j = 0; j < ansSize; j++){
            sum = sum + p[j] * Ap[j];
        }
        alpha = rsold / sum;
        for (j = 0; j < ansSize; j++){
            x[j] = x[j] + alpha * p[j];
        }
        for (j = 0; j < ansSize; j++){
            r[j] = r[j] - alpha * Ap[j];
        }
        rsnew = 0;
        for (j = 0; j < ansSize; j++){
            rsnew = rsnew + r[j] * r[j];
        }
        if (sqrt(rsnew) < epsilon){
            return x;
        }
        for (j = 0; j < ansSize; j++){
            p[j] = r[j] + (rsnew/rsold)  * p[j];
        }
        rsold = rsnew;
    }
    free(r);
    free(p);
    free(Ap);
    return NULL;
}