Параллельная обработка данных 14.12.2000
В.А. Фисун /ВМиК МГУ/
Рекомендую просматривать через Netscape, поскольку Explorer не умеет показывать многоэтажные индексы
СОДЕРЖАНИЕ
Параллельные алгоритмы

1. Введение
2. Измерения вычислительных алгоритмов
2.1. Мера параллелизма
3. Параллельный алгоритм суммирования
4. Параллельные формы рекурсивных алгоритмов
5. Каскадное суммирование
6. Согласованность параллельных алгоритмов
7. Циклическая редукция для решения уравнения рекурсии 1 порядка
8. Параллельные алгоритмы и производительность ЭВМ

Вычисления на ЭВМ

1. Источники вычислительных погрешностей
2. Погрешности исходных данных
3. Формы представления чисел в ЭВМ
4. Числа с плавающей запятой
4.1. Представление систем с плавающей точкой в ЭВМ
5. Особенности арифметики с плавающей запятой
5.1. Нарушение алгебраических законов в плавающей арифметике
6. Точность плавающей арифметики
7. Погрешности округления при вычислениях
7.1. Погрешности выполнения арифметических операций
8. Оптимизационные преобразования программ
1. Балансировка дерева вычислений
2. Исключение общих подвыражений
3. Разворачивание циклов
9. Блокировка оптимизационных преобразований
10. Особенности вычисления выражений в Фортране 90
11. Погрешности асинхронного распараллеливания

Параллельные алгоритмы

1. Введение

Одно из формальных определений алгоритма следующее. Алгоритм - процесс последовательной обработки структур данных, протекающий в дискретном времени так, что в каждый следующий момент времени структура данных получается по заданным правилам по структуре данных величин, имевшихся в предыдущий момент.
Структура элементарных шагов в определении алгоритма не оговаривается, поэтому, правила преобразования величин могут допускать или предписывать параллельную их обработку. Алгоритм может быть параллельным, когда некоторые шаги обработки могут выполняются параллельно. Так, алгоритм сложения векторов может быть представлен как последовательный, в виде: Ci = Ai + Bi для i=1,2..n, или в матричной записи: C = A + B, семантика которого допускает параллельную обработку элементов векторов.

2. Измерения вычислительных алгоритмов

Основной единицей измерения вычислительных алгоритмов является число операций над числами с плавающей запятой "флопов" (flop). Для реализации алгоритмов на параллельных вычислительных системах особый интерес представляет измерения параллелизма алгоритмов.

2.1. Мера параллелизма

"Реальному миру присущ параллельный характер, в то же время наше математическое восприятие в течении 300 лет находилось под воздействием, главным образом, последовательной математики, 50 лет - под воздействием создания последовательных алгоритмов, 15 лет - под воздействием программирования на последовательном Фортране."
/К. Тербер 1976 г./
При реализации на современных вычислительных системах вычислительных алгоритмов доля их параллельных вычислений непосредственно влияет на время их работы.
"Степень параллелизма" /параллельная сложность/ численного алгоритма называется число его операций, которые можно выполнить параллельно. Степень параллелизма алгоритма сложения векторов равна n.
Операции, которые всегда могут быть выполнены на идеальной ЭВМ параллельно независимо от их количества, иногда называют операциями фиксированной глубины.
Последовательный алгоритм сложения чисел: S = A1, S = S + Ai, i=2,..n имеет нулевую степень параллелизма.

3. Параллельный алгоритм суммирования

