Сравнение методов решения систем линейных алгебраических уравнений

[schulzr]

Iterative Methods for Sparse Linear Systems Second Edition, by Y Saad · 2003 

[Tiriet]

после 2003 года несколько продвинулись итерационные методы, не требующие транспонирования матриц (например, Journal of Computational Physics 291 (2015) 20–33) и всякие algebraic multigrid методы. сходимость у них обычно получше, чем у рассмотренных в статье.

[RealLazyCat]

Пожалуйста, прячьте портянки кода под спойлер

[Jijiki]

метод вращений можно поидее на кватернионах сделать, я его бегло глянул у вас, система уравнений видимо описывает матрицу - если так и если там компоненты врщений, то можно и кватернион вставить поидее

[Tiriet]

вообще, для тестирования солверов СЛАУ есть стандартные тестовые матрицы, например, https://math.nist.gov/MatrixMarket/collections/hb.html. Их специально придумали, чтоб оценивать полезность и практическую применимость методов. А на каких матрицах тестировались методы статьи- непонятно, поэтому оценивать их плюсы и минусы по приведенным данным- несколько сомнительно.

[Tiriet]

и еще с утра подумалось:

		for (k = m; k <= size - 1; k++)
			*(alf + m*size + k) /= amm;
		*(bet +m) /= amm;

        amm= 1.0/amm
        for (k = m; k <= size - 1; k++)
			*(alf + m*size + k) *= amm;
		*(bet +m) *= amm;

минимальное изменение в коде, которое повышает скорость работы на жирных матрицах разика в два. например, на Zen5 FMUL имеет латентность 6 тактов, и задержку между инструкциями 2 такта, а FDIV имеет латентность 15 и задержку между инструкциями 5 тактов. то есть, даже если FDIVы будут исполняться последовательно без перерывов на другие инструкции (а в этом коде так и должно быть, инкремент счетчиков будет идти в других регистрах и на ФПУ влиять не будет никак), и латентность влиять ни на что не будет (инструкции будут идти длинной серией на подготовленных заранее данных), то вот Reciprocal throughput 2 vs 5 даст ускорение в 2.5 раза. вроде мелочь, но если вам надо перемножать матрицы в пару гигов весом- оно ой как влияет. Надеяться, что компилятор вкурит, что вот тут вот деление можно соптимизировать- пока еще нельзя.

[Deosis]

	// градиент: 2*(Ax0 - b)
	auto Gradient = [A, b, n](const double* x, double *result)
	{
		double* ax = new double[n * sizeof(double)];
		memset(ax, 0, n * sizeof(double));
		for (int i = 0; i < n; i++)
		{
			for (int j = 0; j < n; j++)
				*(ax + i) += *(A + i * n + j) * *(x + j);
			*(result + i) = 2.0 * (*(ax + i) - *(b + i));
		}
		delete[] ax;
	};

Свободной памяти у автора вообще крайне много: везде выделяется в 8 раз больше чем нужно, но одновременно используется только одно временное значение

[Sibnventor]

Для решения СЛАУ в несколько тысяч уравнений с плохо обусловленной матрицей, а также для прецизионного численного решения дифференциальных уравнений движения или уравнений прогибов и других задач, где требуется высочайшая точность вычислений и работа с переменными очень больших астрономических или очень маленьких внутриатомных значений, я разработал и с успехом применяю собственную действительную переменную BigDouble. Ее особенность в том, что разрядность можно самому задавать любую, например, 20, 25, 30 или 36 десятичных знаков, требования к экспоненте без реальных ограничений.

Для примера, решение СЛАУ с 1000 уравнений дает максимальную невязку среди всех найденных решений. при заданной разрядности 28 не более 1.7 10 -23, а для 2000 уравнений менее 5.0*10-23. Существенный недостаток - большое время расчетов. Например, 1000 уравнений считается порядка 3-х минут, зато простота применения и высокая точность решения любых уравнений с плохо и даже очень плохо обусловленной матрицей коэффициентов гарантируется.