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

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

Вопрос в том как склеить эти два компонента так, чтобы optim могла использовать uncprb?
Очевидно,
что мы не можем передать ни uncprb_getobjfcn ни uncprb_getgrdfcn в качестве входных аргументов в optim, поскольку целевая функция должна вычислить и значение функции и градиент. Решение этой задачи заключается в определении промежуточной функции, которая имеет заголовок, требуемый функцией optim,
и который вызывает функции uncprb_getobjfcn и uncprb_getgrdfcn. Единственной про- блемой является управление nprob, n и m, которое зависит от задачи, которую нужно решить.
Чтобы это сделать, возможны два разных способа.
• Мы можем определить функцию с дополнительными аргументами nprob, n и m и ис- пользовать специальные возможности optim для управления ими.
• Мы можем определить функцию, которая динамически определяется с помощью так,
что в конце переменные nprob, n и m являются константами.
Далее мы изучим оба метода.
Первый метод основан на управлении дополнительными аргументами функции обрат- ного вызова, которые были представлены в разделе
4.6.4
. Определим целевую функцию costfun, которая имеет ожидаемую последовательность вызова, с дополнительными аргу- ментами nprob, n и m.
f u n c t i o n [ f , g , ind ]= c o s t f u n ( x , ind , nprob , n , m )
[ n , m , x0 ] = u n c p r b _ g e t i n i t f ( n p r o b )
// Делаем вектор-столбец x = x (:)
f = u n c p r b _ g e t o b j f c n ( n , m , x , n p r o b )
g = u n c p r b _ g e t g r d f c n ( n , m , x , n p r o b )
e n d f u n c t i o n
Затем, мы можем использовать специальные возможности функции optim, которые управ- ляют случаем, в котором для целевой функции нужны дополнительные входные аргументы.
Аргумент costf может также быть списком (myfun,a1,a2,...). В этом случае myfun, пер- вый элемент списка, должна быть функцией и должна иметь заголовок:
[ f , g , ind ]= m y f u n ( x , ind , a1 , a2 , . . . )
Мы используем эту возможность в следующем файле-сценарии:
n p r o b = 1;
[ n , m , x0 ] = u n c p r b _ g e t i n i t f ( n p r o b );
m y o b j f u n = l i s t ( costfun , nprob , n , m );
[ fopt , x o p t ]= o p t i m ( m y o bj f u n , x0 )
Предыдущий файл-сценарий производит следующий вывод:
- - >[ fopt , x op t ]= o p t i m ( m y o b j f u n , x0 )
x o p t
=
1.
1.
f o p t
=
0.
Мы знаем, что глобальный минимум этой задачи равен x
?
= (1, 1), так что предыдущий результат верен.
Второй метод основан на динамическом создании целевой функции с помощью функ- ции execstr.
100

В следующем примере мы определяем функцию objfun, которая будет передана в ка- честве входного аргумента функции optim. Тело функции objfun динамически определяется использованием переменных nprob, m и n, которые были определены ранее. Заголовок, тело и окончание функции затем конкатенируются в матрицу строковых элементов instr.
h e a d e r = " f u n c t i o n [ fout , gout , ind ]= o b j f u n ( x , ind ) " ;
b o d y = [
" f o u t = u n c p r b _ g e t o b j f c n ( " + s t r i n g ( n )+ " , " + s t r i n g ( m )+ " ,x , " + s t r i n g ( n p r o b )+ " ) "
" g o u t = u n c p r b _ g e t g r d f c n ( " + s t r i n g ( n )+ " , " + s t r i n g ( m )+ " ,x , " + s t r i n g ( n p r o b )+ " ) "
];
f o o t e r = " e n d f u n c t i o n "
i n s t r = [
h e a d e r b o d y f o o t e r
]
Предыдущий исходный код довольно абстрактный. Фактически генерируемая функция про- стая, как та, что показана в следующем примере.
- - > i n s t r i n s t r
=
! f u n c t i o n [ fout , gout , ind ]= o b j f u n ( x , ind )
!
! f o u t = u n c p r b _ g e t o b j f c n (2 ,2 , x ,1)
!
! g o u t = u n c p r b _ g e t g r d f c n (2 ,2 , x ,1)
!
! e n d f u n c t i o n
!
Как мы можем видеть, входной аргумент x функции objfun единственный оставшийся:
остальные параметры заменены их действительными значениями для этой конкретной за- дачи. Следовательно, функция objfun явно определена и не использует своё окружение для получения, например, значения nprob. Мы подчёркиваем этот особый момент, кото- рый показывает, что использование области видимости переменных через стек вызовов, как представлено в разделе
4.5.1
, может быть исключено использованием функции execstr.
Для того, чтобы определить целевую функцию, мы, наконец, используем execstr.
e x e c s t r ( i n s t r )
В следующем примере мы вызываем функцию optim и вычисляем решение первой тестовой задачи.
- - >[ fopt , x op t ]= o p t i m ( objfun , x0 )
x o p t
=
1.
1.
f o p t
=
0.
4.8
Заметки и ссылки
Янн Коллетт (Yann Collette) разработал модуль parameters как инструмент для предо- ставления тех же возможностей, что и у функций optimset/optimget в Matlab. Он заметил,
что функции optimset/optimget нельзя было настроить снаружи: нам приходится модифи- цировать функции для того, чтобы добавить новую опцию. Вот почему модуль parameters был создан таким образом, который позволяет позволить пользователю управлять стольки- ми опциями, сколько потребуется.
101

Область видимости переменных представлена в разделе
4.5.1
. Эта тема также пред- ставлена в [
50
].
Проблемы с функциями обратного вызова были представлены в разделе
4.6
. Эта тема была проанализирована в сообщениях о программных ошибках №7102, №7103 и №7104 [
12
,
11
,
9
].
В [
49
] Энрико Сегре (Enrico Segre) описал несколько возможностей, связанных с функ- циями.
В разделе
4.3
мы представили шаблон для разработки устойчивых функций. Исполь- зуя этот шаблон можно сделать код совместимым с соглашением по коду, представленным
Пьером Марешалем (Pierre Mar´
echal) в [
33
].
В разделе
4.3
мы представляем правила написания устойчивых функций. Написание такой устойчивой функции требует большого числа проверок и может привести к дублиро- ванию кода. Модуль «apifun» [
8
] является экспериментальной попыткой предоставить API
для проверки входных аргументов с б´
ольшей простотой.
5
Производительность
В этом разделе мы представляем ряд уловок программирования в Scilab, которые поз- воляют делать быстрые файлы-сценарии. Эти методы известны как векторизация и явля- ются ядром наиболее эффективных функций Scilab.
Мы покажем, как использовать функции tic и toc для измерения характеристик ал- горитма. Во втором разделе мы анализируем исходный алгоритм и проверяем, что векто- ризация может кардинально улучшить характеристики, обычно на множитель от 10 до 100,
но иногда и более. Затем, мы анализируем методы компилирования и представляем связь между проектированием интерпретатора и характеристиками. Мы представляем возмож- ности профилирования Scilab и даём пример векторизации алгоритма метода исключения
Гаусса. Мы представляем принципы векторизации и сравниваем характеристики циклов с характеристиками векторизованных инструкций. В следующем разделе мы представляем различные методы оптимизации, основанные на практических примерах. Мы представля- ем библиотеки линейной алгебры, используемые в Scilab. Мы представляем числовые биб- лиотеки BLAS, LAPACK, ATLAS и Intel MKL, которые предоставляют оптимизированные функции линейной алгебры и можем увидеть существенную разницу относительно характе- ристик матричных операций. В последнем разделе мы представляем различные измерения производительности Scilab’а, основанные на подсчёте операций с плавающей запятой.
На рисунке
27
представлены некоторые функции Scilab, которые используются в кон- тексте анализа характеристик файлов-сценариев Scilab.
tic, toc измерение пользовательского времени timer измерение системного времени
Рис. 27: Функции измерения производительности.
5.1
Измерение производительности
В этом разделе мы описываем функции, которые позволяют измерить время, требуемое на вычисление. Действительно, перед оптимизацией любой программы нам следует точно
102