Параллельно суммирование последовательности n чисел можно произвести так: на первом этапе складываются соседние числа. Полученные суммы также складываются попарно, и т.д. Для n = 2q алгоритм состоит из q = log n этапов, на первом этапе выполняются n/2 сложений (степень параллельности этапа n/2), на втором - n/4 и т. д. Такой алгоритм называется сложение методом сдваивания, он имеет различную степень параллелизма на разных этапах вычислений. Граф, описывающий последовательность операций сложения, граф сдваивания (по Д. Ортега "fan-in graph.") представляет собой двоичное дерево, соответственно, выполняемые операции можно называть операциями на дереве.
Способ реализации процедуры суммирования данным методом зависит от архитектуры вычислительной системы. При наличии n/2 процессоров эту работу можно выполнить так: на первом этапе одновременно получить суммы четных/нечетных соседних элементов последовательности Ai, т.е. (A1+A2), (A3+A4), ...(An-1+An); затем такая процедура повторяется для суммирования полученных частных сумм и так далее. Если n = 2q, то через q = log2n шагов получается искомая сумма. Однако, потери на синхронизацию вычислений, на пересылки частных сумм могут оказаться сравнимы с временем вычисления суммы двух чисел в каждом процессоре. Поэтому, с учетом особенностей характеристик вычислительной системы, дерево вычислений может быть преобразовано, например, с целью увеличения числа операций, выполняемых в узлах, повышения "зернистости" алгоритма.
Алгоритм сдваивания реализуются также в блоках оптимизации компиляторов последовательных ЭВМ для полной загрузки конвейерных вычислителей. Так алгоритм оптимизация "балансировка дерева вычислений" (tree-height reduction or balancing) будет трактовать вычисление суммы вещественных чисел: A+B+C+D+E+F+G+H, как последовательность операций: (((A+B)+(C+D))+((E+F)+(G+H))).
Замечание 1
Алгоритм сдваивания нарушает заданную по-умолчанию последовательность вычислений с накоплением одной частной суммы, что может повлиять на результат.
Замечание 2
Порядком суммирования можно управлять, используя принцип сохранения скобок: компилятор не должен менять последовательность вычислений, предписанных пользователем. Например: запись выражение в виде: (((((((A+B)+C)+D)+E)+F)+G)+H) запрещает суммирование сдваиванием.
Замечание 3
Алгоритм сдваивания является примером широкого класса параллельных алгоритмов, называемых редукционными алгоритмами, иногда, операциями логарифмической глубины.
/ Редукция - упрощение, в биологии, уменьшение размера органа вплоть до его полного исчезновения./
В некоторых работах алгоритм суммирования методом сдваивания называются также алгоритмом каскадных сумм.

4. Параллельные формы рекурсивных алгоритмов

Рекурсия - последовательность вычислений, при котором, значение самого последнего терма в последовательности зависит от одного или несколько ранее вычисленных термов. Пусть группа вычислений может производиться параллельно, использую результаты вычислений, выполненных на предыдущих этапах (полученных в виде начальных данных). Тогда, каждая группа вычислений называется "ярусом" параллельной формы, число групп - "высотой", максимальное число операций в группе "шириной" параллельной формы. Один и тот же алгоритм может иметь несколько представлений в виде параллельных форм, различающиеся как шириной, так и высотой. Редукционный алгоритм сдваивания для суммирования чисел с получением частных сумм может иметь вид:
Данные А1 А2 А3 А4 А5 А6 А7 А8
Ярус 1 А12 А34 А56 А78
Ярус 2 А123 А1234 А567 А5678
Ярус 3 А12345 А123456 А1234567 А12345678
Высота параллельной формы равна трем, ширина - четырем, причем загрузка вычислителей (четырех) полная.
В данном алгоритме производится вычисления пяти 'лишних' чисел по сравнению с последовательным алгоритмом.

5. Каскадное суммирование

Примером параллельных алгоритмов, ориентированных на векторные вычислители, может служить метод вычисления каскадных сумм (алгоритм рекурсивного удвоения) для распараллеливания операций суммирования. Пусть необходимо просуммировать n чисел с сохранением промежуточных сумм: Si = Si-1 + Ai i = 2,..n, S1 = A1. Исходный вектор А поэлементно складывается с вектором Аs, полученный из исходного со сдвигом на один элемент и заполнением позиции элемента А0 нулем. Для вектора результата процедура повторяется, но сдвиг - на 2 позиции. Если n = 2k, то через k операций получается вектор результата. Для i=8:
A1 + 0 = A1 + 0 = A1 + 0 = A1
A2 A1 A12 0 A12 0 A12
A3 A2 A23 A1 A123 0 A123
A4 A3 A34 A12 A1234 0 A1234
A5 A4 A45 A23 A2345 A1 A12345
A6 A5 A56 A34 A3456 A12 A123456
A7 A6 A67 A45 A4567 A123 A1234567
A8 A7 A78 A56 A5678 A1234 A12345678

Алгебра данного метода может быть записана в виде вычисления (возможно, параллельного) частных сумм вида: Si = Ali, где Ali = A(l-1)i + A(l-1)(i-2(l-1)), A0i = Ai для i = 1,2,...n.
Вычисления проводятся l = 0,1,...,log2n раз, причем, если у Ali индекс i выходит из интервала 1<= i <= n то он принимается равным нулю.
Хокни предлагает элегантную векторную форму записи алгоритма каскадного суммирования массива D(n):
     X = D
     DO L = 1,LOG2(N)
        X = X + SHIFTR(X,2**(L-1))
     ENDDO
