Дж. Деммель - Вычислительная линейная алгебра, страница 20

С другой стороны, можно проверить, что число операций с плавающей точкой, выполняемых в строке (1) алгоритма более медленными операциями уровня 2 и 1, составляет для малых Ь примерно нгЬ/2; это число растет с ростом Ь, поэтому не стоит брать слишком большое значение для Ь. Оптимальное значение Ь зависит от машины и может быть найдено для каждого конкретного типа компьютеров.

Часто используют значения Ь = 32 и Ь = 64. Тщательными реализациями алгоритмов 2.9 и 2.10 являются ЬАРАСК- подпрограммы эйест2 и вйесгГ (см. ХЕТЫВ/1арас1«). Дальнейшую информацию о блочных алгоритмах, включая подробные сведения о производительности на ряде различных компьютеров, можно найти в 110] (см. также конспекты лекций на РАВ.А1 ЬЕ1, НОМЕРАОЕ).

2.6.4. Еще раз о параллелизме и других проблемах производительности В этом разделе мы дадим краткий обзор прочих проблем, связанных с максимально эффективной реализацией гауссова исключения (и других алгоритмов линейной алгебры). Параллельный компьютер имеет р(> 1) процессоров, которые можно единовременно использовать для решения одной и той же задачи. Можно надеяться решить любую задачу на таком компьютере в р раз быстрее, чем на обычной однопроцессорной машине.

Однако подобная «идеальная эффективность» достигается редко даже в том случае, когда в любой момент имеется р независимых подзадач; это происходит из-за накладных расходов, связанных с координацией процессов и рассылкой данных из хранящих их процессоров тем, которые в этих данных нуждаются. Это последнее обстоятельство есть еще один пример иерархически организованной памяти: с точки зрения процессора », его собственная память быстра; однако получение данных из памяти другого процессора у' может иногда происходить медленнее в тысячи раз. Гауссово исключение содержит в себе много возможностей для параллелизации, твк как на каждом шаге элементы оставшейся подматрицы могут перевычисляться параллельно и независимо друг от друга.

Однако нужна некоторая осмотрительность, чтобы достигнуть высокой производительности. В этой области имеется два стандартных программных продукта. Описанная в предыдущем разделе ЬАРАСК-программа вбесгг (см. 110]) предназначена для параллельп»лх машин с разделлвмвй памятью в предположении, что они оснащены параллельными реализациями пакета ВЬАБ.

Родственная библиотека, называемая ЯсаЬАРАСК (от Яса1аЬ1е ЬАРАСК [34, 53]), была разработана для параллельнъгх машин с распределенной памятью, т.е. машин, требующих специальных операций для обмена данными между различными процессорами. Все программы доступны в ХЕТЫВ'е в поддиректориях ЬАРАСК и ЯсаЬАРАСК. Более подробно ЯсаЬАРАСК описан в конспектах лекций, помещенных на РАНАЬЬЕЬ НОМЕРАОЕ.

Обширные данные о производительности программ решения линейных систем содержатся в техническом отчете 86 Глава 2. Решение линейных уравнений ЫЬ!РАСК ВепсЬшаг)с [85), а более современную версию обзора производительности можно найти по адресу 5!ЕТЫВ/ЬепсЬшаг1с/рег1огшапсе.ра; см. также РегГогшапсе ЫаЬаЬазе Яегпег'.

Па май 1997 г. наибольшая скорость решения линейной системы посредством гауссова исключения была достигнута для задачи размерности и = 215000 на машине 1пФе1 АБС1 Ориоп Нег! с р = 7264 процессорами. Задача решалась со скоростью, превышающей 1068 гигафлопов при пиковой производительности компьютера 1453 гигафлопа. Встречаются матрицы настолько большие, что они не могут поместиться в основной памяти никакого из имеющихся компьютеров. Такие матрицы хранят на диске и в ходе гауссова исключения считывают в основную память частями.

Конструкция соответствующих программ весьма схожа с описанными выше алгоритмами; программы этого типа тоже включены в ЯсаЬАРАСК. Можно надеяться, что трансляторы станут в конце концов настолько умными, что, исходя из простейшей реализации гауссова исключения посредством трех вложенных циклов, будут способны автоматически «оптимизировать» программу до уровня блочного алгоритма, обсуждавшегося в предыдущем разделе. Хотя в этом направлении в настоящее время ведется большая работа 1см. библиографию в недавно изданном учебнике по трансляторам 12641), все еще не существует гарантированно быстрой альтернативы использованию оптимизированных библиотек типа ЬАРАСК и ВсаЬАРАСК.

2.7. Специальные линейные системы Как уже говорилось в разделе 1.2, важно уметь использовать любую специальную структуру матрицы, чтобы ускорить решение задачи и уменьшить расход памяти. Разумеется, на практике нужно принимать в учет стоимость дополнительных программистских усилий, затрачиваемых на эксплуатацию этой структуры. Пусть, например, нашей единственной целью является минимизация времени вычисления решения, и пусть нужна неделя дополнительной программистской работы, чтобы уменьшить время решения задачи с 10 секунд до одной.