измерить её текущую производительность.
В первом разделе мы представляем функции tic, toc и timer. Затем мы сравнива- ем эти функции и подчёркиваем из отличие на многоядерных машинах. Мы представляем функциональности профилирования в Scilab’е, которые позволяю проанализировать части алгоритма, которые наиболее затратные. Наконец, мы представляем функцию benchfun,
которая позволяет провести простой и надёжный анализ производительности.
5.1.1
Основные применения
Для того, чтобы измерить время, требуемое на вычисление, мы можем использовать функции tic и toc, которые измеряют пользовательское время в секундах. Пользователь- ское время — время настенных часов, то есть, время, которое требуется компьютеру от на- чала задачи до её конца. Это включает и само вычисление, конечно, но также все остальные операции системы, такие как обновление экрана, обновление файловой системы, позволе- ние другим процессам делать часть их работы и т. д. В следующем примере мы используем функции tic and toc для измерения времени вычисления собственных значений случайной матрицы на основе функции spec.
- - > tic (); l a m b d a = s p e c ( r a n d ( 2 0 0 , 2 0 0 ) ) ; t = toc ()
t
=
0.1
Вызов функции tic запускает счётчик, а вызов функции toc — останавливает его, воз- вращая число прошедших секунд. В этом случае мы печатаем инструкции в одной той же строке, используя разделитель ;: действительно, если бы мы напечатали инструкции интер- активно на нескольких строчках, то измеренное время было бы главным образом временем,
требуемым на набор инструкций с клавиатуры.
Предыдущее измерение времени не очень надёжно в том смысле, что если бы мы вы- полнили одни и те же инструкции несколько раз, то не получили бы одного и того же времени. Следующий пример демонстрирует это поведение.
- - > tic (); l a m b d a = s p e c ( r a n d ( 2 0 0 , 2 0 0 ) ) ; t = toc ()
t
=
0 . 1 8
- - > tic (); l a m b d a = s p e c ( r a n d ( 2 0 0 , 2 0 0 ) ) ; t = toc ()
t
=
0 . 1 9
Возможное решение этой проблемы заключается в том, чтобы выполнить одни и то же вычисление несколько раз, скажем, 10 раз, например, и затем распечатать минимальное,
максимальное и среднее время.
for i = 1 : 10
tic ();
l a m b d a = s p e c ( r a n d ( 2 0 0 , 2 0 0 ) ) ;
t ( i ) = toc ();
end
[ min ( t ) m e a n ( t ) max ( t )]
Предыдущий пример даёт следующий вывод.
- - >[ min ( t ) m e a n ( t ) max ( t )]
ans
=
0 . 1 4 1 0 . 2 2 1 4 0 . 3 3 103

Если другая программ на компьютере использует ЦП в то же время, когда Scilab вы- полняет свои вычисления, то пользовательское время может возрасти. Следовательно функ- ции tic и toc не используются в ситуациях, когда имеет значение только время ЦП. В
этом случае мы можем вместо этого использовать функцию timer, которая измеряет вре- мя системы. Оно соответствует времени, требуемому ЦП для выполнения специфических вычислений. требуемых Scilab’ом, а не другими процессами.
Следующий пример представляет функцию timer.
for i = 1 : 10
t i m e r ();
l a m b d a = s p e c ( r a n d ( 2 0 0 , 2 0 0 ) ) ;
t ( i ) = t i m e r ();
end
[ min ( t ) m e a n ( t ) max ( t )]
The previous script produces the following output.
- - >[ min ( t ) m e a n ( t ) max ( t )]
ans
=
0 . 1 0 0 1 4 4 0 . 1 1 6 1 6 7 0 0 . 1 6 0 2 3 0 4 5.1.2
Пользовательское время и время ЦП
В этом разделе мы анализируем разницу между пользовательским временем и време- нем ЦП на одно- и многоядерных системах.
Давайте рассмотрим следующий файл-сценарий, который позволяет выполнить про- изведение двух квадратных матриц. Мы измеряем пользовательское время функциями tic и toc, а время ЦП измеряем функцией timer.
s t a c k s i z e ( " max " )
n = 1 0 0 0 ;
A = r a n d ( n , n );
B = r a n d ( n , n );
t i m e r ();
tic ();
C = A * B ;
t U s e r = toc ();
t C p u = t i m e r ();
d i s p ([ t U s e r t C p u t C p u / t U s e r ])
Следующий пример представляет результат, когда мы запустили этот файл-сценарий на Linux-машине с одним ядром на 2 ГГц.
- - > d i s p ([ t U s e r t C p u t C p u / t U s e r ])
1 . 4 6 9 1 . 3 4 8 0 8 4 0 . 9 1 7 6 8 8 2
Как мы можем видеть два времени довольно близки и показывают, что 91 % времени на- стенных часов потрачено центральным процессором на конкретное вычисление.
Теперь рассмотрим следующий пример, выполненный на Windows-машине с 4-мя яд- рами, где Scilab использует многопотоковую Intel MKL.
- - > d i s p ([ t U s e r t C p u t C p u / t U s e r ])
0 . 2 4 5 1 . 0 2 9 6 0 6 6 4 . 2 0 2 4 7 5 9
Мы видим что есть существенная разница между временем ЦП и временем настенных часов.
Отношение близко к 4, что равно числу ядер. Действительно, на многоядерной машине,
если Scilab использует многопотоковую библиотеку, функция timer сообщает сумму времён,
104