Результат векторной функции SHIFTR(A,l) есть массив (вектор), полученный из А , элементы которого сдвинуты на L позиций вправо, а L левых элементов заполнены нулями. Практическая реализация алгоритма может исключить излишние операции сложения с нулем, однако, и после этого, по сравнению с последовательным алгоритмом, данный - требует лишние операции.
Замечание:
Такие параллельные формы можно строить также для любых ассоциативных операций, например, умножения, поиск минимакса, умножения матриц.

6. "Согласованность" параллельных алгоритмов

При конструировании параллельных алгоритмов следует иметь ввиду, что они не эквивалентны с точки зрения машинных вычислений, устойчивость параллельных алгоритмов хуже устойчивости последовательных.
Для некоторых скалярных алгоритмов аналогичные параллельные могут быть выгодными только для задачи малой размерности, с ростом порядка задачи время работы параллельной версии может превзойти скалярный расчет. Пусть число арифметических операций в параллельном алгоритме Р(n), а в эквивалентном ему наилучшим последовательном - S(n). Тогда, параллельный алгоритм, решающий задачу порядка n, "согласован" с последовательным алгоритмом для той же задачи, если отношений P(n)/S(n) остается ограниченным при n -> @. Параллельный алгоритм не согласован, если P(n)/S(n) -> @ при n -> @. (@-бесконечность).
Алгоритм рекурсивного удвоения (метод каскадных сумм) не согласован, так как число операций в параллельной версии равно n * log n - n + 1 против n-1 последовательных операций, и (n*log n-n)/n -> @ при n -> @. Однако, для некоторых вычислителей данный алгоритм не бесполезен.
Пусть сложение векторов длины n производится за время T = s+yn, а скалярная операция за m.
Тогда, каскадное вычисление n-1 сумм производится за время s * log n + (y(n*log n -n+1), а последовательное - m(n-1).
И если s=1000, y = 10, m = 200 наносекунд, скалярный алгоритм эффективней векторного при n <= 32 и при n > 221. Некоторые времена:
n t=P(n) t=S(n)
8 3170 1400
32 6190 6200
128 15*103 25*103
1024 0.1*106 0.2*106
221 419*106 415*106

7.Циклическая редукция для решения уравнения рекурсии 1 порядка

Редукция - упрощение, в биологии уменьшение размера органа вплоть до его полного исчезновения. Циклическая редукция - алгоритмы численного анализа для распараллеливания последовательных алгоритмов, основанный на последовательном, циклическом применении параллельных вычислений, число которых на каждом этапе уменьшается (делится) пополам.
Линейной рекурсией 1 порядка называется система уравнений вида:
         X1 = D1
         X2 = X1 * A2 + D2
         .................
         Xi = Xi-1 * Ai + Di
         .................
         Xn = Xn-1 * An + Dn
   
в общем виде: Xi = Xi-1 * Ai + Di, i = 2,3,...n, X1 = D1
Это система эквивалентна двухдиагональной системе уравнений Ax=d, где
A=
1 0 0 0
-a2 1 0 0
........
0 0 -an 1

d=
d1
..
..
dn
Последовательный алгоритм вычислений может быть записан так:
      X(1) = A(1) + D(1)
      DO i = 1,n
        X(i) = X(i-1) * A(i) + D(i)
      ENDDO
 
Рекурсивная зависимость итераций цикла не позволяет ускорить вычисления за счет параллельной работы оборудования. Преобразуем данный алгоритм в параллельный методом циклической редукции. Рассмотрим два соседних уравнения:
         Xi-1 = Xi-2 * Ai-1 + Di-1
         Xi = Xi-1 * Ai + Di
и подставив первое во второе, получаем:
Xi = (Xi-2 * Ai-1 + Di-1) * Ai + Di = Xi-2 * A1,i + D1,i ,
где
A1,i = Ai * Ai-1 , D1,i = Ai * Di-1 + Di
Тогда, проведя эту операцию для всей системы уравнений, получим систему уравнений порядка n/2. Если повторить процедуру l раз (если n = 2l), то в результате получается значение: Xn = Dnl. Для получения полного вектора X необходимо модифицировать алгоритм, например, по аналогии с алгоритмами суммирования.
Очевидно, что вычисления Aji и Dji можно проводить параллельно методом каскадных сумм с сохранением частных сумм. Приведенные уравнения для уровня i имеют вид:
     Xi = Ali * Xi-2l + Dli ,
где l = 0,1,..,log2n , i = 1,2,..,n
Al,i = Al-1,i * Al-1,(i-2l-1) Dl,i = Al-1,i * Dl-1,(i-2l-1) + Dl-1,i Начальные данные: A0,i = Ai, D0,i = Di
Если индекс i у любого Al,i, Dl,i и Xi попадает вне диапазона 1 <= i <= n , то он должен быть приравнен к нулю. Тогда , при l = log2n в уравнениях: Xi = Ali * Xi-2l + Dli индекс Xi-2l = Xi-n находится вне диапазона, и, следовательно, решением системы уравнений будет: вектор: Xi = Dli,
Нотация Хокни для данного алгоритма:
     X = D
     DO L = 1,LOG2(N)
        X = A * SHIFTR(X,2**(L-1)) + X
        A = A * SHIFTR(A,2**(L-1))
     ENDDO

8. Параллельные алгоритмы и производительность ЭВМ

В общем случае, ускорение параллельного алгоритма как мера параллелизма может быть определено как отношение Sp = <время выполнения алгоритма на одном процессоре> / <время выполнение алгоритма на N процессорах>. Очевидно, что ускорение алгоритма сложения двух векторов размерности: Ci = Ai + Bi, i = 1,..,N с "идеальным" параллелизмом равно N. Однако, реальные задачи, кроме вычислений, которые могут выполняться параллельно, содержат и чисто последовательный вычисления. С другой стороны, для практики представляет интерес не параллелизм, присущий алгоритму, а ускорение, реально получаемое при реализации алгоритма на конкретном оборудовании. Для грубой оценки ускорения с учетом степени параллелизма программы и числа мультипроцессоров используется формула, известная как закон Амдала. Закон Амдала показывает коэффициент ускорения выполнения программы на параллельных системах в зависимости от степени распараллеливания программы.
Пусть: N - число процессоров системы, P - доля распараллеливаемой части программы, а S = 1-P - доля операций программы, выполняемых без совмещения по времени.( S+P = 1 , S,P >= 0)
Тогда, по Амдалу, общее время выполнения программы на однопроцессорном вычислителе S+P, на мультипроцессоре: S+P/N, а ускорение при этом есть функция P: Sp = (S+P)/(S+P/N) = 1/(S+P/N)
Из формулы закона Амдала следует,что при:
P = 0 S = 1 - ускорения вычисления нет, а
P = 1 S = 0 - ускорение вычислений в N раз
Если P = S = 0.5, то даже при бесконечном числе процессоров ускорение не может быть более чем в 2 раза.
Однако, практика показала, что независимость N и P условна, ибо в ряде случаев можно для повышения P перепрограммировать задачу, например, увеличить ее размерность. Густавсон предложил версию масштабируемого или модифицированного закона Амдала: Sp = S + P*N
Тогда, при P = S = 0.5 Sp = 0.5N + 0.5
Замечание: N - представленное как число процессоров мультисистемы, может быть:числом АЛУ в микропроцессоре, числом конвейерных степеней в АЛУ.
Литература

Хокни Р., Джессхоуп К. Параллельные ЗВМ. М.: Радио и связь, 1986.-392 с.
Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. М. Мир 1991.-367 с.

Вычисления на ЭВМ

1. Источники вычислительных погрешностей

Первоначально вычислительные машины разрабатывались и использовались для автоматизации исключительно вычислительных работ (COMPUTER - вычислитель), т.е. для выполнения над действительными (целыми) числами четырех арифметических действий: +,-,*,/. Сейчас под вычислениями понимаются и обработка текстовой, графической информации, исчисления множеств, работа с графами, а доля чисто арифметических вычислений довольно небольшая. Тем не менее предметная область с арифметическими вычислениями была и остается одной из важнейших сфер применения ЭВМ.
При этом имеет место парадокс /Г.К.Боровин/: современные ЭВМ из всего спектра задач менее всего приспособлены для выполнения арифметических вычислений из за способа представления чисел в памяти ЭВМ. ЭВМ допускают корректное представления только конечных и конструктивных объектов, т.е. таких объектов, которые могут быть представлены последовательностью конечной длины, например, символами некоторого, также конечного, алфавита. Другое дело - работа с числами. Так, в машине нельзя представить не только все трансцендентные и иррациональные числа, но даже все рациональные числа. Даже действительные числа могут быть представлено только конечным, хотя может быть и очень большим подмножеством. Обычно. в ЭВМ используются две формы представления чисел: в форме чисел с плавающей запятой и в форме чисел с фиксированной запятой. Над этими множествами чисел и определен класс вычислений, называемый - "вычислениями на ЭВМ". Некорректные в общем случае, и неубедительные с точки зрения "чистой математики", эти вычисления обеспечивают интерпретацию отображений общих представлений об окружающем мире в виде конкретных моделей и алгоритмов.
В общем случае, источники ошибок, возникающие при решении вычислительных задач на ЭВМ, можно условно разделить на:
А. Ошибки, возникающие при формализации задач предметной области
  • неточности используемой математической модели
  • погрешности используемого приближенного метода
  • устойчивость вычислительных алгоритмов
B. Ошибки вычислительных операций
  • ошибки в исходных данных
  • погрешности выполнения арифметических операций
  • погрешности округления при вычислениях
С. Погрешности параллельных вычислений
  • погрешности параллельного аппаратного оборудования
  • погрешности программного обеспечения

2. Погрешности исходных данных

Исходные данные, выражающие, обычно значения тех или иных физических величин, носят приближенный характер, содержат ошибку измерения. Число x, которое принимается принимается за истинным значением величины X, называется ее приближенным значением, абсолютная величина разности между ними - "абсолютная погрешность" (ошибку): Dx = |X-x|. Так как истинное значение x обычно неизвестно, под Dx понимается предельная абсолютная погрешность, число, не меньшее абсолютной погрешности. "Относительная погрешность" данных определяется как величина: d =Dx/X. Для приближенного числа pi = 3.14 предельная абсолютная погрешность равна 0.0016, относительная - 0.00051 или 0.05%. Абсолютная погрешность приближенного числа определяется числом десятичных знаков в его записи, относительная - количеством значащих цифр в числе. (Значащими цифрами являются все цифры числа, кроме тех нулей, которые употребляются для определения места других цифр в десятичной дроби или же замещают неизвестные цифры. Так, в записи 0.040 значащими цифрами являются 40)

3. Формы представления чисел в ЭВМ

Формы представления чисел в ЭВМ определяются правилами аппроксимации действительных чисел посредством конечных множеств машинных представлений. Формы представления чисел могут отличаться диапазоном и точностью представления чисел. Диапазон представления характеризуется самым большим (M) и самым малым (m) по абсолютной величине из машинных чисел. (В качестве диапазона, иногда, рассматривают величину: v = log(b) (M/m), где b - основание логарифма, равное базе системы счисления. Здесь предполагается использование позиционной системы счисления для представления чисел. ). Точность представления оценивается абсолютной разностью машинных чисел: D = |a2-a1| и относительной: d = |(a2-a1)/a1|, где a1,a2 - значения двух последовательных машинных чисел, таких, что a1!=0 и |a2|>|a1|.

4. Числа с плавающей запятой

Подмножество вещественных чисел, которое может быть представлено в ЭВМ в форме чисел с плавающей запятой, принято обозначать буквой F и определять его элементы для конкретной архитектуры - "машинные числа", (по Форсайту и др.) четырьмя целочисленными параметрами: базой b, точностью t и интервалом значений показателя [L,U]. Множество F содержит число нуль и все f числа вида: f = (+/-).d1d2...dt * be, где е называется показателем, число .d1d2...dt = (d1/b+ ....+dt/(bt)) - дробной частью - мантиссой, причем: 0<=di<b, L<=e<=U. Каноническая или нормализованная форма F определяется дополнительным соотношением d1 != 0 ; это условие позволяет устранить неоднозначность представления одинаковых чисел, дает наивысшую возможную точность представления чисел. Особенности F:
Вывод оценок
Оценка минимального абсолютного значения числа множества F.
min(f= .d1d2...dt * be) = min(f= (.d1d2...dt) * (be)) = min(.d1d2...dt) * min(be) = .00...1 * bL = 1 * b-t * bL = bL-t
Для нормализованных чисел последние выкладки выглядят так:
.... = .10...0 * bL = 1 * b-1 * bL = bL-1

4.1.Представление систем с плавающей точкой в ЭВМ

Примеры систем с плавающей точкой для различных ЭВМ:
Cray-1 - [2,48,-16384,8191],
Сетунь - [3,18,-121,121],
IBM-360/370 - [16,6,-64,63] - обычная точность
IBM-360/370 - [16,14,-64,63] - двойная точность
Обычно, в ЭВМ реализуется работа над числами с плавающей запятой в двух форматах последних: стандартным и длинным (с двойной точностью). Форматы чисел второй группы,обычно, имеют мантиссы с удвоенным числом разрядов, реже, для них отводится и другое, большее число разрядов и для представления порядков. Мантиссы обычно представляется в виде последовательности цифр: d1d2..dt*bt со знаком.
Форматы действительных чисел с плавающей запятой для i486: стандартное (одинарной точности - Single Precision) 32 битовое слово состоит из 23 битовой мантиссы (и знака) и 8 битового порядка, двойное (двойной точности - Double Precision): 64 ~ 52 и 11. Процессоры класса выше 486 могут обрабатывать числа повышенной точности - Extended Precision: 80 ~ 64 бита мантиссы, 15 бит порядок.
Хотя целые числа являются подмножеством вещественных чисел, машинные числа с плавающей и фиксированной запятой рассматриваются как непересекающиеся типы. Вопрос совместного использования этих представлений достаточно деликатен, так как каждое из них имеет свой формат представления и раздельный аппарат (команды) работы с ними. При этом, например, результат интерпретации Фортран-конструкции: (7 .EQ. 7.0) не обязательно 'истинно' (.EQ.- это операция отношения "равно"). Обычно, реализуется набор программно-аппаратных средств преобразования данных.

5. Особенности арифметики с плавающей запятой

В общем случае, арифметические операции над элементами дискретного подмножества вещественных чисел F не корректны.
Результат арифметических операций чисел с плавающей запятой может: Поэтому, на множестве чисел с плавающей запятой определяются и "плавающие" арифметические операции, за результаты которых, если они не выходит за границы множества F, принимается ближайшие по значению элементы F. Точная семантика этих машинных операций, реакция на переполнение определяется схемой арифметико-логического устройства ЭВМ, для унификации предложен стандарт: IEEE-754. Плавающие вычисления отличаются от арифметических и это необходимо учитывать при работе на ЭВМ. Н. Вирт (Вирт Н. Систематическое программирование.Введение. М- Мир, 1977- 183с.), в частности, показал, что в плавающей арифметике нарушаются законы ассоциативности и дистрибутивности.

5.1. Примеры нарушения алгебраических законов в плавающей арифметике

Примеры из четырехразрядной десятичной арифметики по Н. Вирту.
А) Пусть x=9.900 y=1.000 z=-0.999 и тогда:
  1. (x+y)+z = 9.910
  2. x+(y+z) = 9.901
