Программирование в Scilab Микаэль Боден Micha

• Название: Программирование в Scilab Микаэль Боден Micha
• Дата: 18.09.2022
• Размер: 1.35 Mb.
• Имя файла: progscilab-v.0.10_ru.pdf
• Тип: Реферат; #683094
• страница: 13 из 13

• BLAS уровня 2 выполняет матрично-векторные операции.
• BLAS уровня 3 выполняет матрично-матричные операции.
Библиотека LAPACK (Linear Algebra Package — пакет линейной алгебры) [
3
] является коллекцией процедур Фортрана для решения задач линейной алгебры высокого уровня. Это включает в себя решение систем линейных уравнений, решение линейных систем по мето- ду наименьшего квадрата, задачи нахождения собственных значений и задачи нахождения сингулярного значения.
Библиотека ATLAS (Automatically Tuned Linear Algebra Software — автоматически на- страиваемое программное обеспечение линейной алгебры) [
1
] предоставляет оптимизирован- ные процедуры линейной алгебры BLAS и небольшой набор процедур, оптимизированных библиотекой LAPACK. Библиотека концентрируется на применении эмпирических методов для обеспечения характеристики переносимости. Библиотека ATLAS предоставляет во мно- гих случаях чувствительное улучшение производительности по сравнению с системами «ос- новных» BLAS/LAPACK.
128

Библиотека математического ядра фирмы Intel (Intel Math Kernel Library) [
24
] явля- ется библиотекой высокооптимизированных, широко разветвлённых математических проце- дур для научно-технических и финансовых приложений, которые требуют максимальной производительности. В данный момент под Windows Scilab использует только процедуры линейной алгебры от MKL взамен BLAS.
В Windows следует выбирать Intel MKL для установки (это сделано по умолчанию для Scilab версии 5). В Linux-системах следует использовать оптимизированную версию библиотеки ATLAS, если производительность играет роль. Эти вопросы рассмотрены в раз- делах
5.4.3
для Windows и
5.4.4
для Linux.
Библиотеки ATLAS и Intel MKL используют разные методы низкоуровневой оптимиза- ции для получения производительности при вычислениях с плавающей запятой. Этот вопрос ещё раз рассмотрен в следующем разделе.
5.4.2
Методы низкоуровневой оптимизации
В этом разделе кратко представлены методы низкоуровневой оптимизации, использу- емые в численных библиотеках линейной алгебры, таких, как ATLAS и Intel MKL. Пред- ставлены различные уровни памяти в компьютере. Рассмотрены алгоритмы, используемые в ATLAS для улучшения использования различных уровней кэша. Кратко представлены на- боры инструкций процессоров и используемые оптимизированными библиотеками линейной алгебры.
Главной причиной улучшения производительности векторизованных функций являет- ся то, что Scilab способен выполнять вычисления над коллекцией данных вместо обработки каждого вычисления отдельно. Действительно, производительность вычисления значитель- но зависит от того, насколько ЦП может обращаться к данным, которые доступны в раз- личных типах памяти. От наиболее быстрой до наиболее медленной памяти, используемой в современных компьютерах, — это различные уровни кэша (встроенного в сам процессор)
и оперативная память (RAM) (здесь мы не рассматриваем жёсткий диск в качестве запо- минающего устройства, хотя он может быть включён в управление виртуальной памятью).
Рисунок
38
представляет различные уровни памяти, обычно используемой в компьютере.
Регистр процессора
Кеш уровня 1
Кеш уровня 2
Кеш уровня 3
RAM
Рис. 38: Уровни памяти в компьютерной системе.
Теперь рассмотрим следующий пример, где мы выполняем перемножение матрицы на матрицу.
A = r a n d ( 1 0 , 1 0 ) ;
B = r a n d ( 1 0 , 1 0 ) ;
C = A * B ;
В соответствующем шлюзе Scilab использует процедуру BLAS DGEMM для того, чтобы выполнить требуемое вычисление. Если пользователь выбрал библиотеки ATLAS или MKL
во время инсталляции, то вычисления производятся оптимизированным способом.
Поскольку исходный код Intel MKL не является открытым, то мы не знаем какие точно методы использованы в этой библиотеке. В отличие от него, ATLAS является библиотекой с открытым исходным кодом, и мы можем подробно анализировать его.
129

Подход ATLAS [
52
] состоит в изоляции машинно-зависимых параметров процедур, ко- торые имеют дело с выполнением оптимизированного матричного умножения на кристалле с кэшем (т. е. кэшем уровня 1). Перемножение на кристалле автоматически создаётся генера- тором кода, который использует измерения времени для определения факторов правильной блокировки и развёртывания циклов для выполнения оптимизированного умножения на кристалле.
Матричное умножение раскладывается на две части:
1. высокоуровневое, общего размера умножение вне кристалла, которое платформо-незави- симо,
2. низкоуровневое, фиксированного размера умножение на кристалле, которое машинно- зависимое.
Алгоритм высокого уровня основан на умножении блочных матриц. Предположим, что
A — это матрица размером m × n, а B — это матрица размером k × n. Предположим, что целое число NB, количество блоков, делится на m, n и k. Алгоритма умножения блочных матриц мог бы быть написан на языке Scilab следующим образом.
f u n c t i o n C = offchip - b l o c k m a t m u l ( A , B , NB )
[ m , k ] = s i z e ( A )
[ k , n ] = s i z e ( B )
C = z e r o s ( m , n )
for ii = 1: NB : m
I = ii : ii + NB -1
for jj = 1: NB : n
J = jj : jj + NB -1
for kk = 1: NB : k
K = kk : kk + NB -1
C ( I , J ) = C ( I , J ) + A ( I , K )* B ( K , J ) // on - c h i p end end end e n d f u n c t i o n
В реализации ATLAS инструкция C(I, J) = C(I, J) + A(I,K)*B(K, J) выполнена с по- мощью функции умножения на кристалле.
Главным параметром алгоритма вне кристалла является NB, но и другие факторы при- нимаются во внимание. Действительно, ATLAS может использовать другие реализации это- го алгоритма, основанного, например, на другом упорядочивании трёх вложенных циклов.
Во всех этих реализациях цикл по k является всегда самым внутренним циклом. А внешний цикл может быть по m (по строкам в A) и по n (по столбцам в B). Для того, чтобы найти оп- тимальное значение параметров (включая NB, порядок циклов и другие параметры), ATLAS
выполняет автоматизированный эвристический поиск через замеры времени. Для каждого набора исследуемых параметров ATLAS использует генератор кода и измеряет производи- тельность с этими установками. В конечном счёте сохраняется наилучшая реализация и генерируется соответствующий исходный код, затем компилируется.
Умножение на кристалле чувствительно к нескольким факторам, включающим повтор- ное использование кэша, переполнение кэша инструкций, порядок инструкций с плавающей запятой, максимальное количество циклов, проявление возможного распараллеливания и ошибки кэша. Более того, умножение на кристалле может использовать специальные ин- струкции, доступные процессору, такие, как инструкции множественных данных одной ин- струкции (SIMD). Эти инструкции позволяют выполнять ту же работу над множественными
130

данными одновременно и встроить распараллеливание данных. Есть несколько наборов ин- струкций SIMD, включающих
• MMX, разработанный фирмой Intel в 1996,
• 3DNow!, разработанный фирмой AMD в 1998,
• Расширение потоковой SIMD (SSE), разработанное фирмой Intel в 1999,
• SSE2, разработанный фирмой Intel в 2001,
• SSE3, разработанный фирмой Intel в 2004,
• Продвинутые векторные расширения, будущее расширение, предложенное фирмой Intel в 2008.
Например,
инструкция
MULPD
может быть использована для выполнения
SIMD-умножения двух пакетных значений с плавающей запятой двойной точности. Похо- жей инструкцией является ADDPD, которая выполняет суммирование SIMD двух пакетных значений с плавающей запятой двойной точности.
На практике формирование бинарного файла для библиотеки ATLAS может потребо- вать несколько часов времени ЦП перед тем, как оптимизированные настройки автомати- чески идентифицируются системой ATLAS.
Эти низкоуровневые методы не является главной заботой пользователей Scilab. Одного это позволяет понимать почему специфичные библиотеки, поставляемые вместе с Scilab’ом,
являются разными для Pentium III, Pentium 4 или Dual или Quad Core. Более важно то, что это позволяет знать уровень оптимизации, который мы можем получить от вычислений с плавающей запятой и особенно линейной алгебры, которая является главной целью языка
Scilab.
5.4.3
Установка оптимизированных библиотек линейной алгебры для Scilab под
Windows
В этом разделе мы представляем метод установки оптимизированных библиотек ли- нейной алгебры для Scilab на Windows.
На Windows мы обычно устанавливаем Scilab после его скачивания с http://www.
scilab.org
(но некоторые пользователи имеют предустановленный Scilab на своих маши- нах). Когда мы скачиваем установщик Scilab для Windows и запускаем его, то мы можем выбирать из трёх возможностей:
• Полную установку (по умолчанию), где устанавливаются все модули,
• Установку, где некоторые модули отключены,
• Выборочную установку, где пользователь может выбирать включение или отключение модулей.
По умолчанию Scilab устанавливается с Intel MKL, но этот выбор может быть настроен пользователем в диалоге установки, который представлен на рисунке
39
В разделе «CPU Optimization for Scilab» диалога установки мы можем выбирать между тремя следующими пунктами.
131