потраченных на каждое ядро. Вот почему во многих ситуациях нам следует использовать функции tic и toc для измерения производительности алгоритма.
5.1.3
Профилирование функции
В этом разделе мы представляем возможности профилирования в Scilab’е. Профиль функции раскрывает части функции, которые наиболее затратны. Это позволяет нам сосре- доточиться на частях исходного кода, которые наиболее всего требуется улучшить.
Сначала мы рассмотрим исходную, медленную реализацию алгоритма метода исклю- чения Гаусса и проанализируем её профиль. Затем мы используем векторизацию для того,
чтобы ускорить вычисления и проанализируем этот обновлённый профиль.
Функции в таблице
28
позволяют управлять профилированием функций.
add_profiling
Добавляет команды профилирования в функцию plotprofile
Выделяет и отображает профили исполнения profile
Выделение профилей исполнения remove_profiling
Удаляет инструкции профилирования reset_profiling
Сбрасывает счётчики профилирования showprofile
Выделяет и отображает профили исполнения
Рис. 28: Команды Scilab’а для профилирования функций.
В этом разделе мы рассматриваем решение систем линейных уравнений Ax = b, где
A — вещественная матрица размером n × n, а b — вектор-столбец размером n × 1. Одним из наиболее стабильных алгоритмов для этой задачи является метод исключения перемен- ных Гаусса с перестановкой строк. Метод исключения переменных Гаусса с перестановкой строк используется в Scilab’е с помощью оператора обратный слеш \. Для того, чтобы про- иллюстрировать профилирование функции, мы не будем использовать оператор обратного слеша. Вместо этого мы разработаем наш собственный алгоритм метода исключения Гаусса
(который, очевидно, на практике не является хорошей идеей). Давайте вспомним, что пер- вый шаг алгоритма требует разложить A в виде PA = LU, где P — матрица перестановок,
L — нижняя треугольная матрица, а U — верхняя треугольная матрица. Сделав однажды,
алгоритмы продолжают выполнение прямого и, затем, обратного исключения.
Алгоритм и короток и сложен, и фактически довольно интересный. Мы не ставим целью получение этого алгоритма через презентацию в этом документе (см. Голуб и Ван
Лоан (Golub и Van Loan) [
21
] для более полного анализа). Этот алгоритм является хоро- шим кандидатом для профилирования, поскольку он более сложен, чем средняя функция.
Следовательно, от него не получить напрямую хорошую производительность. Этот тот тип ситуации, где профилирование наиболее полезно. Действительно, интуиции иногда не до- статочно, чтобы точно обнаружить «виновные строки», то есть строки, которые являются узким местом для производительности. В общем и целом, измерение производительности часто лучше, чем гадание.
Следующая функция gausspivotalnaive — это простая реализация метода исключе- ния Гаусса с перестановкой строк. Она принимает матрицу A и b в качестве аргумента с правой стороны, и возвращает решение x. Она объединяет разложение PA = LU с пря- мым и обратным ходом в одной единственной функции. Множители матрицы L сохранены в нижней части матрицы A, а множители матрицы U сохранены в верхней части матрицы
A.
105

1
f u n c t i o n x = g a u s s p i v o t a l n a i v e ( A , b )
2
n = s i z e ( A , " r " )
3
// Инициализируем перестановку
4
P = z e r o s (1 , n -1)
5
// Выполняем преобразования Гаусса
6
for k =1: n -1 7
// Поиск ведущего элемента
8
mu = k
9
a b s p i v o t = abs ( A ( mu , k ))
10
for i = k : n
11
if ( abs ( A ( i , k )) > a b s p i v o t ) t h en
12
mu = i
13
a b s p i v o t = abs ( A ( i , k ))
14
end
15
end
16
// Перестановка строк k и mu из столбцов от k до n
17
for j = k : n
18
tmp = A ( k , j )
19
A ( k , j ) = A ( mu , j )
20
A ( mu , j ) = tmp
21
end
22
P ( k ) = mu
23
// Выполнение преобразования для строк от k +1 до n
24
for i = k +1: n
25
A ( i , k ) = A ( i , k )/ A ( k , k )
26
end
27
for i = k +1: n
28
for j = k +1: n
29
A ( i , j ) = A ( i , j ) - A ( i , k ) * A ( k , j )
30
end
31
end
32
end
33
// Выполнение прямого хода
34
for k =1: n -1 35
// Перестановка b ( k ) и b ( P ( k ))
36
tmp = b ( k )
37
mu = P ( k )
38
b ( k ) = b ( mu )
39
b ( mu ) = tmp
40
// Подстановка
41
for i = k +1: n
42
b ( i ) = b ( i ) - b ( k ) * A ( i , k )
43
end
44
end
45
// Выполнение обратного хода
46
for k = n : -1:1 47
t = 0.
48
for j = k +1: n
49
t = t + A ( k , j ) * b ( j )
50
end
51
b ( k ) = ( b ( k ) - t ) / A ( k , k )
52
end
53
x = b
54
e n d f u n c t i o n
Перед тем, как действительно измерять производительность, мы должны проверить корректность. поскольку нет смысла делать быстрее функцию, в которой имеются про-
106

граммные ошибки. Это кажется очевидным, но на практике частот не берётся во внимание.
В следующем примере мы формируем матрицу A размером 10 × 10 со значениями элемен- тов случайно распределёнными на интервале [0, 1). Мы выбрали этот тестовый случай за его очень малое число обусловленности. Мы выбираем ожидаемое решение e с единичными элементами и вычисляем правую сторону b из уравнения b=A*e. Наконец, мы используем нашу функцию gausspivotalnaive для вычисления решения x.
n = 10;
A = g r a n d ( n , n , " def " );
e = o n e s ( n , 1 ) ;
b = A * e ;
x = g a u s s p i v o t a l n a i v e ( A , b );
Затем мы проверяем, что число обусловленности матрицы не слишком большое. В сле- дующем примере мы вычисляем число обусловленности матрицы A по 2-норме, и проверяем,
что оно приблизительно равно 10 1
, что довольно мал´
о (только порядок величины играет зна- чение). Мы также проверим, что относительная ошибка для x имеет ту же амплитуду, что и %eps≈ 10
−16
, что превосходно.
- - > l o g 1 0 ( c o nd ( A ))
ans
=
1 . 3 8 9 8 4 1 7
- - > n o r m ( e - x )/ n o r m ( e )
ans
=
6 . 4 5 5 D -16
Требовалась бы более полная проверка, но мы сейчас допускаем, что мы можем верить нашей функции.
В следующем примере мы вызовем функцию add_profiling для того, чтобы профи- лировать функцию gausspivotalnaive.
- - > a d d _ p r o f i l i n g ( " g a u s s p i v o t a l n a i v e " )
Внимание : переопределение функции : g a u s s p i v o t a l n a i v e
Выполните f u n c p r o t (0) для отключения этого сообщения
Целью функции add_profiling является переопределение функции, которую профилируют,
так чтобы Scilab мог сохранить число выполнений каждой строки. Как указано в предупре- ждении, она в действительности изменяет профилируемую функцию, делая её медленнее.
Это не играет значения, поскольку мы не принимаем слишком серьёзно времена выполнения на графике профиля: мы, вместо этого, будем фокусироваться только на число раз, которое конкретная строчка исходного кода была выполнена.
Для того, чтобы измерить производительность нашей функции, мы должны выполнить её, как в следующем примере:
x = g a u s s p i v o t a l n a i v e ( A , b );
Во время прогона, система профилирования сохранила данные профилирования, так что мы можем теперь анализировать. В следующем примере мы вызываем функцию plotprofile для того, чтобы построить профиль функции.
p l o t p r o f i l e ( g a u s s p i v o t a l n a i v e )
Предыдущий пример формирует фигуру
29
. График сделан из трёх гистограмм. Каж- дая функция имеет 54 строчки, которые представлены на оси X каждой гистограммы. Ось Y
зависит от гистограммы. На первой гистограмме ось Y представляет число вызовов (colls)
каждой соответствующей строки в профилируемой функции. Это измеряется просто об- новлением счётчика каждый раз, когда выполняется строчка. На второй гистограмме ось
107