В) Пусть x=1100. y=-5.000 z=5.001 и тогда:
  1. (x*y)+(x*z) = 1.000
  2. x*(y+z) = 1.100
Здесь операции + и * - плавающие машинные операции. Такие 'чиcленные' процессы называют иногда 'неточными', здесь нарушаются ассоциативный и дистрибутивный законы арифметики..
Поучительный пример соотношений различных форм представлений чисел с плавающей запятой приводит Дж. Форсайт (Дж. Форсайт, М. Малькольм, К. Моулер, Машинные методы математических вычислений. М. Мир 1980 - 277 с.).
Используя десятичное число 0.1 в качестве величины шага в неком виртуальном алгоритме можно расчитывать, что десять таких шагов будет равен одному шагу длины 1.0. Но это не так, на множестве F при b=2, ибо 1/10 не имеет конечного разложения по степеням 1/2:
1/10 = 0/21 + 0/22 + 0/23 + 1/24 + 1/25 + 1/26 + ... Итак, десятичное число 0.1 равно двоичному: 0.0001100011000..., или восьмеричному: 0.063146314..., или шестнадцатиричному: 0.199999.... Если зафиксировать в дробях t разрядов точности, то сложение десяти таких дробей не даст в результате точной единицы.

6. Точность плавающей арифметики

