2

我用 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 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;
}
4

4 回答 4

4

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

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

如果您可以接受许可证，我推荐Harwell Subroutine library。虽然接口 C++ 和 Fortran 并不好玩……

于 2013-03-27T18:46:43.103 回答
3

你的意思是稀疏到什么程度？

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

代码和示例在这里。

这是完成大部分工作的类：

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 映射按行和列存储所有值的位置。对于像您提供的那样非常稀疏的矩阵（密度 < .001），我能够在 1/4 秒内解决一个包含 10000 个方程的系统。

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

于 2013-03-27T18:57:10.033 回答
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;
}
于 2016-03-21T11:41:40.640 回答
-1

试试 PETSC。为此，您需要 CRS（压缩行存储）格式。

于 2013-08-04T11:35:42.800 回答