Y представляет сложность каждой соответствующей строчки в профилируемой функции.
Это измеряет усилия интерпретатора для выполнения соответствующей строчки. Ось Y на третьей гистограмме представляет время ЦП (CPU time), в секундах, требуемое на выпол- нение соответствующей строчки. На практике, измерение времени ЦП ненадёжно, когда оно недостаточно велико.
0 50 100 150 200 250 10 20 30 40 50
Количество вызовов
0 500 1000 1500 2000 2500 3000 3500 4000 4500 10 20 30 40 50
Сложность
0.000 0.002 0.004 0.006 0.008 0.010 10 20 30 40 50
Время процессора (общее 0,01 сек)
Рис. 29: Профиль функции gausspivotalnaive.
Функция plotprofile также открывает редактор и правит профилируемую функцию.
Мышкой мы можем щёлкнуть левой кнопкой на цветные вертикальные полосы гистограммы.
В редакторе это немедленно переместит курсор на соответствующую строчку. Мы сначала щёлкаем левой кнопкой мыши на самую большую полосу, которая соответствует строчке,
которая выполнялась более 250 раз. Редактор затем перемещает курсор на строчки №23
и №24 в нашей функции. Это позволяет немедленно обратиться к следующим строчкам исходного кода.
for i = k +1: n for j = k +1: n
A ( i , j ) = A ( i , j ) - A ( i , k ) * A ( k , j )
В оставшейся части этого раздела мы собираемся постоянно использовать профили- ровщик. Анализируя различные части функции, которые наиболее затратные, мы можем измерить часть, которая имеет наибольшее влияние на производительность. Это позволя- ет применить принципы векторизации только на строчки исходного кода, которые имеют значение, и позволит нам сохранить большое количество времени. Как только мы это сде-
108

лаем, мы щёлкаем на меню «Exit» («Выход») в графическом окне, что закрывает график профиля.
Теперь проанализируем предыдущие вложенные циклы по переменным i и j. Мы мо- жем сделать циклы по j быстрее, используя векторизованные инструкции. Используя опе- ратор двоеточие :, мы заменяем переменную j на инструкцию k+1:n, как в следующем примере.
for i = k +1: n
A ( i , k +1: n ) = A ( i , k +1: n ) - A ( i , k ) * A ( k , k +1: n )
Мы можем распространить принцип далее, и векторизовать также цикл по переменной i,
опять же используя оператор двоеточие. Это приводит к следующему исходному коду.
A ( k +1: n , k +1: n ) = A ( k +1: n , k +1: n ) - A ( k +1: n , k ) * A ( k , k +1: n )
Эффект предыдущей инструкции заключается в
обновлении подматрицы
A(k+1:n,k+1:n) только одной инструкцией. Ключевым моментом является то, что преды- дущая инструкция обновляет множество элементов матрицы только за один вызов Scilab’а.
В данной ситуации, число обновлённых элементов равно (n − k)
2
, что может быть много,
если n велико. В общем, мы заменили примерно n
2
вызовов обработки интерпретатором матрицы 1 × 1 одним вызовом обработки матрицы размером n × n. Заметим, что умноже- ние A(k+1:n,k)*A(k,k+1:n) — это умножение вектора-столбца на вектор-строку, что даёт матрицу размером (n − k) × (n − k) с помощью только одной инструкции.
Аналогично, мы можем векторизовать другие инструкции в алгоритме.
Цикл в следующем исходном коде является очевидным кандидатом на векторизацию.
for i = k +1: n
A ( i , k ) = A ( i , k )/ A ( k , k )
Предыдущий алгоритм математически эквивалентен, но гораздо медленнее, чем следующая векторизованная инструкция.
A ( k +1: n , k ) = A ( k +1: n , k )/ A ( k , k )
Поиск ведущих элементов в следующей части алгоритма кажется менее очевидным кандидатом на векторизацию.
mu = k a b s p i v o t = abs ( A ( mu , k ))
for i = k : n if ( abs ( A ( i , k )) > a b s p i v o t ) t h en mu = i a b s p i v o t = abs ( A ( i , k ))
end end
По существу, этот алгоритм является поиском элемента наибольшей величины в од- ном подстолбце матрицы A. Следовательно, мы можем использовать функцию max, которая возвращает как максимальный элемент её аргумента так и индекс, который получил это максимальное значение. Остаётся комбинировать несколько векторизованных функций для достижения нашей цели, что нетривиально.
Анализируя, как обновляется индекс i во время цикла, мы видим, что рассматривае- мая подматрица — A(k:n,k). Функция abs работает поэлементно, то есть, выполняет своё
вычисление элемент за элементом. Следовательно, инструкция abs(A(k:n,k)) вычисляет абсолютные величины, в которых мы заинтересованы. Мы можем, наконец, вызвать функ- цию max для выполнения поиска, как в следующем примере.
109

