Comments 20
Iterative Methods for Sparse Linear Systems Second Edition, by Y Saad · 2003
Пожалуйста, прячьте портянки кода под спойлер
метод вращений можно поидее на кватернионах сделать, я его бегло глянул у вас, система уравнений видимо описывает матрицу - если так и если там компоненты врщений, то можно и кватернион вставить поидее
вообще, для тестирования солверов СЛАУ есть стандартные тестовые матрицы, например, https://math.nist.gov/MatrixMarket/collections/hb.html. Их специально придумали, чтоб оценивать полезность и практическую применимость методов. А на каких матрицах тестировались методы статьи- непонятно, поэтому оценивать их плюсы и минусы по приведенным данным- несколько сомнительно.
и еще с утра подумалось:
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 раза. вроде мелочь, но если вам надо перемножать матрицы в пару гигов весом- оно ой как влияет. Надеяться, что компилятор вкурит, что вот тут вот деление можно соптимизировать- пока еще нельзя.
// градиент: 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 раз больше чем нужно, но одновременно используется только одно временное значение
Для решения СЛАУ в несколько тысяч уравнений с плохо обусловленной матрицей, а также для прецизионного численного решения дифференциальных уравнений движения или уравнений прогибов и других задач, где требуется высочайшая точность вычислений и работа с переменными очень больших астрономических или очень маленьких внутриатомных значений, я разработал и с успехом применяю собственную действительную переменную BigDouble. Ее особенность в том, что разрядность можно самому задавать любую, например, 20, 25, 30 или 36 десятичных знаков, требования к экспоненте без реальных ограничений.
Для примера, решение СЛАУ с 1000 уравнений дает максимальную невязку среди всех найденных решений. при заданной разрядности 28 не более 1.7 10 -23, а для 2000 уравнений менее 5.0*10-23. Существенный недостаток - большое время расчетов. Например, 1000 уравнений считается порядка 3-х минут, зато простота применения и высокая точность решения любых уравнений с плохо и даже очень плохо обусловленной матрицей коэффициентов гарантируется.
Жаль, что такую интересную разработку не замечают!
1000 уравнений. ээээх. гидродинамика сварочной ванны при сплавлении порошков в 3д постаовке. полмиллиона уравнений на грубой сетке. Расчет НДС полимерной матрицы с добавкой углеволокна- матрица весит... а сколько есть у тебя оперативы- столько и будет матрица весить. хочется загнать в нее сетку 1000*1000*1000, но она не влезает, так как на такой расчет надо под терабайт оперативы. при этом система вроде бы решается вида Ах=b, но элементами векторов x и b сами по себе являются трехмерные векторы, а элементами матрицы А- тензоры 3*3. так че там у Вас? 1000 уравнений три минуты шарашит.... так это неприемлемо низкая скорость счета даже для двумерных сеток 100*100, не говоря уже о 100*100*100.
еще момент- если всеже нужно упороться, и повысить точность представления чисел в СЛАУ, то частенько удобнее не повышать точность работы с отдельными числами, а саму СЛАУ разбить на два кусочка (A+A')(x+x')=b+b', где А- это старшие разряды в обычном FP64, а A'- это ее уточнение, младшие разряды, в том же самом FP64 (естественно, с уменьшенной на 16разрядов экспонентой). решать четыре СЛАУ с обычными умножениями и делениями сильно быстрее, чем решать одну СЛАУ с нестандартными умножениями и делениями. При этом, даже если Вам нужно именно повышенная точность представления элементов самой матрицы (а не только расчетов с ней), то все равно- на вычисление этой матрицы в виде А+А' времени уходит пренебрежимо мало по сравнению с самим решением СЛАУ. А если все-таки припрет, и реально нужна запредельная точность- то тогда можно вообще на длинную арифметику рациональных чисел переходить- там расчеты абсолютно точные. только долгие жутко.
Когда-то давным-давно в начале 70-х мне приходилось решать уравнения Пуассона на матрице 1000*1000 при оперативной памяти допотопной ЭВМ Одра-1204 всего 11 К, однако все получалось, диплом защитил на отлично! Кто Вам запрещает использовать дисковую память, периодически поочередно читая и записывая после обработки фрагменты матрицы, которые тянет оперативка? Можно очень эффективно ускорять решение, распараллеливая процесс по ядрам, кроме того, итерационные методы ускоряются подбором оптимального коэффициента верхней релаксации, если, конечно процесс сходящийся.
Для решения уникальных по сложности задач разрабатываются специальные методы с учетом имеющихся возможностей. Иногда, для получения максимальной заданной точности результата можно пожертвовать временем. А вот длинная арифметика годится только лишь очень коротеньких сверх точных расчетов или для проверки точности новых алгоритмов, когда есть такая необходимость.
Средненький рассчет- 100килошагоа. На нем одно решение СЛАУ весом в гиг. Итого- 100ТБ данных надо прокачать. Какой там у нас ресурс дисков на запись нынче? Как раз на один рассчет и хватит. А потом все. Новый ссд покупать. Дальше. Гиг надо прокачать через диск туда и обратно, напомню, 100тыс раз. Сколько щас запись на диск? Средненькая- тудасюда- 2гб/сек. Чисто на прокачку через диск туда-сюда уйдет 50тысяч секунд. Это если у меня хорошо все блочно разбивается. В общем, не верю я Вам.
Если задача того стоит, надо ее максимально упростить и решать на суперкомпьютере, а на нашем железе даже начинать не стоит. Как говорил один герой известного фильма: "Бери ношу по себе, чтоб не падать при ходьбе".
зачем мне это, если моя задача на 100к шагов с гиговой матрицей вполне спокойно решается на средней мощности десктопном проце (AMD R9/ Xeon2011v3) за вполне вменяемые десяток часов/несколько суток счета. без свопа, без всяких выкрутасов, обычным (ну, тут кому как, кому обычным, а кому и может сложноватым) методом типа TFQMR2. и хоть мне на каждом шаге приходится решать это самое уравнение Пуассона на трехмерной сетке в полмиллиона узлов и с векторными коэффициентами (у нее ширина ленты получается порядка 10к элементов, и ненулевых в ленте под сотню тензоров, а в двумерном Пуассоне на 1000*1000 ширина ленты будет только 3000 элементов и ненулевых строго 9 штук скаляров)- этот TFQMR2 дает физически годное решение в среднем за 30-50 итераций, там где релаксационные методы хорошо если хоть куда-то примерно туда сходятся итераций за 300. это ж все проверено давно уже, методы Крыловского типа не просто так сочинили- они для физики сильно лучше сходятся, чем релаксационные.
В общем это я все к чему- матрица 1000*1000- это очень скромная в текущих реалиях матрица. Это всего миллион коэффициентов, она вся целиком и полностью влезает в кэш процессора несколько раз. поэтому ворочать ее по три минуты- это неприемлемо долго, и даже для FP256 это вызывает сомнения в выбранном способе реализации расчета именно по части реализации базовых операций умножения и сложения чисел повышенной точности.
Что за матрица с плохой обусловленностью у Вас была? не Гильберта, случаем? дайте пример- мне интересно стало поковыряться, че там за задача.
Последний раз я решал уравнение Пуассона, когда рассчитывал трехмерное гравитационное поле внутри Земли с учетом распределения плотности вещества и формы эллипсоида вращения. Но я немного упростил задачу, свел ее к уравнению Лапласа и получил классное решение при разбиении планеты на кольца сечением 1000 х1000 метров. Этих колец (здесь бы рисунок хорошо вставить) получилось несколько миллионов. В итоге рассчитал распределение гравитационного потенциала и напряженности поля в любой точке объема планеты. На поверхности Земли на уровне океана получилось хорошее соответствие измеряемым величинам ускорения силы тяжести на любой широте, а также известным формулам Гельмерта, Киссиниса. Так что мы тоже чего-то можем!
Я не сталкивался с матрицами Гилберта, но проблема плохой обусловленности у меня возникла при решении задач расчета электрических полей методом интегральных уравнений или методом коллокаций. Там, при значительном смещении фиктивных зарядов внутрь металлических электродов, соседние уравнения имеют весьма близкие коэффициенты и при решении больших СЛАУ с переменными обычной разрядности бывают очень большие ошибки. Приходится жертвовать 2-3 минуты для одноразового точного решения используя переменные с повышенной разрядностью.
Кроме того, я слышал, но не уверен, что проблема плохо обусловленных матриц возникает при решении систем уравнений прогибов в трехмерных прочностных расчетах теории сопротивления материалов.
Основное применение переменных повышенной разрядности у меня в решении дифференциальных уравнений движения тел при взаимном воздействии гравитации. Численное решение любой задачи трех тел с ними вообще не проблема. А вот решение задачи шести тел, а тем более 12 тел без них не получается совсем.
у нас с Вами совершенно разные классы задач. Вам приходится жертвовать 2-3минутами, а у меня нередки рассчеты на пару недель молотилки отдельно стоящего сервака. И в них я не стремаюсь ассемблерными вставками пользоваться (даже тут статью на хабре как-то запостил по этому поводу), чтобы время обращения матриц сократить.
про движение 12 тел: а есть че почитать- накой ляд там повышенная точность представления чисел? у Вас 16 знаков после запятой- это погрешность в 1 атом водорода на 1млн км расстояния. Вы с подобной точностью даже начальную массу своих движущихся тел задать не можете. поэтому если у Вас нормальная консервативная схема (ну или хотя бы полуконсервативная)- то точность решения как таковая вообще ни на что не влияет, ибо даже если Вы будете решать эти уравнения абсолютно точно- погрешность задания начальных условий вам все равно точность решения обломает. Мне небесную механику читали профессиональные астрономы из НИИ из отдела расчета движения тел в космосе. Я с тех пор запомнил, что у них хорошо все считалось Рунге-Куттой 24го порядка на обычных двойных числах. Хорошо все считалось- это десятки небесных тел на интервалах в миллионы лет. а у Вас тут 12 тел и уже требуют FP128, при этом ни слова про метод интегрирования ОДУ- интрига, однако. потому что если Вы загнали FP128- то у Вас погрешность численной схемы должна быть ну хотя-бы сравнима с 2-3 ULP, а построить такую численную схему- это надо ой как постараться.
Спасибо за ответ, хотелось бы все таки чем-то быть полезным, хотя мои работы скорее всего безнадежно устарели. Например, решая уравнение Пуассона на большой трехмерной матрице методом итераций, я задавал грубое начальное приближение в каждой точке исходя из приближенного решения уравнения Лапласа, причем совершенно другим методом. После этого требовалось всего несколько итераций решения уравнения Пуассона с применением оптимального параметра верхней релаксации. В итоге получалось самосогласованное электрическое поле с учетом объемного распределения зарядов в интенсивных потоках заряженных частиц в объеме. Это все, больше мне добавить нечего. В настоящее время такие задачи по-видимому никому не интересны.
Что касается космических задач, то я специально придумал тяжелые тестовые задачи 6-ти или 12-ти тел одинаковой массы (порядка массы звезд) для отладки точности алгоритмов численного интегрирования дифференциальных. уравнений движения при взаимном гравитационном влиянии. Характерно и это давно доказано теоретически, что при одинаковой массе эти решения крайне не устойчивы и очень чувствительны и к начальным значениям и к методу интегрирования и к накопленным округлениям переменных на каждом шаге.
Согласен, самый устойчивый метод - это метод Рунге- Кутты, я применяю распространенные методы 2-го или 4-го порядка, и пытаюсь оптимизировать до более высоких порядков, а про метод 24-го порядка слышу впервые. Очень интересно узнать о нем поподробнее.
После того, как мне удалось разработать переменную BigDouble, я первым делом проверил на тестовой задаче точность решения системы уравнений движения при различных заданных значениях разрядности и при прочих равных условиях. Количество шагов интегрирования по времени - 150 000. Результаты примерно такие:
n=15, dr =1*10^-2;
n=20, dr =1*10^-7;
n=25, dr =1*10^-12;
n=30, dr =1*10^-16;
n=35, dr =1*10^-16;
n - разрядность, dr - порядок погрешности расчета координат.
Вывод простой: разрядность влияет до значений порядка 30, дальнейшее увеличение разрядности смысла не имеет, поэтому длинная арифметика не годится от слова совсем.
не верю я Вам. как можно проверить точность решения? только в одном случае- если есть точное решение и приближенное можно с ним сравнить. Откуда у Вас для численной схемы для 12-ти тел точное решение? ниоткуда. его нет. то, что с ростом разрядности у Вас там че-то перестало меняться- скорее всего вызвано тем, что на высокой разрядности погрешность представления числе стала пренебрежимо мала по сравнению с погрешностью схемы, и поэтому погрешность результата перестала зависеть от погрешности самих вычислений, а стала зависеть только от погрешности схемы. Но это все ваще не так тестируется, а путем пересчета с уменьшенным шагом по времени. Берете не 150к шагов с dt, а 300к шагов с dt/2, а потом- 600к с dt/4, и устанавливаете факт наличия или отсутствия схемной сходимости. и только если схемная сходимость достигается- только после этого можно анализировать- а как на эту сходимость влияет изменение разрядности.
у Вас 12 тел. Вы определяете погрешность координаты. координаты какой? центра масс? так на него влияет не разрядность чисел, а консервативность схемы. Погрешность координаты определяете, ок, а погрешность полной энергии системы при n=15 определяется? нет? тогда опять бессмысленно повышать разрядность, если мы не уверены в том, что интегралы решения сохраняются.
Про Рунге-Кутты 24го порядка: не понял вопроса. что про него читать? рекурсивная формула для вычисления коэффициентов схемы для любого числа стадий, хоть для явного, хоть для неявного случая. 24й порядок точности достигается на неявной 12ти стадийной схеме Рунге-Кутты, любой учебник по численным схемам это описывает. задачи жесткие, поэтому для них нужно именно неявные РК использовать, они поустойчивее.
Мне совершенно нет никакого интереса вас обманывать. я ничего не продаю и не собираюсь продавать, просто делюсь информацией по результатам собственных исследований. Задача 12 тел, так же как и задача 6-ти тел в одном частном случае имеет точное аналитическое решение. Это решение возможно, когда тела расположены полностью симметрично относительно центра масс, имеют одинаковые массы и одинаковые по модулю скорости, при которых центробежная сила, действующая на каждое тело равна сумме гравитационных сил от других тел (формула не сложная). Тогда, если задать равновесную скорость все тела будут двигаться по кругу с постоянным радиусом. Если задать скорость меньше равновесной, то тела будут двигаться по эллипсам вокруг центра масс.
Не надо применять метод деления шага интегрирования пополам, достаточно сравнить аналитический радиус орбиты вращения и радиус, полученный как корень квадратный из суммы квадратов приращений координат. Можно сравнить результаты даже на первом шаге интегрирования.
Кстати в результате подобных расчетов получился весьма интересный вывод, что шаг интегрирования уравнений движения по времени имеет оптимальное значение. То-есть при некотором шаге интегрирования дальнейшее его уменьшение приводит к повышению погрешности, за счет увеличения ошибок округления при малых шагах. При этом рекомендуемое применение метода деления шага интегрирования часто приводит к прямо противоположному результату.
И еще, опыт расчетов движения систем тел во взаимных гравитационных полях однозначно показал, что при высокой точности интегрирования, все законы сохранения и теорема движения центра масс выполняются автоматически. К примеру, даже при самой низкой точности задания начальных эфемерид всех планет и Солнца (16 разрядов в таблицах NASA), расчет их траекторий приводит к смещению вычисленного значения центра масс Солнечной системы всего на 14.5 метра на интервале 100 лет, хотя метод интегрирования был далеко не самый точный (Рунге-Кутты 2-го порядка).
у нас с Вами какое-то сильно разное образование, походу.
Кстати в результате подобных расчетов получился весьма интересный вывод, что шаг интегрирования уравнений движения по времени имеет оптимальное значение. То-есть при некотором шаге интегрирования дальнейшее его уменьшение приводит к повышению погрешности, за счет увеличения ошибок округления при малых шагах.
это базовый момент теории численных схем- два источника погрешностей- погрешность аппроксимации уравнений разностной схемой и погрешность аппроксимации чисел. и делать первую меньше второй- нет никакого смысла. Кроме того, оптимальный шаг по времени- это не какое-то "интересный вывод", а рутинное и бытовое действие- выбор оптимального шага по времени позволяет поднять порядок точности схемы как правило- на единичку, а иногда и на двоечку- если в разложении погрешности в ряд отсутствуют четные слагаемые. и это- повторю- рутинная процедура в теории разностных схем- расчет с оптимальным шагом по времени имеет более высокий порядок точности, а потому его можно делать на более грубой сетке и значит- быстрее.
Сохранение энергии (и вообще интегралов уравнения- они же разные бывают, не только энергия, но импульс, момент импульса, координаты центра тяжести, всякие интегралы Римана и прочее)- это следствие того, что в схеме есть аппроксимация, сходимость и устойчивость. Но есть схемы, которые принципиально сохраняют эти интегралы, хоть че ты делай и хоть с каким шагом по времени считай, а есть схемы- которые сохраняют эти интегралы примерно. с хорошей точностью, но примерно. Первые- консервативные- как правило сильно стабильнее вторых.
Центр масс Солнечной системы на интервале 100 лет уплывал на 10-8. ну, че я могу сказать- хреновый результат, однако. 100 лет- это 8 оборотов Юпитера, устойчивость планетных систем рассчитывается на интервалах в десятки и сотни миллионов лет. то есть, на интервалах в 10^6 степени раз больше. за такое время уплывание центра тяжести системы будет уже порядка 1%. а с такой точностью там ваще ловить нечего.
поэтому я и говорю, что излагаемые Вами мысли вызывают у меня удивление и некоторый скепсис. И поэтому мне интересно- какая была постановка задачи и начальные условия и что именно тестировалось- потому что пока из Ваших слов у меня складывается впечатление, что эти все результаты элементарно достигаются на обычных fp64 с какой-нибудь явной рунге-куттой 8го порядка даже без заморочек с неявными повышенной точности и с заморочками с консервативностью.
Сравнение методов решения систем линейных алгебраических уравнений