Рис. 39: Выбор библиотеки линейной алгебры под Windows (Scilab v5.2.2). Мы можем вы- бирать между BLAS, ATLAS и Intel MKL.
• Download Intel Math Kernel Library for Scilab. Эта установка по умолчанию Эта биб- лиотека содержит оптимизированные библиотеки BLAS и набор LAPACK, предостав- ленный фирмой Intel. Эта библиотека является не полной Intel MKL, а только набором функций, которые используются в Scilab для BLAS и LAPACK.
• Atlas library for Windows. Эта библиотека скомпилирована консорциумом Scilab Con- sortium для обеспечения оптимальной производительности в широком диапазоне со- временных компьютеров.
• BLAS, LAPACK Reference library for Windows. Эта библиотека является опорной ре- ализацией, предоставленной http://www.netlib.org/blas и
http://www.netlib.org/
lapack
Эти опции имеют прямое влияние на две динамических библиотеки, которые хранятся в директории bin Scilab’а:
• bin \blasplus.dll : динамическая библиотека для BLAS,
• bin \lapack.dll : динамическая библиотека для LAPACK.
На практике можно иметь одну и ту же версию Scilab’а, установленную с разными биб- лиотеками линейной алгебры. Действительно, несложно настроить директорию установки
132

и затем выбрать конкретную библиотеку, которую желают использовать. Следовательно,
можно установить Scilab в три разные директории, в каждой свою конкретную библиотеку,
так, что мы можем сравнить их поведение и производительность. Этот метод используется в разделе
5.4.5
, где представлена разница производительности между этими тремя библио- теками в Windows на классическом тесте производительности.
5.4.4
Установка оптимизированных библиотек линейной алгебры для Scilab под
Linux
В этом разделе представлен метод установки оптимизированных библиотек линейной алгебры для Scilab под операционной системой Gnu/Linux.
Под Linux сложнее описать процесс установки, поскольку есть много разных дистри- бутивов Linux. Тем не менее есть два основных способа получения Scilab под Linux:
• скачать Scilab с http://www.scilab.org
,
• использовать систему пакетов дистрибутива Linux, например, Debian, Ubuntu, и т. д. . .
Какой бы источник ни был, Scilab идёт с Reference BLAS и LAPACK, собранными из http:
//www.netlib.org/blas и
http://www.netlib.org/lapack
Чтобы установить оптимизированную библиотеку ATLAS, мы должны собрать нашу собственную версию ATLAS. Поскольку это несколько сложно и занимает время ЦП (обычно один или два часа сборки), то это рассмотрено в конце этого раздела.
Более простое решение в том, чтобы использовать бинарные файлы из дистрибутива
Debian. Бинарные файлы ATLAS в Ubuntu приходят из дистрибутива Debian, так что можно так же получить бинарные файлы ATLAS на Ubuntu. Под Ubuntu мы можем использовать,
например, менеджер пакетов Synaptic для установки и удаления бинарных пакетов. В этом случае мы можем использовать бинарный файл libatlas3gf-base «ATLAS generic shared», который предоставляется Debian’ом.
Использование бинарного файла оптимизированной библиотеки ATLAS под Linux тре- бует несколько простых изменений вручную, которые мы собираемся описать. В директории scilab/lib/thirdparty мы найдём следующие файлы.
$ ls scilab - 5 . 3 . 0 / lib / t h i r d p a r t y
- rw - r - - r - - 1 mb mb
6.3 M 2010 -09 -23 1 4 : 4 5 l i b l a p a c k . so .3 gf .0
- rw - r - - r - - 1 mb mb
539 K 2010 -09 -23 1 4 : 4 5 l i b b l a s . so .3 gf .0
Когда мы удалим файл libblas.so.3gf.0 из директории установки Scilab’а, то Scilab будет использовать библиотеку BLAS из системы Linux. Следовательно, чтобы убедиться, что
Scilab использует библиотеку ATLAS, которую мы уже установили в систему, мы просто удаляем файл libblas.so.3gf.0 как в следующем примере.
$ rm scilab - 5 . 3 . 0 / lib / t h i r d p a r t y / l i b b l a s . so .3 gf .0
Для того, чтобы получить хорошую производительность, нам не следует использовать бинарный файл, предоставленный Debian’ом. Вместо этого, следует собрать свою собствен- ную библиотеку ATLAS на конкретную машину, которую мы используем. Это из-за обна- руженных во время сборки конкретных настроек, которые позволяют ATLAS’у улучшить свою производительность. Следовательно, двоичный файл ATLAS, который поставляется в
Debian разумно оптимизирован для машины, которую использовал для сборки заведующий пакетом. Он довольно хорошо работает и его легко использовать, но на другой машине не даёт максимальной производительности, которую можно получить.
133

К счастью, исходные коды доступны на http://math-atlas.sourceforge.net
, так что кто угодно в Linux’е может настроить собственную библиотеку ATLAS, чтобы получить хо- рошую производительность линейной алгебры. Более того, в Debian, процесс сделан гораздо легче, поскольку этот конкретный дистрибутив предоставляет все необходимые инструмен- ты. Более подробно эта тема рассматривается в разделе
5.6 5.4.5
Пример улучшения производительности
В этом разделе представлен простой эксперимент в Windows, который демонстрирует выгоду использования оптимизированной библиотеки линейной алгебры в Scilab’е. В каче- стве программы проверки производительности системы, мы рассматриваем перемножение двух квадратных, плотных, действительных матриц значений типа double.
Определим в следующем файле-сценарии функцию mymatmul, которая принимает раз- мер матриц n в качестве входного аргумента и выполняет перемножение двух матриц, за- полненных единичными элементами.
f u n c t i o n m y m a t m u l ( n )
A = o n e s ( n , n )
B = o n e s ( n , n )
C = A * B
e n d f u n c t i o n
Оператор умножения, в данном случае, связан с функцией библиотеки BLAS, которая на- зывается DGEMM. Эта функция часто используется для измерения производительности оптимизированных библиотек линейной алгебры.
Следующий набор команд использует функцию benchfun, которую мы уже представля- ли в разделе
5.1.4
. Мы устанавливаем размер стека в максимум, затем вычисляем умножение плотных матриц увеличивающегося размера.
s t a c k s i z e ( " max " )
b e n c h f u n ( " m a t m u l " , m y m a t m u l , l is t ( 1 0 0 ) , 0 , 10 );
b e n c h f u n ( " m a t m u l " , m y m a t m u l , l is t ( 2 0 0 ) , 0 , 10 );
b e n c h f u n ( " m a t m u l " , m y m a t m u l , l is t ( 4 0 0 ) , 0 , 10 );
b e n c h f u n ( " m a t m u l " , m y m a t m u l , l is t ( 1 0 0 0 ) , 0 , 10 );
Использовался Scilab версии 5.2.2 под Windows XP 32 бита. ЦП — AMD Athlon 3200+
на 2 ГГц, и система работает с 1 ГБ памяти. В следующем примере мы выполняем файл- сценарий программы проверки производительности с библиотекой Atlas.
- - > e x e c b e n c h _ m a t m u l . sce ;
m a t m u l : 10 i t e r a t i o n s , m e a n = 0 . 0 0 3 0 0 , min = 0 . 0 0 0 0 0 , max = 0 . 0 1 0 0 1
m a t m u l : 10 i t e r a t i o n s , m e a n = 0 . 0 3 3 0 4 , min = 0 . 0 2 0 0 2 , max = 0 . 0 8 0 1 1
m a t m u l : 10 i t e r a t i o n s , m e a n = 0 . 2 7 9 4 0 , min = 0 . 2 6 0 3 7 , max = 0 . 3 7 0 5 3
m a t m u l : 10 i t e r a t i o n s , m e a n = 4 . 1 1 8 9 2 , min = 4 . 0 8 5 8 7 , max = 4 . 1 9 6 0 3
Рисунок
40
представляет результаты, которые мы получили от трёх оптимизированных библиотек. Наилучшая производительность получена от Intel MKL, за которой близко сле- дует библиотека Atlas. Производительность библиотеки Reference BLAS-LAPACK очевидно ниже.
Используя оптимизированные библиотеки линейной алгебры, пользователи Scilab’а мо- гут получить быстрые вычисления линейной алгебры. На Windows Scilab поставляется с Intel
MKL, которая даёт в большинстве случаев великолепную производительность. Заметим, что эта библиотека коммерческая и предоставляется бесплатно пользователям Scilab, поскольку
134