[ a b s p i v o t , m u r e l ]= max ( abs ( A ( k : n , k )))
mu = m u r e l + k - 1
Заметьте, что функция max выполняет только поиск в части столбца матрицы A. Следо- вательно переменная murel хранит строку, достигнувшую максимума, относительно блока
A(k:n,k). Вот почему мы используем инструкцию mu = murel + k - 1, которая позволяет передать глобальный индекс строки в mu.
Следующий цикл позволяет переставлять строки k и mu матрицы A.
for j = k : n tmp = A ( k , j )
A ( k , j ) = A ( mu , j )
A ( mu , j ) = tmp end
Сейчас для нас очевидно, что можно заменить предыдущий набор команд на следующий,
векторизованный.
tmp = A ( k , k : n )
A ( k , k : n ) = A ( mu , k : n )
A ( mu , k : n ) = tmp
Это уже огромное улучшение предыдущего кода. Заметьте, что он создаст промежуточный вектор tmp, что не обязательно, если мы используем следующую инструкцию.
A ([ k mu ] , k : n ) = A ([ mu k ] , k : n )
Действительно, матрица [mu k] является корректной матрицей индексов и может непосред- ственно использоваться для определения части матрицы A.
В общем следующая функция gausspivotal собирает все предыдущие улучшения и использует только векторизованные инструкции.
f u n c t i o n x = g a u s s p i v o t a l ( A , b )
n = s i z e ( A , " r " )
// Инициализация перестановок
P = z e r o s (1 , n -1)
// Выполнение преобразований Гаусса for k =1: n -1
// Поиск ведущего элемента
[ a b s p i v o t , m u r e l ]= max ( abs ( A ( k : n , k )))
// Сдвиг mu , поскольку max возвращает относительно ( k : n , k )
mu = m u r e l + k - 1
// Перестановка строк k и mu из столбцов от k до n
A ([ k mu ] , k : n ) = A ([ mu k ] , k : n )
P ( k ) = mu
// Выполнение преобразования для строк от k +1 до n
A ( k +1: n , k ) = A ( k +1: n , k )/ A ( k , k )
A ( k +1: n , k +1: n ) = A ( k +1: n , k +1: n ) - A ( k +1: n , k ) * A ( k , k +1: n )
end
// Выполнение прямого прохода for k =1: n -1
// Перестановка b ( k ) и b ( P ( k ))
mu = P ( k )
b ([ k mu ]) = b ([ mu k ])
// Замена b ( k +1: n ) = b ( k +1: n ) - b ( k ) * A ( k +1: n , k )
end
// Выполнение обратного прохода
110

for k = n : -1:1
b ( k ) = ( b ( k ) - A ( k , k +1: n ) * b ( k +1: n )) / A ( k , k )
end x = b e n d f u n c t i o n
Ещё раз, мы можем проверить что наша функция работает, как показано в следующем примере.
- - > x = g a u s s p i v o t a l ( A , b );
- - > n o r m ( e - x )/ n o r m ( e )
ans
=
6 . 3 3 9 D -16
Мы можем обновить профиль, запустив функцию plotprofile с улучшенной функцией gausspivotal. Она формирует фигуру
30 0
2 4
6 8
10 5
10 15 20 25 30
Количество вызовов
0 50 100 150 200 250 300 350 5
10 15 20 25 30
Сложность
0e+00 1e-07 2e-07 3e-07 4e-07 5e-07 6e-07 7e-07 8e-07 9e-07 5
10 15 20 25 30
Время процессора (общее 0 сек)
Рис. 30: Профиль функции gausspivotal.
Мы видим, что профиль числа вызовов гладкий, что означает, что нет строк, которые вызываются чаще других. Поскольку это великолепный результат, то мы можем остановить здесь наш процесс векторизации.
В упражнении
5.1
, мы сравниваем действительные производительности этих двух функ- ций. Мы видим, что для матрицы размером 100 × 100 (что скромно), среднее улучшение производительности может быть более 50. В общем, использование векторизации позволяет улучшить производительность на порядки.
111

5.1.4
Функция benchfun
В этом разделе мы представляем функцию, которая может использоваться для анализа производительности алгоритмов. Мы представляем изменчивость времени ЦП, требуемого для вычисления верхней треугольной матрицы Паскаля и представляем более надёжный способ измерения производительности алгоритма.
Попытаемся измерить производительность следующей функции pascalup_col, кото- рая вычисляет верхнюю матрицу Паскаля.
f u n c t i o n c = p a s c a l u p _ c o l ( n )
c = eye ( n , n )
c (1 ,:) = o n e s (1 , n )
for i = 2:( n -1)
c (2: i , i +1) = c ( 1 : ( i -1) , i )+ c (2: i , i )
end e n d f u n c t i o n
Алгоритм, представленный здесь, основан на постолбцовом доступе к матрице, который особенно быстр, поскольку это именно тот способ, которым Scilab хранит элементы матрицы значений типа double. Эта тема представлена более глубоко в разделе
5.3.6
В следующем примере мы измеряем время, требуемое для вычисления верхней тре- угольной матрицы Паскаля размером 1000 × 1000.
- - > tic (); p a s c a l u p _ c o l ( 1 0 0 0 ) ; toc ()
ans
=
0 . 1 7 2 0 1 1
- - > tic (); p a s c a l u p _ c o l ( 1 0 0 0 ) ; toc ()
ans
=
0 . 1 5 6 0 1
Как мы можем видеть, есть некоторая случайность в измерении затраченного времени. Эта случайность может быть объяснена несколькими фактами, но какие бы причины ни были,
нам нужен более надёжный способ измерения времени, требуемого алгоритму. Общий ме- тод заключается в запуске алгоритма несколько раз и усреднении измерений. Этот метод используется в следующем примере. где мы запускаем функцию pascalup_col 10 раз. За- тем мы распечатываем минимальное, среднее и максимальное времена, требуемые во время эксперимента.
- - > for k = 1 : 10
- - >
tic ();
- - >
p a s c a l u p _ c o l ( 1 0 0 0 ) ;
- - >
t ( k ) = toc ();
- - > end
- - > d i s p ([ min ( t ) m e a n ( t ) max ( t )])
0 . 1 0 8 0 0 6 0 . 1 1 6 8 0 7 2 0 . 1 6 8 0 1
Предыдущий метод может быть автоматизирован, что ведёт к функции benchfun, ко- торая представлена ниже. Эта функция запускает алгоритм для которого мы хотим по- лучить надёжные измерения производительности. Входными аргументами являются имя функции name, которую тестируют, функция fun, которую запускают, список iargs, содер- жащий входные аргументы fun, число выходных аргументов nlhs и число итераций для выполнения kmax. Тело benchfun основано на выполнении функции fun, которая выполня- ется kmax раз. На каждой итерации измеряется производительность fun на основе функций tic() и toc(), которые измеряют затраченное функцией время (настенные часы). Функция benchfun возвращает матрицу t данных типа double размером kmax1 × 1, которая содержит
112

