Решение СЛАУ методом ГАУССА - C#
Формулировка задачи:
По роду деятельности зачастую необходимо решать СЛАУ вида A*x=b A[N,N] x[N] b[N]
До недавних пор использовал библиотеку ALGLIB для .NET Но(!!!) буквально на днях понадобилось решать несколько тысяч систем при N=20k. Половина ресурсов сервака проработала весть день и ночь исключительно на мою задачу )))) Стал разбираться, почему же с ростом размерность так деградирует скорость - оказалось все в высшей степени тривиально. NET не умеет оптимизировать 2х мерные массивы - т.е. при запросе A[i,j] вызывается функция GetElement(i,j) которая внутри себя берет адрес начала массива, добавляет к нему i*(число столбцов)+j, извлекает адрес и возвращает значение. (!) Кстати, при работе с одномерными массивами скорость также деградирует при запросе больших(длинных) индесков, хотя и не так сильно, как для многомерных.Решение 1.
Сказал нафиг, скачал BLAS(LAPACK) и стал юзать ее посредством PInvoke. Довольно долго разбирался, почему решает неправильно - выяснил, что выравнивание матриц происходит в памяти по столбцам, т.е. A(0,0),A(1,0),A(2,0),... В С# же построчно - A(0,0),A(0,1),A(0,2),... В итоге транспонирование матрицы в C# перед передачей в PInvoke решило задачу. И о чудо - скорость возросла в 20-50 раз.Решение 2.
Подумал, а почему бы не написать самому, по умному алгоритм на C#, который будет быстро перемножать матрицы. Весь день убил и в итоге написал алгоритм, который по скорости оказался быстрее(!) LAPACK и не требовал транспонирования матрицы. Вот финальный результат:// result - указатель на буффер правой части СЛАУ // mtx - указатель на первый элемент матрицы // n - размер системы static bool _Solve_square_system(double* result, double* mtx, int n) { double* mtx_u_ii, mtx_ii_j; double a; double* mtx_end = mtx + n * n, result_i, mtx_u_ii_j = null; double mx; int d = 0; // mtx_ii - указатель на 'ведущий' диагональный элемент матрицы for (double* mtx_ii = mtx, mtx_ii_end = mtx + n; mtx_ii < mtx_end; result++, mtx_ii += n + 1, mtx_ii_end += n, d++) { // Смотрим под ведущим элеменом, значение максимальное по модулю { mx = System.Math.Abs(*(mtx_ii_j = mtx_ii)); for (mtx_u_ii = mtx_ii + n, result_i = result + 1; mtx_u_ii < mtx_end; mtx_u_ii += n, result_i++) { if (mx < System.Math.Abs(*mtx_u_ii)) { mx = System.Math.Abs(*(mtx_ii_j = mtx_u_ii)); mtx_u_ii_j = result_i; } } // если максимальный по модулю элемент равен 0 - система вырождена и не имеет решения if (mx == 0) return false; // если максимальный элемент не является ведущим, делаем перестановку строк, чтобы он стал ведущим else if (mtx_ii_j != mtx_ii) { a = *result; *result = *mtx_u_ii_j; *mtx_u_ii_j = a; for (mtx_u_ii = mtx_ii; mtx_u_ii < mtx_ii_end; mtx_ii_j++, mtx_u_ii++) { a = *mtx_u_ii; *mtx_u_ii = *mtx_ii_j; *mtx_ii_j = a; } } } //'обнуляем' элементы над ведущим for (mtx_u_ii = mtx_ii - n, result_i = result - 1; mtx_u_ii > mtx; mtx_u_ii -= n) { a = *(mtx_u_ii) / *mtx_ii; for (mtx_ii_j = mtx_ii + 1, mtx_u_ii_j = mtx_u_ii + 1; mtx_ii_j < mtx_ii_end; *(mtx_u_ii_j++) -= *(mtx_ii_j++) * a) ; *(result_i--) -= *(result) * a; } //'обнуляем' элементы под ведущим for (mtx_u_ii = mtx_ii + n, result_i = result + 1; mtx_u_ii < mtx_end; mtx_u_ii += d) { a = *(mtx_u_ii++) / *mtx_ii; for (mtx_ii_j = mtx_ii + 1; mtx_ii_j < mtx_ii_end; *(mtx_u_ii++) -= *(mtx_ii_j++) * a) ; *(result_i++) -= *(result) * a; } } //матрица приведена к дигональному виду //приводим ее к единичному, получая права решение for (; mtx_end > mtx; *(--result) /= *(--mtx_end), mtx_end -= n) ; return true; } public static bool SolveSquare(double[] outResult, double[,] A, double[] x) { int n = outResult.Length; if (n == x.Length && n == A.GetLength(0) && n == A.GetLength(1)) { if (n > 0) { Array.Copy(x, outResult, n); fixed(double* presult = &outResult[0]) fixed (double* pmtx = &A[0,0]) return _Solve_square_system(presult, pmtx, n); } else return true; } throw new IndexOutOfRangeException(); }
Решение задачи: «Решение СЛАУ методом ГАУССА»
textual
Листинг программы
static public void _SMul(double* rmDestination, double* rmX, int xn, int xm, double* rmY, int ym) { for (double* _dest_end = rmDestination + xn * ym; rmDestination < _dest_end; rmDestination += ym) for (double* _x_i_end = rmX + xm, _dest_i_end = rmDestination + ym, _y_tmp = rmY; rmX < _x_i_end; rmX++) { //сохраняем текущее значение для многократного использования double cache_val = *rmX; for (double* _dest_ij = rmDestination; _dest_ij < _dest_i_end; _dest_ij++, _y_tmp++) *_dest_ij += cache_val * *_y_tmp; } } static public double[,] SMul(double[,] X, double[,] Y) { int xm = X.GetLength(1); if (xm == Y.GetLength(1)) { int xn = X.GetLength(0); int ym = Y.GetLength(1); double[,] result = new double[xn, ym]; if (result.Length > 0) fixed (double* p_result = &result[0, 0]) fixed (double* p_x = &X[0, 0]) fixed (double* p_y = &Y[0, 0]) _SMul(p_result, p_x, xn, xm, p_y, ym); return result; } else throw new ArgumentException(); }
ИИ поможет Вам:
- решить любую задачу по программированию
- объяснить код
- расставить комментарии в коде
- и т.д