Библиотека
N
Среднее (с)
REF-LAPACK
100 0,003004
REF-LAPACK
200 0,033048
REF-LAPACK
400 0,279402
REF-LAPACK
1000 4,118923
ATLAS
100 0,001001
ATLAS
200 0,011016
ATLAS
400 0,065094
ATLAS
1000 0,910309
Intel MKL
100 0,001001
Intel MKL
200 0,010014
Intel MKL
400 0,063091
Intel MKL
1000 0,834200
Рис. 40: Чувствительность производительности перемножения матриц в зависимости от биб- лиотеки линейной алгебры. Проверкой производительности является умножение матрицы на матрицу. Тест выполнялся в Scilab версии 5.2.2 под Windows XP 32 бита. ЦП — AMD
Athlon 3200+ на 2 ГГц, и система работает с 1 ГБ памяти.
консорциум Scilab Consortium имеет лицензию на Intel MKL. Под Linux Scilab поставляет- ся с библиотекой reference BLAS-LAPACK. Оптимизированные библиотеки Atlas доступны для из любого дистрибутива GNU/Linux. Более того, на GNU/Linux мы можем собрать свою собственную библиотеку Atlas и создать настроенную библиотеку, которая специально настроена на наш компьютер.
5.5
Измерение числа операций с плавающей запятой за секунду
(флопс)
Общая практика в информатике заключается в вычислении числа операций с пла- вающей запятой, которое может быть выполнено в течение одной секунды. Это даёт аб- бревиатуру flops, что означает FLoating Point Operations by Second (количество операций с плавающей запятой за секунду). В этом разделе мы представляем два стандартных те- ста производительность в Scilab, включающих в себя вычисление произведения двух плот- ных матриц и вычисления решения линейных уравнений с помощь LU-разложения. В обоих случаях мы сравниваем производительности библиотек линейной алгебры Reference BLAS,
ATLAS и Intel MKL, предоставляемых в Scilab под Windows.
5.5.1
Произведение матрицы на матрицу
В этом разделе мы рассматриваем производительность произведения матрицы на мат- рицу на машине с Windows, с несколькими библиотеками линейной алгебры.
Давайте рассмотрим произведение A и B, двух плотных матриц размером n × n, исполь- зуя оператор *. Число операций с плавающей запятой вычисляется напрямую по формуле:
C
ij
=
n
X
k=1
A
ik
B
kj
,
(3)
135

для i, j = 1, 2, . . . , n. Чтобы вычислить один элемент C
ij матрицы C, приходится выполнять n умножений и n суммирований. В результате, количество операций с плавающей запятой
2n. Поскольку количество элементов в матрице C n
2
, то общее число операций с плавающей запятой равно 2n
3
В этом случае Scilab использует процедуру DGEMM библиотеки BLAS, которая явля- ется частью BLAS уровня 3. Основы этого набора BLAS описаны в [
15
], где авторы подчёр- кивают, что общий руководящий принцип заключается в том, что эффективные реализации,
вероятно, достигаются уменьшением отношения обмена памяти к арифметическим операци- ям. Это позволяет полностью использовать векторные операции (если возможно) и позво- ляет встраивать распараллеливание (если возможно). В случае DGEMM число обращений к памяти равно 3n
2
, а число операций с плавающей запятой равно 2n
3
. Следовательно, отно- шение равно 2n/3, что является одним из наивысших в библиотеки BLAS. Например, отно- шение скалярного произведения (уровень 1) равно 1, а отношение для матрично-векторного произведения (уровень 2) равно 2. Следовательно, мы ожидаем получить хорошую произ- водительность для процедуры DGEMM.
Более того, в [
52
] авторы подчёркивают, что многие процедуры BLAS уровня 3 мо- гут быть эффективно реализованы, давая эффективное произведение матрицы на матрицу.
Поэтому хорошая производительность DGEMM является ключевым моментом в общей про- изводительности BLAS. В общем, матрицы A, B и C слишком велики, чтобы уместиться в кэше процессора. Поэтому авторы ATLAS’а [
52
] предлагают использовать алгоритмы с разделением на блоки. Действительно, можно расположить операторы так, чтобы б´
ольшая часть операций была выполнена с данными в кэше, с помощью алгоритмов с разделением на блоки.
Для того, чтобы измерить производительность произведения матрицы на матрицу в
Scilab, мы используем следующий файл-сценарий. Мы используем генератор случайных чи- сел rand для получения матриц, чьи элементы вычисляются из функции нормального рас- пределения. Мы знаем, что генератор случайных чисел в основе rand является плохого статистического качества, и что функция grand может быть использована в качестве заме- ны. Но это не имеет значения на производительность произведения матрицы на матрицу,
так что мы оставим его для простоты. Множитель 1.e6 в вычислении переменной mflops преобразует флопсы в мегафлопсы.
s t a c k s i z e ( " max " );
r a n d ( " n o r m a l " );
n = 1 0 0 0 ;
A = r a n d ( n , n );
B = r a n d ( n , n );
tic ();
C = A * B ;
t = toc ();
m f l o p s = 2* n ^3/ t /1. e6 ;
d i s p ([ n t m f l o p s ])
Программа проверки производительности использовалась в Scilab версии 5.2.2 под
Windows XP 32 бита. ЦП — AMD Athlon 3200+ на 2 ГГц, и система работает с 1 ГБ памяти.
Мы сравниваем здесь производительности трёх библиотек линейной алгебры, предоставлен- ных в Scilab, то есть Reference BLAS, ATLAS и Intel MKL. Для того, чтобы гарантировать выполнимость замеров времени, нам пришлось выбирать отдельно размеры матриц для раз- личных библиотек. Результаты представлены на рисунке
41
В этом эксперименте самыми быстрыми библиотеками являются Intel MKL и ATLAS
136

Reference BLAS
Intel MKL
ATLAS
500 1000 1500 2000 500 1000 1500 2000
Умножение матрицы на матрицу
Порядок матрицы (n)
Мега фло п
сы
Рис. 41: Производительности умножения действительных, квадратных, плотных матриц в
Scilab с различными библиотеками под Windows. Выше — лучше.
(свыше 2000 мегафлопсов), что намного быстрее библиотеки Reference BLAS (около 500
мегафлопсов). Размеры матрицы n, которые мы использовали в эксперименте, зависят от действительной производительности библиотек по причине, которую мы сейчас анализиру- ем. Вы выбрали размер n для которого библиотека показывает время от 0,1 секунды до 8
секунд. Главным преимуществом создания как можно более быстрых программ определения производительности является сохранение хорошей воспроизводимости. Действительно, если n слишком мало, то может случиться, что измеренное время t равно нулю. В противополож- ность этому, библиотека Reference BLAS настолько медленна, что измерение времени растёт так быстро (согласно формуле 2n
3
), что запуск алгоритма для матриц большого размера приводит к очень большим временам, делая процесс замера времени затруднительным на практике.
Теперь сравним действительную производительность библиотек с пиковой производи- тельностью, которую мы можем ожидать для данного конкретного процессора. Процессор —
AMD Athlon 3200+, кодовое имя Venice
2
, для сокета 939. Частота процессора 2 ГГц. Если процессор способен выполнить одну операцию с плавающей запятой за цикл, то частота
2 ГГц подразумевает, что пиковая производительность должна быть более 2 гигафлопсов,
то есть 2000 мегафлопсов. Предыдущий численный эксперимент показывает, что ATLAS и
Intel MKL способны достичь эту производительность и даже чуть больше. Это доказывает,
2
Венеция (прим. перев.)
137