Точность плавающей арифметики можно характеризовать посредством
машинного эпсилона. Максимальное число Е такое, что 1.+ Е = 1. является мерой точности представления чисел на данной ЭВМ (машинное эпсилон). Грубая схема вычисления эпсилона:
     EPS = 1.0
  1  EPS = 0.5 * EPS
     EPS1 = EPS + 1.0
     IF (EPS1 .GT. 1.0) GO TO 1
Задача. Написать программу, определяющую количество разрядов, используемых для представления мантиссы чисел с плавающей запятой. (Пусть на испытываемой ЭВМ мантисса числа хранится в нормализованном виде 1A2A3...An).

7. Погрешности округления при вычислениях

При выполнении арифметических операций на ЭВМ могут появляться числа с большим количеством значащих цифр, чем у исходных чисел. и большим, чем может быть представимо на данной ЭВМ. Например, при умножении, число значащих цифр может удвоиться. Поэтому необходимо проводить округление результата вычислений.
Простейший способ округления состоит в отбрасывания младших разрядов. Пусть при вычиcлении получен неотрицательный результат (q-ичная дробь): X = 0.An An-1 An-2....As As-1 As-2....A1 .
Тогда,округленное число есть: Xs= 0.An An-1An-2....As , а ошибка округления (абсолютная):Dx=|Xs-X|<=q-[n-(s-1)], где n-(s-1) число значащих цифр в Xs. Максимально возможная относительная ошибка при этом будет равна: dx = Dx/X <= (q-[n-(s-1)])/(1/q) = q-(n-s).
Другим общепринятым способом округления считается "симметричное" округление, при котором:
Xs = Xs + q-(n-s), если |Xs-X| >= 0.5*q-(n-s) и Xs = Xs, если |Xs-X| < 0.5*q-(n-s)
При этом способе округления максимально возможное значение относительной ошибки: dx <= 0.5*q-(n-s)+1.