различные замеры времени, записанные во время исполнения fun. Она также возвращает строковую переменную msg, которая содержит дружественное к пользователю сообщение,
описывающее производительность тестируемого алгоритма.
f u n c t i o n [ t , msg ] = b e n c h f u n ( n a m e , fun , i a r g s , n l h s , k m a x )
// Вычисляем строку инструкции, которую надо запустить ni = l e n g t h ( i a r g s )
i n s t r = " "
// Кладём аргументы левой стороны i n s t r = i n s t r + " [ "
for i = 1 : n l h s if ( i > 1 ) t h e n i n s t r = i n s t r + " , "
end i n s t r = i n s t r + " x " + s t r i n g ( i )
end i n s t r = i n s t r + " ] "
// Кладём аргументы правой стороны i n s t r = i n s t r + " = fun ( "
for i = 1 : ni if ( i > 1 ) t h e n i n s t r = i n s t r + " , "
end i n s t r = i n s t r + " i a r g s ( " + s t r i n g ( i )+ " ) "
end i n s t r = i n s t r + " ) "
// Цикл по тестам for k = 1 : k m a x tic ()
i e r r = e x e c s t r ( i n s t r , " e r r c a t c h " )
t ( k ) = toc ()
end m s g f m t = " %s : %d i t e r a t i o n s , m e a n = %f , min = %f , max = %f \ n "
msg = m s p r i n t f ( msgfmt , name , kmax , m e a n ( t ) , min ( t ) , max ( t ))
m p r i n t f ( " %s \ n " , msg )
e n d f u n c t i o n
В следующем примере мы прогоняем функцию pascalup_col 10 раз для того, чтобы получить верхнюю треугольную матрицу Паскаля размером 2000 × 2000.
- - > b e n c h f u n ( " p a s c a l u p _ c o l " , p a s c a l u p _ c o l , l i s t ( 2 0 0 0 ) , 1 , 10 );
p a s c a l u p _ c o l : 10 i t e r a t i o n s , m e a n = 0 . 4 3 1 6 2 , min = 0 . 4 1 0 5 9 , max = 0 . 4 7 0 6 7
Хотя функцию benchfun можно было бы сделать более устойчивой и более гибкой,
этого достаточно для многих случаев. Мы будем часто использовать её в этом разделе.
5.2
Основы векторизации
В этом разделе мы представляем основы векторизации. Этот метод программирования является лучшим способом достичь хорошей производительности в Scilab. Первый раздел представляет основы интерпретатора, а также связь между примитивом и шлюзом. Мы сравниваем производительность циклов for с векторизованными инструкциями. В конце мы представляем пример анализа производительности, где мы преобразуем простой, медленный алгоритм в быстрый, векторизованный алгоритм.
113

5.2.1
Интерпретатор
В этом разделе мы кратко представляем методы компилирования так, что мы смо- жем иметь чёткие представления о значении интерпретатора на вопросы производительно- сти. Мы также представляем связь между интерпретатором и шлюзами, которые являются функциями, которые связывают Scilab с нижележащими библиотеками.
Для того, чтобы понять что именно случается, когда выполняется такая инструкция как y=f(x), мы должны взглянуть на структуру компилятора. Действительно, это позволит нам увидеть связь между инструкциями на уровне интерпретатора и выполнением сносно скомпилированного исходного кода на уровне библиотеки.
Компилятор [
5
] — это программа, которая читает программу на одном языке и перево- дит её в эквивалентную программу. Если целевая программа является исполнимой, написан- ной на машинном языке, она может быть вызвана пользователем для обработки входных значений и получения выходных значений. Обычные компиляторы этого типа являются компиляторами Fortran и C/C++ (среди других).
В отличие от предыдущей схемы, Scilab не является компилируемым языком: от интер- претируемый язык. Вместо производства программы на машинном языке, Scilab использу- ет как программу в исходном коде, так и входные данные, и непосредственно производит выходные данные. Рисунок
31
представляет интерпретатор. Распространённые интерпрета- торы этого типа — Scilab, Matlab и Python (среди прочих).
Исходный код программы
Выходные данные интерпретатор
Входные данные
Рис. 31: Scilab — это интерпретатор.
Программа на машинном языке, сформированная компилятором вообще быстрее, чем программа интерпретатора в преобразовании входных данных в выходные. Интерпретатор,
однако, может обычно предоставить лучшую диагностику, чем компилятор, поскольку он исполняет исходный код программы инструкцию за инструкцией. Более того, Scilab исполь- зует высокооптимизированные числовые библиотеки, которые в некоторых ситуациях могут быть быстрее, чем простая реализация на основе компилируемого языка. Эта тема будет рассмотрена в разделе
5.4.1
, где мы представляем числовые библиотеки линейной алгебры,
используемые в Scilab.
Компилятор обычно структурирован как последовательность операций. Первая опе- рация включает в себя чтение символов исходного кода и выполняет лексический анализ.
Этот шаг состоит в группировании символов в смысловые последовательности, именуемые лексемами. Лексемы затем передаются на анализатор синтаксиса (или парсер), который использует лексемы и создаёт древовидное представление исходного кода, который ассоции- руется с грамматической структурой исходного кода. В семантическом анализе компилятор использует дерево синтаксиса и проверяет код на соответствие определению языка. Этот шаг включает в себя проверку типов операндов каждого оператора.
114

Затем могут быть включены несколько промежуточных операций, ведущих к финаль- ному вызову библиотечных процедур. В Scilab этот финальный шаг реализован в шлюзах.
Действительно, каждая функция Scilab’а или оператор неодинаково связан с вычислитель- ными процедурами, внедрёнными в компилированный язык, обычно на C или Fortran. Шлюз читает входные аргументы, проверяет их типы и соответствие, выполняет необходимые вы- числения и затем выталкивает обратно интерпретатору выходные значения. Такая структу- ра представлена на рисунке
32
Шлюз
Интерпретатор
Файл-сценарий
Библиотека
Рис. 32: Scilab связывает интерпретатор с коллекцией библиотек через шлюзы.
5.2.2
Циклы в сравнении с векторизацией
В этом разделе мы представляем простой пример векторизации и анализируем почему циклы for являются хорошими кандидатами на векторизацию.
Рассмотрим следующую матрицу значений типа double, состоящую из n = 5 элементов.
x =[1 2 3 4 5];
Предположим, что мы хотим вычислить сумму 5-ти чисел, хранимых в переменной x.
Мы сравним два файла-сценария, которые дают один и тот же верный результат, но имеют очень разную производительность.
В первом файле-сценарии мы вызываем функцию sum и устанавливаем переменную y.
y = sum ( x );
Это требует только один вызов интерпретатора и включает в себя только один вызов к одному отдельному шлюзу. Внутри шлюза численная библиотека выполняет циклы по n = 5
в компилированном и оптимизированном исходном коде.
Во втором файле-сценарии мы используем циклы for.
y = 0;
n = s i z e ( x , " * " );
for i = 1 : n y = y + x ( i );
end
Предыдущий файл-сценарий привлекает интерпретатор по крайней мере 2 + n раз, где n —
размер матрицы x. При n = 5 это 7 вызовов интерпретатора, но линейно возрастает с ростом размера x. Внутри шлюза численная библиотека выполняет одно единственное суммирова- ние каждый раз, когда она вызывается. Когда число циклов большое, то дополнительные затраты на использование инструкции for вместо векторизованной инструкции являются слишком большими.
Поэтому векторизованные инструкции гораздо быстрее циклов. В общем, нам следует делать минимальное количество вызовов интерпретатора и позволить ему работать с доста- точно большими наборами данных.
Теперь мы подчеркнём заключение предыдущего обсуждения.
115