что процессор был способен выполнить от 1 до 2 операций с плавающей запятой на цикл про- цессора. Этот процессор — одноядерный, который поддерживает MMX, 3DNow!, SSE, SSE2
и SSE3. Тот факт, что процессор поддерживает эти наборы инструкций может объяснить почему действительная производительность была слегка выше 2000 мегафлопсов.
Давайте рассмотрим чувствительность производительности к размеру кэша. Может случиться, что библиотека линейной алгебры покажет повышенную производительность для матриц, которые могут быть полностью сохранены в кэше. Цель библиотек ATLAS и Intel
MKL в том, чтобы этого не случилось. Следовательно, производительность не должна за- висеть от от размера кэша. Давайте проверим, что это действительно происходит в этом конкретном эксперименте. Кэш данных L1 равен 64 КБ, а кэш L2 равен 512 КБ. Поскольку число операций с плавающей запятой двойной точности требует 8 байт (т. е. 64 бита), то чис- ло чисел двойной точности, которое может быть сохранено в кэше L1 равно 64000/8 = 8000
чисел двойной точности. Это соответствует плотной, квадратной матрице чисел двойной точ- ности размером 89 × 89. Кэш L2, с другой стороны, может хранить плотную, квадратную матрицу чисел двойной точности размером 256 × 256. Поскольку размеры наших матриц варьируются от 500 × 500 до 2000 × 2000, то мы видим, что общее число чисел двойной точности не может быть (по большей мере) целиком сохранено в кэше. Следовательно, мы приходим к заключению, что алгоритмы, используемые в Intel MKL и ATLAS не слишком чувствительны к размеру кэша.
5.5.2
Обратный слэш
В этом разделе мы измеряем производительность оператора обратный слеш в Scilab’е под Linux.
Этот тест производительности часто называется
«тест производительности
LINPACK», поскольку этот тест был создан при разработке проекта LINPACK. Теста произ- водительности LINPACK используется и сейчас, например, в качестве меры производитель- ности для высокопроизводительных компьютеров, таких, как машины из топа 500 (Top 500)
[
4
], например. Эта библиотека вытеснена BLAS’ом и LAPACK’ом, которые используются в
Scilab’е.
Давайте рассмотрим квадратную, плотную, действительную матрицу A размером n × n и действительный вектор b размером n × 1. Оператор « \» позволяет вычислить решение x линейного уравнения A*x=b. Внутри, если матрица хорошо обусловлена, Scilab производит
LU-разложение матрицы A, используя перемещение строк. Затем, мы выполняем прямую и обратную подстановки, используя LU-разложение. Эти операции основаны на библиотеке
LAPACK, которая, внутри, использует библиотеку BLAS. Точные имена процедур LAPACK,
используемых оператором обратного слэша, представлены в разделе
5.6
Решение систем линейных уравнений является типичным случаем использования биб- лиотек линейной алгебры. Поэтому этот конкретный тест производительности очень часто используется для измерения производительности вычислений с плавающей запятой на кон- кретной машине [
39
,
16
,
42
]. Число операций с плавающей запятой для всего процесса равно
2 3
n
3
+ 2n
2
. Поскольку это число растёт по закону n
3
, а число элементов матрицы лишь n
2
,
то мы можем ожидать хорошей производительности в этом тесте.
Следующий файл-сценарий измерить производительность оператора обратного слэша в Scilab’е. Мы настраиваем генератор случайных чисел в функции rand таким образом, что он формирует матрицы с элементами, вычисленными по функции нормального распреде- ления. Наконец, мы вычисляем число операций с плавающей запятой и делим на 1.e6 для
138

того, чтобы получить мегафлопсы.
n = 1 0 0 0 ;
r a n d ( " n o r m a l " );
A = r a n d ( n , n );
b = r a n d ( n , 1 ) ;
tic ();
x = A \ b ;
t = toc ();
m f l o p s = ( 2 / 3 * n ^3 + 2* n ^ 2 )/ t /1. e6 ;
d i s p ([ n t m f l o p s ])
Мы выполняем этот тест в scilab-5.3.0-beta-4 на Linux Ubuntu 32 бита. Процессор Intel
Pentium M с частотой 2 ГГц, оперативная память 1 ГБ. Размер кэша 2 МБ. Процессор под- держивает наборы инструкций MMX, SSE и SSE2. Рисунок
42
представляет производитель- ности Reference BLAS и ATLAS. Версия ATLAS’а, которую мы используем, предоставлена
Ubuntu. Из-за практических ограничений эти эксперименты остановлены, когда время, тре- буемое на выполнение одного оператора обратного слэша, стало больше 8 секунд.
Reference BLAS
ATLAS
600 700 800 900 1000 1100 1200 1300 1400 1500 0
500 1000 1500 2000 2500 3000 3500
Обратный слеш
Порядок матрицы (n)
Мега фло п
сы
Рис. 42: Производительность решения вещественной, квадратной, плотной, линейной систе- мы уравнений оператором обратного слэша с помощью различных библиотек на Linux.
Как мы можем видеть, библиотека ATLAS улучшает производительность оператора обратного слеша примерно в два раза. Более того, мегафлопсы сильно возрастают с разме- ром матрицы, а производительность библиотеки Reference BLAS кажется остаётся одной и той же. Частота процессора 2 ГГц, это предполагает, что мы можем ожидать пиковую произ- водительность свыше 2000 мегафлопсов. Действительная производительность ATLAS’а на
139

этих матрицах в пределах [1000,1400], что несколько разочаровывает. Поскольку эта про- изводительность возрастает, то пиковая производительность может быть достигнута для б´
ольших матриц, но это здесь не проверялось, поскольку времена возрастают очень быстро,
делая эксперименты долгими для выполнения.
5.5.3
Многоядерные вычисления
В этом разделе мы покажем, что Scilab может выполнять многоядерные вычисления с плавающей запятой, если мы используем, например, библиотеку линейной алгебры Intel
MKL, которая в Windows предоставляется вместе с Scilab’ом.
Мы использовали для этого теста Scilab v5.3.0-beta-4 на операционной системе Windows
Vista Ultimate 32 бита. Процессор Intel Xeon E5410 с 4 ядрами на 2,33 ГГц [
25
,
53
]. Процессор может управлять 4 потоками (то есть один поток на ядро) и имеет L2 кэш на 12 МБ. Опе- рационная система может иметь только 4 ГБ физической памяти. Рисунок
43
показывает результаты проверки производительности на перемножении матрицы на матрицу, которое мы уже представили в разделе
5.5.1
Библиотека
N
Пользовательское время (с)
мегафлопсы
Intel MKL
4766 8,172 26494
ATLAS
1595 3,698 2194
Ref. BLAS
444 0,125 1400
Рис. 43: Производительность произведения матрицы на матрицу на 4-х ядерной системе под
Windows с различными библиотеками.
Учитывая, что этот процессор имеет 4 ядра на 2,33 ГГц, мы ожидаем производитель- ность больше 9,33 гигафлопсов. Действительная производительность Intel MKL равна 26,494
гигафлопса, что указывает на то, что эта библиотека использовала 4 ядра, выполняя 2,8 опе- рации с плавающей запятой за цикл. С другой стороны, библиотеки ATLAS и Ref-BLAS не используют 4 ядра этой машины и дают значительно меньшие результаты.
5.6
Замечания и ссылки
В разделе
5.2.3
мы анализировали производительность алгоритма суммирования. Мож- но было бы указать, что точность суммы зависит от размера данных и конкретных значений в наборе данных. Это проанализировано, например, Уилкинсоном (Wilkinson)[
54
] и Хайэмом
(Higham) [
23
]. Точнее, можно указать, что набор данных может быть плохообусловлен по отношению к сумме, например, если значения смешанных знаков (то есть, некоторые по- ложительные, а некоторые отрицательные) и очень разных амплитуд. В этом случае мы можем наблюдать, что относительная ошибка суммы может быть большой.
В этом документе мы сфокусировались главным образом на производительности и ссы- лаемся на предыдущие ссылки на вопросы точности. На практике сложно полностью разде- лить эти два вопроса. Действительно, гораздо легче получить лучшую производительность,
если мы полностью проигнорируем вопросы точности. Поэтому можно рассматривать опти- мизацию только после того, как мы будем иметь достаточно широкий базис тестов, включая и корректность и точность. Это гарантирует, что во время процесса «оптимизации» мы не введём программную ошибку или, хуже того, ненужную ошибку в результате.
140

