Как оптимизировать процедуру Гаусса-Зейделя в C ++ для разреженных матриц?

Question

Как оптимизировать процедуру Гаусса-Зейделя в C ++ для разреженных матриц?

Я написал подпрограмму на 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 |

Используя метод Гаусса-Зейделя, мы получим следующую итерационную формулу для x1:

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

Как видите, решатель тратит время, умножая ноль элементов. Поскольку я работаю с большими матрицами (например, 10 ^ 5 на 10 ^ 5), это отрицательно повлияет на общее время процессора. Интересно, есть ли способ оптимизировать решатель так, чтобы он пропускал те части вычислений, которые связаны с умножением нулевого элемента.

Обратите внимание, что форма матрицы «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;
}

2

c++numerical-methods sparse-matrix

Решение

Другие решения

Писать разреженный линейный системный решатель сложно. ОЧЕНЬ СЛОЖНО.

Я бы просто выбрал одну из существующих реализаций. Любой разумный решатель LP имеет разреженный линейный решатель системы, см., Например, lp_solve, GLPK, и т.п.

Если лицензия приемлема для вас, я рекомендую Библиотека подпрограмм Harwell. Взаимодействие C ++ и Фортран не весело, хотя …

5

В последнее время я сталкиваюсь с той же проблемой.
Мое решение — использовать векторный массив для сохранения разреженной матрицы.
Вот мой код:

#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;
}

0

Попробуйте PETSC. Для этого вам нужен формат CRS (сжатое хранилище строк).

-1

Источник

Accepted Answer

Как редко вы имеете в виду?

Вот дрянная редкая реализация, которая должна хорошо работать для решения систем линейных уравнений. Вероятно, это наивная реализация, я очень мало знаю о структурах данных, обычно используемых в разреженных матрицах промышленного уровня.

Код и пример здесь.

Вот класс, который выполняет большую часть работы:

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. Я могу решить систему из 10000 уравнений всего за 1/4 секунды для очень разреженной матрицы, подобной той, которую вы предоставили (плотность < +0,001).

Моя реализация должна быть достаточно универсальной, чтобы поддерживать любой интегральный или определенный пользователем тип, который поддерживает сравнение, 4 арифметических оператора (+, -, *, /), и это может быть явным образом приведено из 0 (пустым узлам присваивается значение (T) 0).

2