如何在c++中优化稀疏矩阵的Gauss-Seidel例程

How to optimize Gauss-Seidel routine in C++ for sparse matrices?

本文关键字：Gauss-Seidel 例程 c++ 优化      更新时间：2023-10-16

我用c++写了一个用Gauss-Seidel方法求解方程组Ax = b的例程。但是，我想将此代码用于稀疏的特定"A"矩阵(大多数元素为零)。这样，求解器的大部分时间都在忙于将一些元素乘以0。

例如，对于以下方程组:

| 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 |

使用高斯-赛德尔方法，我们将得到x1的以下迭代公式:

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

如你所见，解算器通过将0个元素相乘来浪费时间。由于我使用的是大矩阵(例如，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;
}

编写一个稀疏线性系统求解器是很困难的。很努力。

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

如果许可对你来说是可以接受的，我推荐Harwell子例程库。c++和Fortran的接口并不有趣…

你的意思是多稀疏?

这是一个蹩脚的稀疏实现，应该可以很好地解决线性方程系统。这可能是一个幼稚的实现，我对工业强度稀疏矩阵中通常使用的数据结构知之甚少。

代码和示例在这里。

下面是完成大部分工作的类:

template <typename T>
class SparseMatrix
{
private:
    SparseMatrix();
public:
    SparseMatrix(int row, int col);
    T Get(int row, int col) const;
    void Put(int row, int col, T value);
    int GetRowCount() const;
    int GetColCount() const;
    static void GaussSeidelLinearSolve(const SparseMatrix<T>& A, const SparseMatrix<T>& B, SparseMatrix<T>& X);
private:
    int dim_row;
    int dim_col;
    vector<map<int, T> > values_by_row;
    vector<map<int, T> > values_by_col;
};

其他方法定义包含在ideone中。我不测试收敛性，而只是简单地循环任意次数。

稀疏表示使用STL映射，按行和列存储所有值的位置。我能够在1/4秒内解决一个包含10000个方程的系统对于一个非常稀疏的矩阵，就像你提供的那样(密度<措施)。>

我的实现应该足够通用，以支持任何整数或用户定义的类型，支持比较，4个算术运算符(+, -, *, /)，并且可以从0显式强制转换(空节点被给定值(T) 0)。

最近，我也面临同样的问题。我的解决方案是使用向量数组来保存稀疏矩阵。下面是我的代码:

#define PRECISION   0.01
inline bool checkPricision(float x[], float pre[], int n) {
    for (int i = 0; i < n; i++) {
        if (fabs(x[i] - pre[i]) > PRECISION) return false;
    }
    return true;
}
/* mx = b */
void gaussIteration(std::vector< std::pair<int, float> >* m, float x[], float b[], int n) {
    float* pre = new float[n];
    int cnt = 0;
    while (true) {
        cnt++;
        memcpy(pre, x, sizeof(float)* n);
        for (int i = 0; i < n; i++) {
            x[i] = b[i];
            float mii = -1;
            for (int j = 0; j < m[i].size(); j++) {
                if (m[i][j].first != i) {
                    x[i] -= m[i][j].second * x[m[i][j].first];
                }
                else {
                    mii = m[i][j].second;
                }
            }
            if (mii == -1) {
                puts("Error: No Solution");
                return;
            }
            x[i] /= mii;
        }
        if (checkPricision(x, pre, n)) {
            break;
        }
    }
    delete[] pre;
}

试试PETSC。您需要CRS(压缩行存储)格式。

相关文章：