Тема, которая тесно связана с векторизацией алгоритмов Scilab’а — это индексирова- ние матриц, то есть, выбор набора элементов в матрице. Статья Эддинса (Eddins) и Шура
(Shure) [
18
] фокусируется на этой теме в контексте Matlab
R
, где авторы представляют индексацию векторов, индексацию матриц с помощью двух подпрограмм, линейную индек- сацию и логическую индексацию на основе логических выражений. Matlab является зареги- стрированной торговой маркой Mathworks.
В разделе
5.2.1
мы представили несколько методов, связанных с принципами компи- ляции. Классической ссылкой по вопросу компиляторов является «Книга дракона
3
» [
5
]. Эта книга вводит в проектирование компиляторов, включая лексический анализ, регулярные выражения, машины конечных состояний и методы проверки синтаксиса.
В разделе
5.4.1
мы кратко представили проект ATLAS. Уоли (Whaley) и Донгарра
(Dongarra) [
51
] описывают метод автоматического формирования и оптимизации числового программного обеспечения для процессоров с глубокой иерархией памяти и конвейерных блоков реализации функций.
В разделе
5.1.3
мы представили алгоритм метода исключения Гаусса и заявили, что этот алгоритм (но не конкретная реализация, которую мы представили) используется в
Scilab’е в основе оператора обратного слеша. Действительное поведение оператора обратно- го слэша в Scilab’е следующее. Сначала мы выполним вызов процедуры DGETRF из LAPACK,
которая раскладывает матрицу A на PA = LU, где P — матрица перестановок, L — ниж- няя треугольная матрица, а U — верхняя треугольная матрица. Диагональные элементы L
являются единицами. Затем мы вызываем процедуру DGECON из LAPACK, которая аппрокси- мирует инвертирование числа обусловленности 1-нормы матрицы A. В этом месте есть два варианта. Если число обусловленности матрицы не слишком велико, то система линейных уравнений не является плохобусловленной (это не точно, поскольку обусловленность только аппроксимирована, но алгоритмы делают гипотезу о том, что оценка может быть достовер- ной). В этом случае мы вызываем процедуру DGETRS из LAPACK, которая решает систему линейных уравнений, используя разложение PA = LU. Точнее, этот шаг использует пере- становку строк, а затем прямой и обратный ход для нахождения решения x в зависимости от правостороннего b. В другом варианте, то есть, если матрица является плохообусловленной,
мы используем модифицированную версию процедуры DGELSY из LAPACK, которая решает соответствующую задачу наименьших квадратов. Этот выбор может уменьшить точность матрицы промежуточного условия [
10
].
В [
39
] Клив Молер (Cleve Moler) анализирует производительность Matlab’а в различ- ных ситуациях, включая решение системы уравнений с помощью оператора обратного сл- эша. Чтобы измерить производительность вычисления, он подсчитывает число флопсов,
измеряя время выполнения в секундах с помощью функций Matlab’а tic и toc и оценивая число операций с плавающей запятой для этого конкретного вычисления. Автор использует матрицу случайных чисел с элементами, полученными из функции нормального распределе- ния и использует функцию randn в Matlab’е. Увеличивая размер матрицы, он наблюдает пик производительности, который соответствует размеру кэша компьютера SPARC-10, который он использует для этого теста. Эта машина имеет кэш в один мегабайт, что соответствует квадратной матрицы порядка
- - > f l o o r ( s q rt (1. e6 / 8 ) )
ans
=
3 5 3 .
3
Название «Книга дракона» происходит от того, что на обложке книги изображены сражающиеся рыцарь и дракон. (прим. перев.)
141

Точнее, Молер использует определение, которое связывает один мегабайт с 2 20
байтами, но результат грубо тот же. Пиковая скорость в мегафлопсах получена для матрицы, которая точно укладывается в кэш. Он заявляет, что этот эффект кэша является следствием библио- теки LINPACK. Он подчёркивает, что это одна из главных мотиваций для проекта LAPACK.
Этот документ был написан в 1994: Matlab использует библиотеку LAPACK с 2000 [
40
]. В
[
16
] Джек Донгарра (Jack Dongarra) собирает различные производительности на этом тесте производительности. Источники по этому вопросу доступны на netlib.org [
42
].
Общие источники, связанные с улучшением производительности Matlab’а, доступны на [
38
].
В разделе
5.4.3
мы обсудили установку Scilab’а на Windows. Эта тема обсуждается подробнее на [
14
], где Аллан Корнэ (Allan Cornet), разработчик из Scilab Consortium, пред- ставляет полный процесс установки Scilab 5.2.0 на Windows.
В разделе
5.4.4
мы представили установку двоичной версии ATLAS для Scilab’а на
Linux/Ubuntu. Поскольку ATLAS оптимизирован на время сборки, то лучше компилировать собственный ATLAS на той машине, которую планируется использовать. Это гарантирует,
что двоичный файл, который мы используем, действительно оптимизирован в соответствии с параметрами текущей машины (например, размер кэша). В [
48
] Сильвестр Ледрю (Sylvestre
Ledru), который является разработчиком из Scilab Consortium и поддерживает пакет Debian
Science, предоставляет информацию об установке ATLAS’а на Debian. Этот вопрос пред- ставлен более глубоко в [
30
].
Есть другие ссылки, которые могут помочь пользователям Scilab’а под Linux. В [
31
]
Ледрю описывает процесс установки Scilab 5.1.1 на Linux. В [
32
] Ледрю обсуждает обнов- ление управления пакетами BLAS и LAPACK в Debian. В [
29
] Ледрю представляет метод,
который позволяет переключаться между различными оптимизированными библиотеками линейной алгебры.
В разделе
5.1.4
мы представили функцию benchfun, которая измеряет производи- тельность указанной функции. В качестве альтернативы мы можем использовать модуль
Scibench [
6
], который предоставляется в ATOMS. Целью проекта Scibench является предо- ставление инструментов для измерения производительности функции. Более того, модуль предоставляет коллекцию программ измерения производительности для измерения произ- водительности Scilab’а. Внутри модуля Scibench функция scibench_benchfun имеет ту же цель, что и benchfun, которая была представлена в разделе
5.1.4 5.7
Упражнения
Упражнение 5.1 (Метод исключения Гаусса с перестановкой строк ) В разде- ле
5.1.3
мы представили две функции gausspivotalnaive и gausspivotal, которые включа- ют в себя алгоритм метода исключения Гаусса с перестановкой строк. Используйте функцию benchfun и сравните действительную производительность этих двух функций.
Упражнение 5.2 (Комбинации ) Рассмотрим две матрицы чисел типа double размером
1 × n x и Y. Допустим, что нам хотелось бы создать функцию, формирующую матрицу c размером 2 × m, содержащую все комбинации векторов x и y, с m = n
2
. То есть, мы хоте- ли бы получить пары [x(1);y(2)], [x(1);y(3)], . . . , [x(n);y(n)]. Например, рассмотрим следующие матрицы x и y.
- - > x = [10 11 1 2 ] ;
- - > y = [13 14 1 5 ] ;
Сформируем следующую матрицу c:
142

