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

Question

我用C++编写了一个例程,用Gauss-Seidel方法求解方程组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乘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子程序库.虽然接口C++和Fortran并不好玩......