7.1. Погрешности выполнения арифметических операций

Формулы оценки абсолютной и относительной погрешности арифметических операций.
  1. Сложение: X=X1+X2, X1>0 , X2>0
    Абсолютная погрешность: Dx = Dx1 + Dx2 , относительная: dx = dx1 + dx2
  2. Вычитание: X=X1-X2, X1>X2>0
    Абсолютная погрешность:Dx=Dx1 + Dx2,
    относительная:dx=(X1*dx1 + X2*dx2)/X
    Если X1 >> (много больше) X2, то dx (почти равно) dx1.
    Если X1 (почти равно) X2, то dx будет очень велико. При вычитании близких по величине чисел получается большая потеря верных знаков.
  3. Произведение: X = X1X2
    Погрешности: Dx = (X/X1)*Dx1 + (X/X2)*Dx2, dx = dx1dx2
  4. Частное двух величин: X = X1/X2
    Погрешности:Dx = (|X2|*Dx1 + |X1|*Dx2)/X22, dx = dx1dx2
Формулы дают выражения для определения ошибки результата каждого из 4 арифметических действий как функции от величины чисел, участвующих в операциях, и абсолютных ошибок для них (например, погрешностей исходных данных). Для определения полной ошибки результата нужно к этим ошибкам отдельно добавить ошибки округления. Пример расчета полной ошибки для суммирования положительных чисел (Г.К. Боровин).
Формула полной ошибки для суммирования положительных чисел Ai(i=1,..,n) имеет вид:
Ds = A1*da1 + A2*da2 +...+ An*dan + d1*(A1+A2) +..+ dn-1*(A1+..+An) + dn , где
dai - относительные ошибки представления чисел в ЭВМ, а di - относительные ошибки округления чисел при каждой следующей операции сложения. Пусть: все dai = da и di = d , a Ks = A1+A2+..+An, тогда:
Ds = da*Ks + d*[(n-1)*A1+(n-1)*A2 +...+ 2*An-1 + An]
Очевидно, что наибольший "вклад" в сумму ошибок вносят числа, суммируемые вначале. Следовательно, если суммируемые положительные числа упорядочить по возрастанию, максимально возможная ошибка суммы будет минимальной. Изменяя порядок суммирования чисел можно получать различные результаты. Но если даже слагаемые отличаются друг от друга незначительно, на точность результата может оказать влияние способ суммирования. Пусть суммируются 15 положительных чисел, тогда ошибка результата: Ds = da*Ks + d*(14*A1+14*A2+13*A3+....+2*A14+A15). Слагаемое da*Ks не зависит от способа суммирования, и далее не учитывается. Пусть слагаемые имеют вид: Ai = А0+ei, где i=1,...,15, тогда:
Dss = 199*(A0+em)*d, где em = max(ei), d - ошибка округления при выполнении арифметической операции сложения. Если провести суммирование этих чисел по группам (три группы по четыре числа и одна группа из трех чисел), то ошибки частных сумм имеют вид:
Ds1 = d*(3*A1+3*A2+2*A3+A4) <= 9*d*(A0+em)
Ds2 = d*(3*A5+3*A6+2*A7+A8) <= 9*d*(A0+em)
Ds3 = d*(3*A9+3*A10+2*A11+A12) <= 9*d*(A0+em)
Ds4 = d*(3*A13+2*A14+A15) <= 5*d*(A0+em)
а полная оценка ошибок округления будет Ds <= 32*d*(A0+em), что меньше Dss. Итак суммирование по группам дает меньшую ошибку результата.
Например, разделив процесс суммирования массива положительных чисел на параллельные процессы, и затем получив сумму частных сумм, можно получить результат, отличный от последовательного суммирования в одном процесс.

