2

ガウス・ザイデル法を使用して連立方程式 Ax = b を解くルーチンを C++ で作成しました。ただし、スパースな特定の "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

スパース線形システム ソルバーを作成するのは困難です。とても厳しい。

既存の実装の 1 つを選択するだけです。合理的な LP ソルバーには、内部にスパース線形システム ソルバーがあります。たとえば、 lp_solve、GLPKなどを参照してください。

ライセンスが受け入れられる場合は、Harwell Subroutine ライブラリをお勧めします。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;
};

他のメソッド定義はイデオンに含まれています。収束をテストするのではなく、単純に任意の回数ループします。

スパース表現は、STL マップを使用して、すべての値の位置を行と列ごとに格納します。あなたが提供したような非常にまばらな行列 (密度 < .001) について、10000 の方程式のシステムをわずか 1/4 秒で解くことができます。

+私の実装は、比較、4 つの算術演算子 ( 、-、*、/) をサポートし、0 から明示的にキャストできる (空のノードには値 が与えられる)整数型またはユーザー定義型をサポートするのに十分な汎用性を持つ必要があります(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 に答える