• Если Scilab имеет встроенную функцию, то нам следует рассмотреть её использова- ние до того, как создавать свою собственную. С точки зрения производительности в большинстве случаев они хорошо спроектированы, то есть, они как только возмож- но используют векторизованные инструкции. Более того, в большинстве случаев они принимают во внимание проблемы чисел с плавающей запятой, так что мы можем получить и быстрое вычисление и хорошую точность.
• Если мы должны создать свой собственный алгоритм, нам следует использовать век- торизованные инструкции где только возможно. Действительно, ускорение векторизо- ванных инструкции обычно 10–100, но может быть иногда 1000. В оставшейся части этого документа мы представим несколько примеров улучшения производительности.
Более того, поскольку Scilab может вызывать оптимизированные библиотеки, которые эффективно используют процессор, когда он должен выполнять вычисление с плавающей запятой. Файлы-сценарии Scilab могут быть такими же эффективными как компилирован- ный исходный код, но они могут быть даже быстрее, чем нехитрый компилированный исход- ный код. Например, библиотеки линейной алгебры, используемые в Scilab’е, берут в расчёт размер кэша процессора и, более того, может выполнять многопотоковые операции. Это центральный момент в Scilab’е, который является преимущественно матричным языком. В
разделе
5.4
представлены оптимизированные библиотеки линейной алгебры, используемые в Scilab’е. Не все программы могут получить выгоду от операций линейной алгебры, вот по- чему мы фокусируемся на методах векторизации, которые можно применять в большинстве ситуаций. Действительно, наш опыт говорит, что векторизация требует некоторой практи- ки. Следующий раздел представляет практический пример анализа производительности на основе векторизации.
5.2.3
Пример анализа производительности
В этом разделе мы сравниваем производительность простого и векторизованного ал- горитма. Мы проверяем, что ускорение в данном конкретном случае может быть равно 500.
Предположим, что мы имеем вектор-столбец, и что мы хотим вычислить среднее зна- чение этого вектора. У нас две возможности выполнить эту задачу:
• использовать функцию mean, имеющуюся в Scilab,
• создать собственный алгоритм.
Функция mean может работать с входными матрицами различных размеров (векторы- строки, векторы-столбцы или общие матрицы) и с различными нормами (1-норма, 2-норма,
∞-норма и норма Фробениуса). Вот почему, в общем, нам не следует переопределять нашу собственную функцию, если в Scilab’е она уже реализована. Более того, производительность функции mean превосходна, как мы увидим. Но в целях иллюстрации вопросов производи- тельности мы разработаем наш собственный алгоритм.
Следующая функция mynaivemean вычисляет среднее значение входного вектора v.
Она выполняет цикл по данным с прямым накоплением суммы.
f u n c t i o n y = m y n a i v e m e a n ( v )
y = 0
n = l e n g t h ( v )
for i = 1: n y = y + v ( i )
116

end y = y / n e n d f u n c t i o n
Перед оптимизацией нашего алгоритма, мы проверяем, что он выполняет свою задачу пра- вильно: действительно, было бы потерей времени делать быстрее алгоритм, содержащий ошибку. . . Алгоритм, кажется, правильный, по крайней мере, с точки зрения следующего примера.
- - > y = m y n a i v e m e a n ( 1 : 9 )
y
=
5.
Для того, чтобы проверить наш алгоритм, мы сформируем большие входные векторы v. Для этого, мы используем функцию rand, поскольку мы знаем, что производительность нашего алгоритма не не зависит от действительных значений входного вектора. Мы подчёр- киваем эту деталь, поскольку практические алгоритмы могут показывать более или менее быструю производительность в зависимости от своих входных параметров. Сейчас не тот случай, потому что наш анализ упрощённый.
Мы формируем векторы размером n = 10
j с j = 1, 2, . . . , 6. Чтобы сформировать целые числа n, мы используем функцию logspace, которая возвращает значения на основе функ- ции логарифма по основанию 10. Следующий файл-сценарий выполняет циклы от n = 10 1
до n = 10 6
и измеряет время, требуемое нашим простым алгоритмом.
n _ d a t a = l o g s p a c e (1 ,6 ,6);
for i = 1:6
n = n _ d a t a ( i );
v = r a n d ( n , 1 ) ;
tic ();
y m e a n = m y n a i v e m e a n ( v );
t = toc ();
m p r i n t f ( " i = %d , n = %d , m y n a i v e m e a n = %f , t = %f \ n " , i , n , y m e a n , t );
end
Далее представлена распечатка производимая файлом-сценарием.
i =1 , n =10 , m y n a i v e m e a n = 0 . 5 0 3 6 9 4 , t = 0 . 0 0 0 0 0 0
i =2 , n =100 , m y n a i v e m e a n = 0 . 4 7 6 1 6 3 , t = 0 . 0 0 0 0 0 0
i =3 , n =1000 , m y n a i v e m e a n = 0 . 5 0 5 5 0 3 , t = 0 . 0 1 0 0 1 4
i =4 , n =10000 , m y n a i v e m e a n = 0 . 5 0 0 8 0 7 , t = 0 . 0 9 0 1 3 0
i =5 , n = 1 0 0 0 0 0 , m y n a i v e m e a n = 0 . 4 9 9 5 9 0 , t = 0 . 5 7 0 8 2 1
i =6 , n = 1 0 0 0 0 0 0 , m y n a i v e m e a n = 0 . 5 0 0 2 1 6 , t = 5 . 2 5 7 5 6 0
Мы видим, что время растёт очень быстро и является довольно большим для n = 10 6
На самом деле, можно разработать гораздо более быстрый алгоритм. В предыдущем файле-сценарии мы заменим функцию mynaivemean на встроенную функцию mean. Это даёт следующую распечатку.
i =1 , n =10 , m e a n = 0 . 4 9 3 5 8 4 , t = 0 . 0 0 0 0 0 0
i =2 , n =100 , m e a n = 0 . 5 1 7 7 8 3 , t = 0 . 0 0 0 0 0 0
i =3 , n =1000 , m e a n = 0 . 4 8 4 5 0 7 , t = 0 . 0 0 0 0 0 0
i =4 , n =10000 , m e a n = 0 . 5 0 1 1 0 6 , t = 0 . 0 0 0 0 0 0
i =5 , n = 1 0 0 0 0 0 , m e an = 0 . 5 0 2 0 3 8 , t = 0 . 0 0 0 0 0 0
i =6 , n = 1 0 0 0 0 0 0 , m ea n = 0 . 4 9 9 6 9 8 , t = 0 . 0 1 0 0 1 4
Как мы видим, функция mean гораздо быстрее нашей функции mynaivemean. В данном случае мы видим, что для n=1000000 отношение по скорости приблизительно 5,25/0,01 =
525.
117