- - > c c
=
10.
10.
10.
11.
11.
11.
12.
12.
12.
13.
14.
15.
13.
14.
15.
13.
14.
15.
• Создайте простую реализацию, формирующую матрицу c.
• Используйте произведение Кронекера .*. и создайте векторизованную функцию для получения матрицы c.
• Используйте выделения матриц, чтобы создать векторизованную функцию для фор- мирования матрицы c.
6
Благодарности
Я хочу выразить благодарность Аллану Корнэ (Allan Cornet), Сильвестру Ледрю
(Sylvestre Ledru), Бруно Жофрэ (Bruno Jofret) и Бернарду Хугуэнею (Bernard Hugueney)
за помощь в подготовке этого документа. Я благодарю Сильвестра Ледрю за информацию,
которой он поделился относительно пакета ATLAS на Debian. Бруно Жофрэ ввёл меня в управление памятью в стеке Scilab’а. Я благодарю Аллана Корнэ за информацию, которой он поделился относительно пакетов ATLAS и Intel MKL на Windows. Я благодарю Бенуа
Жепфера (Benoit Goepfert) за его комментарии к этому документу. Я благодарю Сержа
Стира (Serge Steer) и Майка Пэйджеса (Mike Pages) за их обсуждения гиперматриц.
143

7
Ответы к упражнениям
7.1
Ответы к разделу
2
Ответ к упражнению
2.1
(Максимальный размер стека) Результат алгоритма, который устанав- ливает размер стека, зависит от операционной системы. Это результат на Windows-машине:
- - > f o r m a t ( 2 5 )
- - > s t a c k s i z e ( " max " )
- - > s t a c k s i z e ()
ans
=
1 1 0 0 5 4 0 9 6 .
3 4 8 8 1 .
а это результат на Linux-машине
- - > f o r m a t ( 2 5 )
- - > s t a c k s i z e ( " max " )
- - > s t a c k s i z e ()
ans
=
2 8 1 7 6 3 8 4 .
3 5 0 7 7 .
Ответ к упражнению
2.2
(who_user ) Следующий пример — результат на Linux-машине.
- - > w h o _ u s e r ()
Пользовательские переменные:
h o m e
Используется 7 элементов из 4 9 9 0 7 9 8
- - > A = o n e s ( 1 0 0 , 1 0 0 ) ;
- - > w h o _ u s e r ()
Пользовательские переменные:
A
h o m e
Используется 1 0 0 0 9 элементов из 4 9 9 0 7 9 6
- - > h o m e h o m e
=
/ h o m e / mb
На самом деле могут быть иные переменные, определённые пользователем, которые могут быть рас- печатаны функцией who_user. Например, я установил модуль ATOMS «uncprb», и вот что я получаю при запуске:
- - > w h o _ u s e r ()
Пользовательские переменные:
u n c p r b l i b
Используется 217 элементов из 4 9 9 0 7 8 5
Это из-за того, что переменная uncprblib содержит библиотеку, связанную с этим конкретным модулем.
Ответ к упражнению
2.3
(whos ) Следующий пример представляет вывод функции whos, которая распечатывает имя, тип, размер и байты, требуемые всеми переменными в текущей среде.
- - > w h o s
N a m e
Тип
Размер
Байт
$
p o l y n o m i a l
1 на 1 56
% d r i v e r N a m e s t r i n g *
1 на 1 40
%e c o n s t a n t
1 на 1 24
% e p s c o n s t a n t
1 на 1 24
[ . . . ]
g u i l i b l i b r a r y
488
h e l p t o o l s l i b l i b r a r y
728 144

h o m e s t r i n g
1 на 1 56
i n t e g e r l i b l i b r a r y
1 4 1 6
i n t e r p o l a t i o n l i b l i b r a r y
336
[ . . . ]
7.2
Ответы к разделу
3
Ответ к упражнению
3.1
(Поиск файлов) Следующая функция searchSciFilesInDir ищет эти файлы в заданной директории.
f u n c t i o n f i l e m a t r i x = s e a r c h S c i F i l e s I n D i r ( d i r e c t o r y , f u n n a m e )
f i l e m a t r i x = []
l s m a t r i x = ls ( d i r e c t o r y ) ’;
for f = l s m a t r i x i s s c i = r e g e x p ( f , " / ( . * ) . sci / " );
if ( i s s c i < > [] ) t he n s c i f i l e = f u l l f i l e ( d i r e c t o r y , f )
f u n n a m e ( s c i f i l e )
f i l e m a t r i x ( $ +1) = s c i f i l e end end e n d f u n c t i o n
В нашем примере мы просто распечатываем имя файла. Следующая функция mydisplay использует функ- цию mprintf распечатывает имя файла в консоли. Для того, чтобы напечатать короткое сообщение, мы убираем имя директории из отображаемой строки. Для того, чтобы выделить название директории из име- ни файла, мы используем функцию fileparts.
f u n c t i o n m y d i s p l a y ( f i l e n a m e )
[ path , fname , e x t e n s i o n ]= f i l e p a r t s ( f i l e n a m e )
m p r i n t f ( " % s% s \ n " , fname , e x t e n s i o n )
e n d f u n c t i o n
Далее мы используем функцию searchSciFilesInDir для распечатки файлов-сценариев Scilab из дирек- тории Scilab’а, содержащую макросы, связанные с многочленами. Переменная SCI содержит название ди- ректории, в которой находится директория Scilab’а. Для того, чтобы определить абсолютное имя файла,
содержащее имя директории, содержащей макросы, мы используем функцию fullfile, которая собирает свои входные аргументы в единую строку, отделяя директории разделителем, соответствующим текущей операционной системе (слеш «/» под Linux и обратный слеш «\» под Windows).
- - > d i r e c t o r y = f u l l f i l e ( SCI , " m o d u l e s " , " p o l y n o m i a l s " , " m a c r o s " )
d i r e c t o r y
=
/ h o m e / u s e r n a m e / scilab - 5 . 2 . 2 / s h a r e / s c i l a b / m o d u l e s / p o l y n o m i a l s / m a c r o s
- - > f i l e m a t r i x = s e a r c h S c i F i l e s ( d i r e c t o r y , m y d i s p l a y );
i n v r . sci p o l f a c t . sci c m n d r e d . sci
[ . . . ]
h o r n e r . sci r o w c o m p r . sci h t r i a n r . sci
Ответ к упражнению
3.2
(Запрос типизированных списков)
Мы уже видели, что функция definedfields возвращает целые числа с плавающей запятой, связанные с расположением полей в струк- туре данных. Это не совсем то, что нам здесь нужно, но легко построить нашу собственную функцию на этом строительном блоке. Следующая функция isfielddef принимает типизированный список tl и строку
145

fieldname в качестве входных аргументов и выдаёт логическую переменную bool, которая равна истине,
если поле определено, и лжи, если нет. Сначала мы вычисляем ifield, которая является индексом поля fieldname. Чтобы это сделать, мы ищем строку fieldname в полном списке полей, определённом в tl(1) и используем find для поиска требуемого индекса. Затем мы используем функцию definedfields для полу- чения матрицы определённых индексов df. Затем мы ищем индекс ifield в матрице df. Если поиск удался,
то переменная k не будет пустой, что значит, что поле определено. Если поиск не удался, то переменная k будет пустой, то означает, что поле не определено.
f u n c t i o n b o o l = i s f i e l d d e f ( tl , f i e l d n a m e )
i f i e l d = f i n d ( tl ( 1 ) = = f i e l d n a m e )
df = d e f i n e d f i e l d s ( tl )
k = f i n d ( df == i f i e l d )
b o o l = ( k < > [] )
e n d f u n c t i o n
Несмотря на то, что программа довольно абстрактна, она остаётся простой и эффективной. Её главное преимущество состоит в том, чтобы показать как типизированные списки могут вести к очень динамичным структурам данных
7.3
Ответы к разделу
5
Ответ к упражнению
5.1
(Метод исключения Гаусса с перестановкой строк )
Мы используем следующий файл-сценарий, который создаёт матрицу размером 100 × 100 и использует её для сравнения производительности функций gausspivotalnaive и gausspivotal.
s t a c k s i z e ( " max " );
n = 1 0 0 ;
A = g r a n d ( n , n , " def " );
e = o n e s ( n , 1 ) ;
b = A * e ;
b e n c h f u n ( " n a i v e " , g a u s s p i v o t a l n a i v e , l i s t ( A , b ) ,1 ,10);
b e n c h f u n ( " f a s t " , g a u s s p i v o t a l , l i s t ( A , b ) ,1 ,10);
Этот файл-сценарий производит следующий вывод:
n a i v e : 10 i t e r a t i o n s , m e a n = 1 . 5 4 1 3 0 0 , min = 1 . 4 8 2 0 0 0 , max = 2 . 0 3
f a s t : 10 i t e r a t i o n s , m e a n = 0 . 0 2 7 0 0 0 , min = 0 . 0 1 0 0 0 0 , max = 0 . 0 7 0 0 0 0
В этом случае для этого размера матрицы среднее улучшение производительности равно 1,5413/0,027 =
57,08. Фактически, оно могло бы быть больше, при условии, что мы используем достаточно большие матри- цы. Для n = 200 мы наблюдали отношение производительностей больше 100.
Ответ к упражнению
5.2
(Комбинации)
Мы рассматриваем два вектора x и y, с количеством элементов n. Целью этого упражнения является создание функции, способной генерировать матрицу c раз- мером 2 × m, содержащую все комбинации векторов x и y, с количеством элементов m = n
2
Наша первая идея заключается в использовании двух вложенных циклов: эта идея реализована в следующей функции. Первый цикл по ix перебирает элементы в x, в то время, как цикл по iy перебирает элементы в y. Алгоритм устанавливает начальное значение k=1. Каждый раз, когда мы создаём новую комбинацию, мы сохраняем её в k-том столбце матрицы c, а затем мы увеличиваем k.
f u n c t i o n c = c o m b i n a t e v e c t o r s S l o w ( x , y )
cx = s i z e ( x , " c " )
cy = s i z e ( y , " c " )
c = z e r o s (2 , cx * cy )
k = 1
for ix = 1 : cx for iy = 1 : cy c (1:2 , k ) = [ x ( ix ); y ( iy )]
k = k + 1
end
146

