c++ - 如何在 C++ 中针对稀疏矩阵优化 Gauss-Seidel 例程？

Question

我用 C++ 编写了一个例程，使用高斯-赛德尔方法求解方程组 Ax = b。但是，我想将此代码用于稀疏的特定“A”矩阵(大多数元素为零)。这样，该求解器花费的大部分时间都忙于将一些元素乘以零。

例如，对于以下方程组:

| 4 -1  0  0  0 | | x1 |   | b1 |
|-1  4 -1  0  0 | | x2 |   | b2 |
| 0 -1  4 -1  0 | | x3 | = | b3 |
| 0  0 -1  4 -1 | | x4 |   | b4 |
| 0  0  0 -1  4 | | x5 |   | b5 |

使用Gauss-Seidel方法，我们将得到x1的以下迭代公式:

x1 = [b1 - (-1 * x2 + 0 * x3 + 0 * x4 + 0 * x5)]/4

如您所见，求解器通过乘以零元素来浪费时间。由于我使用大矩阵(例如 10^5 x 10^5)，这会对总 CPU 时间产生负面影响。我想知道是否有一种方法可以优化求解器，使其忽略与零元素乘法相关的那些计算部分。

请注意，上例中“A”矩阵的形式是任意的，求解器必须能够处理任何“A”矩阵。

代码如下:

void GaussSeidel(double **A, double *b, double *x, int arraySize)
{   
const double tol = 0.001 * arraySize;
double error = tol + 1;

for (int i = 1; i <= arraySize; ++i)
    x[i] = 0;

double *xOld;
xOld = new double [arraySize];
for (int i = 1; i <= arraySize; ++i)
    xOld[i] = 101;

while (abs(error) > tol)
{

    for (int i = 1; i <= arraySize; ++i)
    {
        sum = 0;
        for (int j = 1; j <= arraySize; ++j)
        {
            if (j == i)
                continue;
            sum = sum + A[i][j] * x[j];
        }
        x[i] = 1 / A[i][i] * (b[i] - sum);
    }

    //cout << endl << "Answers:" << endl << endl;
    error = errorCalc(xOld, x, arraySize);

    for (int i = 1; i <= arraySize; ++i)
        xOld[i] = x[i];

cout << "Solution converged!" << endl << endl;
}

Answer 1

编写稀疏线性系统求解器很困难。 非常难。

我只会选择一个现有的实现。任何合理的 LP 求解器内部都有一个稀疏线性系统求解器，参见示例 lp_solve , GLPK等

如果您可以接受许可证，我推荐 Harwell Subroutine library .尽管连接 C++ 和 Fortran 并不有趣...