Это имеет смысл делать лишь в том случае, если мы собираемся пользоваться полученной программой более чем (1 неделя * 7 дней/в неделе * 24 час/в день * 3600 сек/в часе) /(10 сек " — 1 сек) = 67200 раз. К счастью, существуют достаточно часто встречающиеся специальные структуры, для которых имеются стандартные методы решения, и мы, конечно, должны использовать зти методы. Мы рассмотрим здесь следующие типы специальных матриц: 1. симметричные положительно определенные матрицы (с.п.о, матрицы); 2. симметричные незнакоопределенные матрицы; 3.

ленточные матрицы; 4. разреженные матрицы общего вида; 5. плотные матрицы, зависящие от менее чем из независимых параметров. Будут рассматриваться лишь вещественные матрицы; обобщения методов на случай комплексных матриц очевидны. г Ыср://регГогшапсе.пес!1Ь.огх/регГогшапсе/Ыш1/Р1эзсор.мш! 2.7. Специальные линейные системы 87 2.7.1. Вещественные симметричные положительные определенные матрицы Напомним, что вещественная матрица А называется с.п.о матрицей, если А = Ат и хтАх > О для всех х ~ О.

В этом разделе мы покажем, как решить систему Ах = Ь с с.п.о. матрицей А за половину времени и с использованием вдвое меньшей памяти по сравнению с гауссовым исключением. Предложение 2.2. 1. Если Х вЂ” невыроэкденная матрица, то А тогда и только тогда является с.п.о. матрицей, когда ХтАХ есть с.п.о. матрица.

2 Если А — с.п.о. матрица, а Н вЂ” произвольная главная подматрица в А (т. е. Н = А(у: 7с,у': й) для некоторого г' < )г), то Н также с.п.о. матрица. 3 А тогда и пголько тогда является с.п.о. матрицей, когда А = Ат и все собственньге значения матрицы А положительны. 4 Если А — с.п.о. матрица, то все ан > О и шах,. ~аг ~ = шах, ан > О. 5 А тогда и только тогда является с.п.о. матрицей, когда существует единственная невырожденпая нижетреугольная матрица Е с положительными диагональными элементами, такая, что А = ЕЬт. Представление А = Е1т называется разложением Холесского, а Š— множителем Холесского матрицы А. Доказательство. 1. Невырожденность матрицы Х означает, что Хх ~ О для всех х ~ О, поэтому хтХтАХх > О, если х ф О. Таким образом, если А — с.п.о.

матрица, то и ХтАх — с.п.о. матрица. Обратная импликация доказывается с помощью матрицы Х 2. Предположим вначале, что Н = А(1: т,1: т). Тогда для любого твектора у вектор х = )у~,О) размерности п удовлетворяет соотношению утНу = хтАх. Поэтому, если хтАх > О для любого ненулевого вектора х, то утНу > О для любого ненулевого вектора у; следовательно, Н есть с.п.о. матрица. Если Н не является ведущей главной подматрицей, то пусть Р— перестановка, переводящая Н в левый верхний угол матрицы РтАР; к последней применяем утверждение 1. 3. Пусть Х вЂ” вещественная ортогональная матрица, составленная из собственных векторов матрицы А, тогда ХтАХ = Л есть диагональная матрица, содержащая на диагонали (вещественные) собственные значения А,. Поскольку хтЛх = 2,'г Л;хг, то Л тогда и только тогда является с.п.о. матрицей, когда все Лг > О. Остается применить утверждение 1.

4. Пусть е, — столбец с номером г' единичной матрицы. Тогда етАе, = ан > О для всех 1. Если ~аы~ = шах, (агг(, но )г ~ 1, положим х = ег — зйп(ам)еь Тогда хтАх = ам + аа — 2~ам ~ < О, что противоречит положительной определенности матрицы А. 5.

Предположим, что А = ЕЬ', где Š— невырожденная матрица. Тогда хтАх = (хг Ь)(Ьтх) = ()г" тхЯ > О для всех х ф О, поэтому А — с.п.о. матрица. Если А — с.п.о. матрица, то мы докажем существование (треугольной) матрицы Е индукцией по порядку и. Если в нашем построении выбирать 88 Глава 2. Решение линейных уравнений положительное значение для каждого 1и, то оно определит Ь единственным образом. При и ='1 полагаем 1ы = 1/аы, что возможно, поскольку ам > О. Как и в случае гауссова исключения, достаточно разобраться с блочными 2 х 2-матрицами.

Имеем аы Агг '=::.;][. а[-: -', ~ 12 22 ап А= [ Ат А таким образом, (и — 1) х (и — 1) матрица А22 —— А22 — — ',2 — ™ — симметрична. О Согласно утверждению 1, ~ - есть с.п.о. матрица, но тогда, согласно 22 утверждению 2, и А22 есть с.п.о. матрица. По предположению индукции, т существует (треугольная) матрица б, такая, что А22 = Х,Ь и ,/аы О ~ Ат о'а~ ~ а=[ .'й1[:" ~! с на А22- ~1 и 11 /аам — т ат Т О Ь 1/ам О А 212 /ап Мы можем оформить это индуктивное рассуждение в виде следующего алгоритма.

Алгоритм 2.11. Алгоритлч Холесского: /от 1' = 1 2о 11 111 = (а11 — 2„1 121)112 /от1=1+1 2оп 1„. = (п,1 — 2 ~1, 1121 ь)/1 епЫ /от ЕП11 /От г 21+ ~ ~' 21 = - г+ о(пг), 1=1 1=1+1 Если А не является положительно определенной матрицей, то (в точной арифметике) алгоритм прекратит работу досрочно при попытке извлечь квадратный корень из отрицательного числа либо разделить на нуль; это наблюдение указывает самый экономичный способ проверки свойства положительной определенности симметричной матрицы. Как и в случае гауссова исключения, Ь может быть записана на место нижней половины матрицы А. В алгоритме происходят обращения только к элементам из нижней половины А, поэтому, в действительности, достаточно памяти п(п+ 1)/2 вместо пг.