8. Оптимизационные преобразования программ

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

1. Балансировка дерева вычислений

Балансировка дерева вычислений (tree-height reduction or balancing) выражений позволяют использовать конвейерное АУ без пропуска рабочих тактов. Так, вычисление суммы вещественных чисел:
A+B+C+D+E+F+G+H, будет запрограммировано как последовательность операций:
(((A+B)+(C+D))+((E+F)+(G+H))); это нарушает заданную по-умолчанию последовательность вычислений с накоплением одной частной суммы и может повлиять на результат.

2. Исключение общих подвыражений

Алгоритмы исключения общих подвыражений (Common subexpession elimination) также могут изменить порядок вычислений.
Если компилятор распознает в выражениях повторяющееся вычисление, то это вычисление производятся один раз, его результат сохраняется на регистре, и в дальнейшем используется этот регистр. Тем самым исключается избыточность вычислений.
X = A + B + C + D ----> REG = B + C
Y = B + E + C X = A + D + REG
Y = E + REG

3. Разворачивание циклов

Разворачивание циклов (loop unrolling) - расписывание цикла последовательностью операторов присваивания: либо полностью, либо размножение тела цикла с некоторым коэффициентом (фактором) размножения.
Производится частичное или полное разворачивание цикла в последовательный участок кода. При частичном разворачивании используется так называемый фактор разворачивания (который можно задавать в директиве компилятору).
     DO I=1,100                  DO I=1,100,4
       A(I) = B(I) + C(I)          A(I) = B(I) + C(I)
     ENDDO                         A(I+1) = B(I+1) + C(I+1)
                                   A(I+2) = B(I+2) + C(I+2)
                                   A(I+3) = B(I+3) + C(I+3)
                                 ENDDO