Для того, чтобы улучшить наш простой алгоритм, мы можем использовать функ- цию sum, которая возвращает сумму элементов её входного аргумента. Следующая функция myfastmean выполняет вычисление на основе функции sum.
f u n c t i o n y = m y f a s t m e a n ( v )
y = sum ( v );
n = l e n g t h ( v );
y = y / n ;
e n d f u n c t i o n
Следующая распечатка представляет производительность нашей улучшенной функции.
i =1 , n =10 , m e a n = 0 . 4 1 9 3 4 4 , t = 0 . 0 0 0 0 0 0
i =2 , n =100 , m e a n = 0 . 5 0 8 3 4 5 , t = 0 . 0 0 0 0 0 0
i =3 , n =1000 , m e a n = 0 . 5 0 1 5 5 1 , t = 0 . 0 0 0 0 0 0
i =4 , n =10000 , m e a n = 0 . 5 0 1 6 0 1 , t = 0 . 0 0 0 0 0 0
i =5 , n = 1 0 0 0 0 0 , m e an = 0 . 4 9 9 1 8 8 , t = 0 . 0 0 0 0 0 0
i =6 , n = 1 0 0 0 0 0 0 , m ea n = 0 . 4 9 9 9 9 2 , t = 0 . 0 1 0 0 1 4
Производительность нашей быстрой функции теперь приблизительно та же, что и произво- дительность встроенной функции mean. В действительности, функция mean является мак- росом, который выполняет преимущественно тот же самый алгоритм, что и наша функция,
то есть он использует функцию sum.
Почему же функция sum гораздо быстрее алгоритма на основе цикла for? Причина в том, что функция sum является примитивом, основанным на компилированном исходном коде. Это легко проверить, получив тип этой функции, как в следующем примере.
- - > t y p e o f ( sum )
ans
=
f p t r
Следовательно, цикл, выполняемый Scilab’ом внутри функции sum компилированным ис- ходным кодом. Это действительно быстрее, чем цикл for, выполненный интерпретатором,
и объясняет разницу производительностей между нашей и двумя предыдущими реализаци- ями.
5.3
Уловки оптимизации
В этом разделе мы представляем общие уловки для получения более быстрых алгорит- мов и обратим внимание на использование векторизованных инструкций. Мы представляем опасное использование матриц, чей размер растёт динамически во время работы алгоритма.
Мы представляем функции общего назначения, такие, как векторизованные функции сор- тировки и поиска, которые могут быть объединены для создания эффективных алгоритмов.
Мы также обсуждаем ориентацию матриц и даём пример улучшения производительности на основе этой уловки.
5.3.1
Опасность динамических матриц
В этом разделе мы анализируем вопросы производительности, которые вызваны ис- пользованием матриц, чей размер растёт динамически. Мы показываем числовой экспери- мент, где предопределённая матрица известного размера гораздо быстрее матрицы, которой позволено расти динамически.
118

Давайте рассмотрим ситуацию, когда мы должны вычислить матрицу, но не знаем заранее конечный размер этой матрицы. Это случается когда, например, матрица прочита- на из файла данных или интерактивно вводится пользователем программы. В этом случае мы можем использовать то, что размер матрицы может расти динамически. Но это тре- бует затрат, которые могут быть чувствительными, когда число динамических обновлений велико. Эти затраты связаны с внутренним поведением обновления размера матрицы. Дей- ствительно, увеличение размера матрицы подразумевает, что интерпретатору приходится обновлять свою внутреннюю структуру данных. Следовательно, всё содержимое матрицы должно быть скопировано, что может потребовать значительное количество времени. Бо- лее того, пока циклы обрабатываются, размер матрицы на самом деле растёт, что в свою очередь подразумевает и более медленное копирование матрицы.
Для того, чтобы проанализировать вопросы производительности, мы сравним три функ- ции с разными реализациями, но с одинаковым выходом. Целью является получение мат- рицы A=[1,2,3,...,n], где n — большое целое число. Чисто в целях эффективности мы могли бы использовать оператор двоеточия «:» и инструкцию A=(1:n). Для того, чтобы проиллюстрировать конкретный вопрос, который является темой этого раздела, мы не бу- дем использовать оператор двоеточия, а, наоборот, мы используем следующие реализации.
Функция growingmat1 принимает n в качестве входного аргумента и выполняет ди- намическое вычисление матрицы A. Сначала алгоритм инициализирует матрицу A в виде пустой матрицы. Затем он использует синтаксис
$+1 для динамического увеличения раз- меров матрицы и сохранения нового элемента i.
f u n c t i o n g r o w i n g m a t 1 ( n )
A = []
for i = 1 : n
A ( $ +1) = i end e n d f u n c t i o n
Мы не возвращаем матрицу A в качестве выходного аргумента, поскольку это необязательно для нашей демонстрации.
Следующая функция growingmat2 также выполняет динамическое вычисление мат- рицы A, но использует другой синтаксис. В этот раз матрица динамически обновляется с помощью синтаксиса [A,i], который создаёт новую матрицу-строку с элементом i, добав- ленным к концу.
f u n c t i o n g r o w i n g m a t 2 ( n )
A = []
for i = 1 : n
A = [ A , i ]
end e n d f u n c t i o n
Это создаёт матрицу-строку. Вариант этого алгоритма создавал бы матрицу-столбец син- таксисом [A;i]. Этот вариант не изменит фундаментально поведение алгоритма.
Следующая функция growingmat3 полностью отличается от предыдущих. Вместо ини- циализации матрицы пустой матрицей, она предварительно создаёт матрицу-столбец нулей.
Затем элементы заполняются один за другим.
f u n c t i o n g r o w i n g m a t 3 ( n )
A = z e r o s ( n ,1)
for i = 1 : n
A ( i ) = i
119

end e n d f u n c t i o n
В следующем файле-сценарии мы сравниваем производительности трёх предыдущих функций с n=20000 и повторяем замер времени 10 раз для получения более надёжной оцен- ки производительности. Мы используем функцию benchfun, которая была представлена в разделе
5.1.4
s t a c k s i z e ( " max " )
b e n c h f u n ( " g r o w i n g m a t 1 " , g r o w i n g m a t 1 , l i s t ( 2 0 0 0 0 ) , 0 , 10 );
b e n c h f u n ( " g r o w i n g m a t 2 " , g r o w i n g m a t 2 , l i s t ( 2 0 0 0 0 ) , 0 , 10 );
b e n c h f u n ( " g r o w i n g m a t 3 " , g r o w i n g m a t 3 , l i s t ( 2 0 0 0 0 ) , 0 , 10 );
В следующем примере мы получаем замер времени трёх алгоритмов.
- - > e x e c b e n c h _ d y n a m i c m a t r i x . sce ;
g r o w i n g m a t 1 : 10 i t e r a t i o n s , m e a n = 1 . 5 6 8 2 5 5
g r o w i n g m a t 2 : 10 i t e r a t i o n s , m e a n = 0 . 9 0 1 2 9 6
g r o w i n g m a t 3 : 10 i t e r a t i o n s , m e a n = 0 . 0 2 3 0 3 3
Мы видим, что предварительное создание матрицы A с её известным размером n гораздо быстрее. Отношение от самого медленного до самого быстрого более 50.
5.3.2
Линейная индексация матрицы
В этом разделе мы покажем как можно получить доступ к обычной матрице с по- мощью одного-единственного линейного индекса k вместо пары подындексов (i,j). Эта возможность работает как для матриц, так и для гиперматриц, хотя мы в этом разделе рассматриваем только обычных матрицы.
11