end e n d f u n c t i o n
Показанная функция работает, но она медленная, поскольку требуется выполнить n
2
циклов. Следовательно,
целью является удалить эти два вложенных цикла и сформировать комбинации с помощью векторизованных инструкций.
Нашей второй идеей является использование произведения Кронекера для формирования всех ком- бинаций одной-единственной инструкцией. Следующий пример показывает основную идею формирования всех комбинаций векторов [10 11 12] и [13 14 15].
- - >[10 11 12] .*. [1 1 1]
ans
=
10.
10.
10.
11.
11.
11.
12.
12.
12.
- - >[1 1 1] .*. [13 14 15]
ans
=
13.
14.
15.
13.
14.
15.
13.
14.
15.
Следующая реализация является прямым обобщением предыдущего примера.
f u n c t i o n c = c o m b i n a t e v e c t o r s K r o n ( x , y )
cx = s i z e ( x , " c " )
cy = s i z e ( y , " c " )
c (1 ,:) = x .*. o n e s (1 , cy )
c (2 ,:) = o n e s (1 , cx ) .*. y e n d f u n c t i o n
Наша третья идея заключается в использовании операций выделения матрицы. Действительно, пер- вый ряд в c — это просто утроенная копия вектора x. Вот почему мы можем создавать матрицу повторяемых индексов на основе выделения матрицы, и выделять соответствующие элементы за одну-единственную ин- струкцию. Повторенные индексы также созданы с помощью индексов, повторенных выделением матрицы.
Эти идеи представлены в следующем примере, который вычисляет первую строку матрицы c. В этом при- мере мы сначала создаём индексы i1 от 1 до 3. Затем мы создаём индексы j1 многократно выделяя строку.
Затем мы создаём k1 используя функцию matrix, которая преобразует матрицу j1 в вектор-строку индексов.
Наконец, мы используем k1 в качестве индексов столбцов в x.
- - > x = [10 11 1 2 ] ;
- - > i1 = 1:3
i1
=
1.
2.
3.
- - > j1 = i1 ( [ 1 ; 1 ; 1 ] , : )
j1
=
1.
2.
3.
1.
2.
3.
1.
2.
3.
- - > k1 = m a t r i x ( j1 ,1 ,9)
j1
=
1.
1.
1.
2.
2.
2.
3.
3.
3.
- - > x (: , k1 )
ans
=
10.
10.
10.
11.
11.
11.
12.
12.
12.
Мы можем использовать те же принципы, чтобы формировать вторую строку в c. Но эти элементы не просто повторяются, они, к тому же, чередуются. Для того, чтобы это сделать, мы используем ту же самую идею, что и прежде, но транспонируем индексы перед тем, как преобразовать матрицу индексов в вектор индексов.
- - > y = [13 14 1 5 ] ;
- - > i2 = 1:3
i2
=
1.
2.
3.
- - > j2 = i2 ( [ 1 ; 1 ; 1 ] , : )
147

j2
=
1.
2.
3.
1.
2.
3.
1.
2.
3.
- - > k2 = m a t r i x ( j2 ’ ,1 ,9)
k2
=
1.
2.
3.
1.
2.
3.
1.
2.
3.
- - > y (: , k2 )
ans
=
13.
14.
15.
13.
14.
15.
13.
14.
15.
В предыдущем примере важно заметить оператор одиночная скобка в выражении matrix(j2’,1,9). Следу- ющий файл-сценарий является прямым обобщением предыдущих идей.
f u n c t i o n c = c o m b i n a t e v e c t o r s E x t r ( x , y )
// Результаты всех комбинаций двух векторов x и y cx = s i z e ( x , " c " )
cy = s i z e ( y , " c " )
//
i1 = 1: cx j1 = i1 ( o n e s ( cy ,1) ,:)
k1 = m a t r i x ( j1 ,1 , cx * cy )
//
i2 = 1: cy j2 = i2 ( o n e s ( cx ,1) ,:)
k2 = m a t r i x ( j2 ’ ,1 , cx * cy )
//
c = [
x (: , k1 )
y (: , k2 )
]
e n d f u n c t i o n
В следующем файле-сценарии мы используем функцию benchfun для сравнения производительности предыдущих функций.
n = 3 0 0 ;
x = ( 1 : n );
y = ( 1 : n );
b e n c h f u n ( " S l o w " , c o m b i n a t e v e c t o r s S l o w , l i s t ( x , y ) ,1 ,10);
b e n c h f u n ( " K r o n . " , c o m b i n a t e v e c t o r s K r o n , l i s t ( x , y ) ,1 ,10);
b e n c h f u n ( " E x t r . " , c o m b i n a t e v e c t o r s E x t r , l i s t ( x , y ) ,1 ,10);
Этот файл-сценарий обычно производит следующий вывод:
S l o w : 10 i t e r a t i o n s , m e a n = 0 . 2 4 0 1 0 , min = 0 . 2 3 7 0 0 , max = 0 . 2 5 0 0 0
K r o n .: 10 i t e r a t i o n s , m e a n = 0 . 0 1 4 4 0 , min = 0 . 0 1 3 0 0 , max = 0 . 0 2 2 0 0
E x t r .: 10 i t e r a t i o n s , m e a n = 0 . 0 0 9 6 0 , min = 0 . 0 0 7 0 0 , max = 0 . 0 1 5 0 0
Мы видим, что реализация на основе произведения Кронекера больше, чем в 10 раз быстрее в среднем.
Заметим, что нет значительной разницы между функциями combinatevectorsKron и combinatevectorsExtr,
хотя combinatevectorsExtr кажется немного быстрее. Это можно объяснить тем, что combinatevectorsKron выполняет много умножений с плавающей запятой, которые не нужны поскольку один из операндов является единицей.
148

Предметный указатель clear,
16
det,
25
fullfile,
13
getos,
13
horner,
24
list,
31
poly,
23
regexp,
22
roots,
25
stacksize,
6
string,
19
tlist,
33
typeof,
17
type,
17
who,
12
MSDOS,
13
SCIHOME,
13
SCI,
13
TMPDIR,
13
[],
67
add_param,
76
argn,
59
check_param,
81
deff,
54
error,
71
get_function_path,
54
get_param,
76
gettext,
71
init_param,
76
tic,
103
timer,
103
toc,
103
typeof,
54
type,
54
varargin,
59
varargout,
59
warning,
71
ATLAS,
128
BLAS,
128
callback,
57
,
58
cell,
48
DGEMM,
129
flops,
135
funcprot,
56
,
57
Intel,
128
interpreter,
114
LAPACK,
128
macros,
54
MKL,
128
mlist,
44
PCRE,
19
polynomials,
22
primitive,
54
regular expression,
19
struct,
45
ТАД,
38
векторизация,
102
макрос,
54
многочлен,
22
примитив,
54
профиль функции,
105
регулярное выражение,
19
тип абстрактных данных,
38
флопс,
135
функции обратного вызова,
58
число операций с плавающей запятой за се- кунду,
135 149

