Метод Гаусса – просто, надежно, точно!

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

Давно известно, что грамотное увеличение объема использованной оперативной памяти приводит к ускорению работы программы. Сегодня этот объем настолько велик, что о нем уже не стоит задумываться. Поэтому воспользуемся старым испытанным методом – построим индекс строк и индекс столбцов матрицы. Трогать столбец со свободными членами нет необходимости, поэтому размеры индексов равны числу строк матрицы. К элементам матрицы будем обращаться не напрямую, а через индекс. Затраты процессора возрастут на 2-3 машинные команды, но это несоизмеримо с затратами на перестановку строк или столбцов. Теперь, при необходимости выполнить перестановки, нам достаточно переставить элементы индексов, чтобы программа выбирала элементы матрицы в нужном порядке. Отдельные индексы для строк и столбцов нужны потому, что главный член матрицы может находится и не на диагонали.

Теперь о поиске главного члена матрицы. Необходимо выбирать наибольший по модулю коэффициент из ВСЕХ элементов квадратной части системы. Хорошо хоть на свободные члены уравнений делить не надо. Выполнять этот выбор приходится на каждом шаге прямого хода метода Гаусса. На это огромное число сравнений, если выполнять его для чисел с плавающей запятой, требуется большое время. Для задач реального времени такие затраты процессорных ресурсов недопустимы. Но есть выход!

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

Для высокой точности расчетов воспользуемся расширенным форматом чисел с плавающей точкой сопроцессора TBYTE. У него мантисса 64 бит и характеристика 16 бит, из которых 15 бит – порядок и 1 бит – знак числа. Для загрузки массивов данных такого формата и такой точности не подходят типовые решения, использующие формат DOUBLE с мантиссой 53 бит. Единственный разумный способ – загрузка массивов из текстовых файлов, где число знаков после запятой может быть доведена до 20-ти под размер мантиссы TBYTE. Если матрица небольшого размера, можно записать все коэффициенты одной строки матрицы в одну строку текста, разделив числа пробелом. Но такой способ плохо пригоден для больших матриц, да еще, как правило, со слабым наполнением, когда более половины коэффициентов нули. Поэтому применим векторный способ записи каждого элемента: N строки, N столбца, число. Именно в таком виде хранятся большие таблицы в реляционных базах данных типа ACCESS или SQL. Строки в этой форме исходных данных не зависят от размера матрицы и несут только полезную информацию. Первой строкой данных нужно указать число строк, чтобы заранее можно было зарезервировать память под матрицу, распределить индексы и контролировать ввод на предмет случайного указания запредельных индексов строки или столбца.

После ввода размера матрицы, создаем индексы строк и столбцов. Начальное заполнение индексов, на первый взгляд тривиально – 0, 1, …, N-1. Но мы же охотимся за производительностью, значит надо заполнить индексы сразу смещением, от начала матрицы. Учитывая, что размер переменной TBYTE=10 байт, смещения для индекса столбцов будут: для первого столбца = 0, второго = 10, третьего = 20 и так далее. Для индекса строк числа станут еще более экзотические – 0, (N+1)*10, , (N+1)*20 и так далее. Выполнив индекс в таком стиле, мы сразу сэкономим массу времени при проведении вычислительных шагов метода Гаусса.

Алгоритм выполнения расчета по шагам будет следующим:

1. Ввод размера матрицы;
2. Построение индексов;
3. Ввод данных для матрицы, с размещением переменных по индексам;
4. Нахождение элемента матрицы с максимальным порядком числа (ЭМП);
5. Перестановка элементов индекса так, чтобы ЭМП стал в левом верхнем углу матрицы;
6. Деление всех элементов верхней строки (ЭВС) на ЭМП, за исключением самого ЭМП. Он будет нужен на обратном ходу метода Гаусса.
7. Вычитание ЭВС из соответствующих элементов нижележащих строк.
8. Возврат на п.4 пока не дойдем до нижнего правого угла матрицы.
9. Значение в нижнем правом углу – это значение последнего искомого элемента X(n). А коэффициент левее него – то, на что нужно умножать X(n), перед вычитанием из значения правого столбца предыдущей строки, чтобы получить X(n-1).
10. Возврат на 9, пока не определим все X.

В следующей статье обсудим конкретные участки программы, вызывающие трудности при работе с индексами и работе переменными типа TBYTE.

Спасибо, что дочитали до конца.