При этом преобразовании снижается количество анализов итерационной переменной. Данный алгоритм также может привести к нарушению предписанного первоначально порядка вычислений. Например:
     DO I=1,10        DO I=1,10,2
     S = S + A(I)       S = S + A(I)
     ENDDO             S1 = S1 + A(A+1)
                      ENDDO
                      S = S + S1
Здесь, суммирование проводится отдельно для четных и нечетных элементов с последующем сложением частных сумм.

9. Блокировка оптимизационных преобразований

Для блокировки оптимизационных преобразования в компиляторах, обычно, предусматриваются опции (режимы), "флаги" , которые можно использовать при получении объектного кода. Некоторые предупредительные указания предусматриваются в стиле написания программ, так, в Фортране 90 концепция сохранения скобок оговорена стандартом языка.

10. Особенности вычисления выражений в Фортране 90 (стандарт)

7.1.7. Выполнение операций

"При вычислении выражения А+В, где А,В - массивы, процессор может выполнять по-элементные операции сложения в любом порядке"
7.1.7.1. Вычисление операндов

"Процессор не обязательно вычисляет все операнды выражения или полностью каждый операнд, если значения можно определить без этого. Так в выражении: X .GT. Y .OR. L(Z) не надо вычислять функцию L(Z), если X больше Y."
7.1.7.2. Целостность скобок

"Для выражения: A+(B+C) процессор на должен вычислять математически эквивалентное выражения: (A+B)+C."
7.1.7.3. Выполнение числовых встроенных операций

"... процессор может вычислять любое математически эквивалентное выражение при условии, что целостность скобок не нарушается."
В стандарте языка Фортран 90 приводятся допустимые и не допустимые альтернативных формы преобразования выражений:
Для целых i,j;вещественных и комплексных a,b,c; произвольных x,y,z
разрешаются преобразования: и запрещаются:
x+y -> y+x i/2 -> 0.5 *i
x*y -> y*x x*i/j -> x*(i/j)
-x+y -> y-x (x+y)+z -> (x+y+z)
x+z+y -> y+(x+z) i/j/a -> i/(j*a)
x*a/z -> x*(a/z) (x+y)+z -> x+(y+z)
a/b/c -> a/(b*c) (x*y)-(x*z)-> x*(y-z)
a/5.0 -> 0.2*a x*(y-z) -> x*y-x*z

11. Погрешности асинхронного распараллеливания

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

Литература
  1. Форсайт Дж., Малькольм М., Молер К. Машинные методы математических вычислений.-М.:Мир,1980
  2. Голуб Дж., Ван Лоун Ч., Матричные вычисления.-М.:Мир,1999.-548с.
  3. Боровин Г.К., Комаров М.М., Ярошевский В.С. Ошибки-ловушки при программировании на фортране. М.: Наука. 1987.-144с.
  4. Вирт Н. Систематическое программирование. Введение.М.: Мир.1977