Список литературы
[1] ATLAS – Automatically Tuned Linear Algebra Software.
http://math-atlas.
sourceforge.net
[2] Blas – Basic Linear Algebra Subprograms.
http://www.netlib.org/blas
[3] Lapack – Linear Algebra PACKage.
http://www.netlib.org/lapack
[4] The Linpack benchmark.
http://www.top500.org/project/linpack
[5] Alfred V. Aho, Ravi Sethi, and Jeffrey D. Ullman. Compilers: principles, techniques, and tools. Addison-Wesley Longman Publishing Co., Inc., Boston, MA, USA, 1986.
[6] Michael Baudin. Scibench.
http://atoms.scilab.org/toolboxes/scibench
[7] Michael Baudin. How to emulate object oriented programming in Scilab.
http://wiki.
scilab.org/Emulate_Object_Oriented_in_Scilab
, 2008.
[8] Michael Baudin. Apifun - Check input arguments in macros.
http://forge.scilab.org/
index.php/p/apifun/
, 2010.
[9] Michael Baudin. The derivative function does not protect against conflicts between user- defined functions and developper-defined variables.
http://bugzilla.scilab.org/show_
bug.cgi?id=7104
, 2010.
[10] Micha¨
el Baudin. The help page of backslash does not print the singularity level - bug report
#7497.
http://bugzilla.scilab.org/show_bug.cgi?id=7497
, July 2010.
[11] Michael Baudin. The integrate function does not protect against conflicts between user- defined functions and developper-defined variables.
http://bugzilla.scilab.org/show_
bug.cgi?id=7103
, 2010.
[12] Michael Baudin. The neldermead component does not protect under particular user-defined objective functions.
http://bugzilla.scilab.org/show_bug.cgi?id=7102
, 2010.
[13] Micha¨
el Baudin. There is no pascal function - bug report #7670.
http://bugzilla.scilab.
org/show_bug.cgi?id=7670
, August 2010.
[14] Allan Cornet. Scilab installation on Windows.
http://wiki.scilab.org/howto/install/
windows
, 2009.
[15] J. J. Dongarra, Jermey Du Cruz, Sven Hammerling, and I. S. Duff. Algorithm 679: A set of level 3 basic linear algebra subprograms: model implementation and test programs. ACM
Trans. Math. Softw., 16(1):18–28, 1990.
[16] Jack J. Dongarra. Performance of various computers using standard linear equations software.
http://www.netlib.org/benchmark/performance.ps
, 2013.
[17] William H. Duquette. Snit’s not incr tcl.
http://www.wjduquette.com/snit/
[18] Steve Eddins and Loren Shure. Matrix indexing in Matlab.
http://www.mathworks.com
150

[19] Jean-Luc Fontaine. Stooop.
http://jfontain.free.fr/stooop.html
[20] Jeffrey Friedl. Mastering Regular Expressions. O’Reilly Media, Inc., 2006.
[21] Gene H. Golub and Charles F. Van Loan. Matrix computations (3rd ed.). Johns Hopkins
University Press, Baltimore, MD, USA, 1996.
[22] Philip Hazel. Pcre – Perl Compatible Regular Expressions.
http://www.pcre.org/
[23] Nicholas J. Higham. Accuracy and Stability of Numerical Algorithms. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, second edition, 2002.
[24] Intel. Intel Math Kernel Library.
http://software.intel.com/en-us/intel-mkl/
[25] Intel. Intel Xeon Processor E5410.
http://ark.intel.com/Product.aspx?id=33080
, 2010.
[26] D. E. Knuth. The Art of Computer Programming, Volume 2, Seminumerical Algorithms.
Third Edition, Addison Wesley, Reading, MA, 1998.
[27] Donald E. Knuth. Fundamental Algorithms, volume 1 of The Art of Computer Programming.
Addison-Wesley, Reading, Massachusetts, second edition, 1973.
[28] Sylvestre Ledru. Localization.
http://wiki.scilab.org/Localization
[29] Sylvestre Ledru. Handle different versions of BLAS and LAPACK.
http://wiki.debian.
org/DebianScience/LinearAlgebraLibraries
, 2010.
[30] Sylvestre Ledru.
Linear algebra libraries in Debian.
http://people.debian.org/
sylvestre/presentation-linear-algebra.pdf
, August 2010.
[31] Sylvestre Ledru. Scilab installation under linux.
http://wiki.scilab.org/howto/install/
linux
, 2010.
[32] Sylvestre Ledru. Update of the linear algebra libraries in Debian.
http://sylvestre.ledru.
info/blog
, April 2010.
[33] Pierre Mar´
echal. Code conventions for the Scilab programming language.
http://wiki.
scilab.org
, 2010.
[34] The Mathworks. Matlab - Product Support - 1106 - Memory Management Guide.
http:
//www.mathworks.com/support/tech-notes/1100/1106.html
[35] The Mathworks. Matlab - Product Support - 1107 - Avoiding Out of Memory Errors.
http:
//www.mathworks.com/support/tech-notes/1100/1107.html
[36] The Mathworks. Matlab - Product Support - 1110 - Maximum Matrix Size by Platform.
http://www.mathworks.com/support/tech-notes/1100/1110.html
[37] The Mathworks.
Matlab - Technical Solution - What are the benefits of 64-bit

Matlab versus 32-bit Matlab?
http://www.mathworks.com/support/solutions/en/data/
1-YVO5H/index.html
[38] The Mathworks. Techniques for improving performance.
http://www.mathworks.com/help/
techdoc/matlab_prog/f8-784135.html
151

[39] Cleve Moler. Benchmarks: Linpack and Matlab - Fame and fortune from megaflops.
http:
//www.mathworks.com/company/newsletters/news_notes/pdf/sumfall94cleve.pdf
,
1994.
[40] Cleve Moler. Matlab incorporates LAPACK.
http://www.mathworks.com
, 2000.
[41] Cleve Moler. The origins of Matlab, December 2004.
http://www.mathworks.fr
[42] netlib.org. Benchmark programs and reports.
http://www.netlib.org/benchmark/
[43] Mathieu Philippe, Djalel Abdemouche, Fabrice Leray, Jean-Baptiste Silvy, and Pierre Lando.
modules/graphics/includes/ObjectStructure.h.
http://gitweb.scilab.org
[44] Bruno Pin¸con. Quelques tests de rapidit´
e entre diff´
erents logiciels matriciels.
http://wiki.
scilab.org/Emulate_Object_Oriented_in_Scilab
, 2008.
[45] pqnelson.
Object oriented
C.
http://pqnelson.blogspot.com/2007/09/
object-oriented-c.html
[46] W. H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical
Recipes in C, Second Edition. Cambridge University Press, 1992.
[47] Axel-Tobias Schreiner. Object-oriented programming with ANSI C.
www.cs.rit.edu/ats/
books/ooc.pdf
[48] Debian
Science.
README.Debian.
http://anonscm.debian.org/viewvc/
debian-science/packages/atlas/trunk/debian
, 2010.
[49] Enrico Segre. Scilab function variables: representation, manipulation.
http://wiki.scilab.
org
, 2007.
[50] Enrico Segre. Scoping of variables in Scilab.
http://wiki.scilab.org/howto/global_and_
local_variables
, 2007.
[51] R. Clint Whaley and Jack J. Dongarra. Automatically tuned linear algebra software. In
Supercomputing ’98: Proceedings of the 1998 ACM/IEEE conference on Supercomputing
(CDROM), pages 1–27, Washington, DC, USA, 1998. IEEE Computer Society.
[52] R. Clint Whaley, Antoine Petitet, R. Clint, Whaley Antoine, Petitet Jack, and Jack J.
Dongarra. Automated empirical optimizations of software and the ATLAS project, 2000.
[53] Wikipedia. Intel Xeon — Wikipedia, the free encyclopedia.
http://en.wikipedia.org/
wiki/Xeon
, 2010.
[54] James Hardy Wilkinson. Rounding errors in algebraic processes. Prentice Hall, 1963.
152