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

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

(1)
12 (2)
13 (3)
14 (4)
15 (5)
16 (6)
17 (7)
18 (8)
19 (9)
20 (10)
21 (11)
22 (12)
Рис. 33: Линейное индексирование матрицы размером 3 × 4. Например, значение «21» рас- положено в ячейке с подындексами (i, j) = (2, 4) и имеет линейный индекс k = 11. Линейные индексы ведут себя так, как если бы элементы были сохранены столбец за столбцом.
Если матрица A имеет m строк, то линейный индекс k, соответствующий подындексам
(i, j), равен k = i + (j − 1)m. В этом случае инструкции A(i,j) и A(k) эквивалентны. Мы можем понимать эту возможность как следствие того, что матрицы хранятся столбец за столбцом (смотри подробности по этой теме в разделе
5.3.6
). В следующем примере мы со- здаём матрицу значений типа double размером 3×4. Затем мы обнуляем элемент с линейным индексом k=11. Этот элемент соответствует элементу с подындексами (i, j) = (2, 4).
- - > A = m a t r i x ( 1 0 + ( 1 : 1 2 ) , 3 , 4 )
A
=
11.
14.
17.
20.
12.
15.
18.
21.
13.
16.
19.
22.
- - > A ( 1 1 ) = 0
A
=
120

11.
14.
17.
20.
12.
15.
18.
0.
13.
16.
19.
22.
Предыдущий пример представлен на рисунке
33
Функции sub2ind и ind2sub делают преобразование из подындексов в линейные ин- дексы и обратно. Эти функции представлены на рисунке
34
. Подчеркнём, что эти функции работают как с матрицами, так и гиперматрицами, но мы остановимся в оставшейся части этого раздела на обычных матрицах.
ind2sub преобразует линейные индексы k в подындексы i1,i2,...
sub2ind преобразует подындексы i1,i2,... линейные индексы k
Рис. 34: Функции Scilab’а для преобразования из подындексов в линейные индексы.
Например, мы можем преобразовать линейный индекс в пару подындексов с помощью инструкции [i,j]=ind2sub(dims,k), где dims является матрицей 1×2, содержащей размеры матрицы. Далее представлен пример функции ind2sub. Рассмотрим матрицу 3 × 4 (как на рисунке
33
) и проверим, что линейный индекс k=11 соответствует подындексам i=2 и j=4.
- - >[ i , j ]= i n d 2 s u b ([3 ,4] ,11)
j
=
4.
i
=
2.
Теперь дадим простой пример как можно использовать функцию sub2ind на практике,
когда используется векторизованное вычисление. Допустим, что у нас есть матрица 3 × 4 и мы хотим обнулить первую диагональ над главной диагональю. В нашем конкретном случае эта диагональ соответствует значениям 14, 18 и 22, которые расположены в строках 1, 2 и
3 и столбцах 2, 3 и 4. Конечно, мы могли бы использовать циклы for, но тогда это не будет векторизованным вычислением. Мы так же не можем напрямую использовать подындексы этих строк и столбцов, поскольку это обнулило бы всю соответствующую подматрицу 3 × 3,
как показано в следующем примере.
- - > A = m a t r i x ( 1 0 + ( 1 : 1 2 ) , 3 , 4 )
A
=
11.
14.
17.
20.
12.
15.
18.
21.
13.
16.
19.
22.
- - > A ([1 2 3] ,[2 3 4]) = 0
A
=
11.
0.
0.
0.
12.
0.
0.
0.
13.
0.
0.
0.
Следовательно, мы используем функцию sub2ind для преобразования этих подындексов в вектор линейных индексов k и затем обнулим эти элементы.
- - > A = m a t r i x ( 1 0 + ( 1 : 1 2 ) , 3 , 4 )
A
=
11.
14.
17.
20.
12.
15.
18.
21.
13.
16.
19.
22.
- - > k = s u b 2 i n d ([3 ,4] ,[1 2 3] ,[2 3 4])
121

k
=
4.
8.
12.
- - > A ( k ) = 0
A
=
11.
0.
17.
20.
12.
15.
0.
21.
13.
16.
19.
0.
Предыдущий пример может использоваться в качестве основы для алгоритмов, которые обращаются к диагоналям матриц, таким, как, например, в функции spdiags.
В следующем разделе
5.3.5
мы пересмотрим пример линейного индексирования. Там мы будем использовать функцию find для обращения к отдельным элементам матрицы.
5.3.3
Обращение через матрицу логических значений
В этом разделе мы покажем как можно обратиться к обычной матрице с помощью индексов, определённых матрицей логических значений.
Операторы, которые создают матрицы логических значений, такие, как == или >, мо- гут использоваться для обращения к матрицам с помощью векторизованных инструкций.
Предположим, что A является матрицей m × n значений типа double, и B является матрицей m × n логических значений. Следовательно, инструкция A(B) обращается ко всем элементам
A, для которых B равно истине. Это показано в следующем примере, где мы вычисляем все элементы A, которые равны 1, и обнуляем их.
- - > A = [1 2 1 2; 2 1 2 1; 1 3 1 3]
A
=
1.
2.
1.
2.
2.
1.
2.
1.
1.
3.
1.
3.
- - > B = ( A = = 1 )
B
=
T F T F
F T F T
T F T F
- - > A ( B ) = 0
A
=
0.
2.
0.
2.
2.
0.
2.
0.
0.
3.
0.
3.
Предыдущая последовательность операций может быть упрощена до одной-единствен- ной инструкции A(A==1)=0. Это показано в следующем примере.
- - > A = [1 2 1 2; 2 1 2 1; 1 3 1 3]
A
=
1.
2.
1.
2.
2.
1.
2.
1.
1.
3.
1.
3.
- - > A ( A = = 1 ) = 0
A
=
0.
2.
0.
2.
2.
0.
2.
0.
0.
3.
0.
3.
В ситуации, когда поиск отдельных элементов A является слишком сложным, чтобы управ- ляться таким методом, мы можем использовать функцию find, как показано в разделе
5.3.5 122

5.3.4
Повтор строк или столбцов вектора
В этом разделе мы покажем как создать копии строк или столбцов вектора с помощью векторизованных инструкций.
Интересным свойством матриц является то, что выделение матрицы может быть ис- пользовано для повтора элементов одной матрицы. Действительно, если A является матри- цей, то мы можем многократно выделять первую строку, что создаёт копии этой строки.
- - > A = 1:3
A
=
1.
2.
3.
- - > B = A ([1 1 1] ,:)
B
=
1.
2.
3.
1.
2.
3.
1.
2.
3.
Интересным здесь является то, что создание B не требует циклов: оно векторизовано.
Эта уловка используется в функции repmat, которая может быть интегрирована в
Scilab 5.4.
Эта тема изучена в упражнении
5.2 5.3.5
Комбинирование векторизованных функций
Для того, чтобы получить хорошую производительность в Scilab’е, мы можем исполь- зовать векторизованные инструкции. В этом разделе мы представляем как комбинировать векторизованные функции.
Векторизованные файлы-сценарии лучше используют оптимизированные библиотеки
Scilab’а. Главное правило — избегать циклы for.
Вторым правилом является использовать матричные операции как можно чаще, по- скольку они используют высокооптимизированные библиотеки линейной алгебры, которые предоставлены Scilab’ом. Однако, есть ситуации, где выполнение векторизованных вычис- лений не очевидно. В этих случаях мы можем комбинировать векторизованные функции:
целью этого раздела в том, чтобы представить относительно продвинутые векторизованные функции и подсказать как их комбинировать.
Рисунок
35
представляет наиболее общие векторизованные функции, которые мы мо- жем использовать для получения хорошей производительности.
+,-,*,/, \
алгебраические операторы
.*.
произведение Кронекера min,max минимум, максимум sum сумма элементов матрицы cumsum накапливаемая сумма элементов матрицы prod произведение элементов матрицы cumprod накапливаемое произведение элементов матрицы find поиск истинных индексов в булевой матрице gsort сортирует элементы матрицы matrix генерирует матрицу из элементов матрицы
Рис. 35: Функции Scilab наиболее часто используемые при векторизации.
123

Алгебраические операторы +, -, *, /, \ представлены здесь из-за того, что они предо- ставляют превосходную производительность в Scilab’е. Для того, чтобы использовать эти операторы, мы должны сформировать матрицы к которым мы можем применить линейную алгебру. Мы можем использовать zeros, ones и другие матрично-ориентированные функ- ции, но полезный оператор для формирования матриц — это произведение Кронекера. В
следующем примере мы перемножаем две матрицы размерами 3 × 2 произведением Кроне- кера, которое создаёт матрицу размером 4 × 9.
- - >[1 2 3;4 5 6] .*. [7 8 9;1 2 3]
ans
=
7.
8.
9.
14.
16.
18.
21.
24.
27.
1.
2.
3.
2.
4.
6.
3.
6.
9.
28.
32.
36.
35.
40.
45.
42.
48.
54.
4.
8.
12.
5.
10.
15.
6.
12.
18.
Произведение Кронекера может быть, например, использовано в комбинаторных алгорит- мах, где мы хотим получить комбинации величин из данного набора. Этот вопрос исследован в упражнении
5.2
Рассмотрим теперь функции find, gsort, min и max, которые имеют специальные вход- ные и выходные аргументы, которые делают их очень полезными, когда мы делаем поиск в быстрых программах.
Использование функции find является общим методом оптимизации, который может заменить алгоритмы поиска. Если B является матрицей m × n логических значений, то инструкция k=find(B) возвращает линейные индексы k, для которых B равна истине.
Если поиск завершился неудачно, то есть, если ни один из элементов не соответствует условию, то функция find возвращает пустую матрицу. Можно подумать, что необходимо явно управлять этим конкретным случаем. Однако инструкция A([])=0 ничего не дела- ет, поскольку пустая матрица [] соответствует пустому набору индексов. Следовательно на выходе функции find в большинстве случаев нет нужды обрабатывать случай неудачи отдельной инструкцией if.
В следующем примере мы обнуляем все элементы матрицы A,которые равны 1.
- - > A = [1 2 1 2; 2 1 2 1; 1 3 1 3]
A
=
1.
2.
1.
2.
2.
1.
2.
1.
1.
3.
1.
3.
- - > B = ( A = = 1 )
B
=
T F T F
F T F T
T F T F
- - > k = fi n d ( B )
k
=
1.
3.
5.
7.
9.
11.
- - > A ( k ) = 0
A
=
0.
2.
0.
2.
2.
0.
2.
0.
0.
3.
0.
3.
Этот пример можно укоротить объединив функции в одну-единственную инструкцию, как показано ниже.
- - > A = [1 2 1 2; 2 1 2 1; 1 3 1 3]
124

A
=
1.
2.
1.
2.
2.
1.
2.
1.
1.
3.
1.
3.
- - > A ( f i n d ( A = = 1 ) ) = 0
A
=
0.
2.
0.
2.
2.
0.
2.
0.
0.
3.
0.
3.
Поэкспериментируем что будет, если значение, которое мы искали, не найдено.
- - > k = fi n d ( A = = 9 9 9 )
k
=
[]
- - > A ( k )=0
A
=
0.
2.
0.
2.
2.
0.
2.
0.
0.
3.
0.
3.
Конечно, как показано в разделе
5.3.3
, мы могли бы использовать более короткую инструкцию A(A==1)=0.
Однако, есть случаи, когда линейные индексы kдолжны быть обработаны перед тем,
как использовать их в установке значений матрицы. Более того, есть у функции опции find,
которые можно использовать для создания более сложных векторизованных операций, как мы увидим далее.
Если используются не все элементы, удовлетворяющие условию, а лишь их часть, то мы можем использовать второй входной аргумент функции find. Например, инструкция k=find(A(:,2)>0,1) возвращает первый индекс второго столбца A, который положитель- ный.
Функции min и max могут быть использованы с двумя выходными аргументами. Дей- ствительно, инструкция [m,k]=min(A) возвращает минимальное значение m и индекс k эле- мента в A, который имеет минимальное значение. Другими словами, индекс k такой, что
A(k)==m. Например, инструкция [m,k]=min(abs(A)) сочетает функцию min и поэлементную функцию abs. Это позволяет получить индекс k того элемента матрицы A, который имеет минимальное абсолютное значение.
Функция gsort имеет также второй аргумент, который может быть полезен для получе- ния хорошей производительности. Действительно, инструкция [B,k]=gsort(A) возвращает сортированную матрицу B и линейные индексы k элементов в исходной матрице. Други- ми словами, матрица B такова, что B(:)==A(k). На практике мы можем применять ту же перестановку во второй матрице C с помощью инструкции C = C(k).
5.3.6
Постолбцовое обращение быстрее
Мы можем получить б´
ольшую производительность с помощью постолбцового обраще- ния к матрице, нежели с помощью построчному. В этом разделе мы представляем пример улучшения производительности с помощью изменения ориентации алгоритма и представим объяснение этому поведению.
Вопрос, который мы собираемся рассмотреть, обсуждался в сообщении об ошибке
№7670 [
13
]. Проблема заключается в вычислении матрицы Паскаля, которая может быть получена из формулы бинома. Ряд отдельных наборов с j элементами, которые могут быть
125

выбраны из набора A с n элементами, является биномиальным коэффициентом и обознача- ется
n j

. Он определяется как
n j

=
n.(n − 1) . . . (n − j + 1)
1.2 . . . j
(1)
Для целых чисел n > 0 и 0 < j < n, биномиальные коэффициенты удовлетворяют
n j

=
n − 1
j

+
n − 1
j − 1

(2)
Следующие два файла-сценария позволяют вычислить матрицу Паскаля. В первой ре- ализации мы вычисляем нижнюю треугольную матрицу Паскаля и выполняем построчный алгоритм.
// P a s c a l l o w e r t r i a n g u l a r : Row by row v e r s i o n f u n c t i o n c = p a s c a l l o w _ r o w ( n )
c = eye ( n , n )
c (: ,1) = o n e s ( n ,1)
for i = 2:( n -1)
c ( i +1 ,2: i ) = c ( i , 1: ( i - 1 ) ) + c ( i ,2: i )
end e n d f u n c t i o n
Предыдущая реализация была предложена мне Каликстом Денизетом (Calixte Denizet) в обсуждении, связанном с сообщением о программной ошибке, но похожая идея была пред- ставлена Самуэлем Гужоном (Samuel Gougeon) в том же потоке. Во второй реализации мы вычисляем верхнюю треугольную матрицу Паскаля и выполняем постолбцовый алгоритм.
// P a s c a l u p p e r t r i a n g u l a r : C o l u m n by c o l u m n v e r s i o n 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
Следующий пример показывает, что обе реализации математически корректны.
- - > p a s c a l u p _ c o l (5)
ans
=
1.
1.
1.
1.
1.
0.
1.
2.
3.
4.
0.
0.
1.
3.
6.
0.
0.
0.
1.
4.
0.
0.
0.
0.
1.
- - > p a s c a l l o w _ r o w ( 5 )
ans
=
1.
0.
0.
0.
0.
1.
1.
0.
0.
0.
1.
2.
1.
0.
0.
1.
3.
3.
1.
0.
1.
4.
6.
4.
1.
Но между двумя этими алгоритмами есть существенная разница в производительности. Ри- сунок
36
представляет производительность этих двух алгоритмов для различных размеров
126

матриц. Рисунок представляет отношение между постолбцовой и построчной реализациями.
Как видим, постоблцовая реализация в целом на 10 % быстрее построчной реализации. За
0.85 0.90 0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30 500 1000 1500 2000 2500 3000 3500
Матрица Паскаля
Порядок матрицы (n)
Врем я постол б
цовое /
В
ремя п о
стро ч
н о
е
Рис. 36: Сравнение постолбцового и построчного алгоритмов для матрицы Паскаля.
исключением редких экспериментов, этот факт верен для матриц всех размеров.
Причиной этого является внутренняя структура данных для матрицы чисел типа double,
которая ориентирована на Фортран. Поэтому элементы матрицы сохраняются столбец за столбцом как показано на рисунке
37
. Следовательно, если мы упорядочим операторы так, что они выполняли постолбцовый алгоритм, то мы будем иметь быстрое обращение к элементам, из-за расположения данных. Это позволяет использовать наилучшим образом процедуру DCOPY (из библиотеки BLAS), которая используется для создания промежу- точных векторов, используемых в алгоритме. Действительно, выражение c(2:i,i) создаёт промежуточный вектор-столбец, а выражение c(i,2:i) создаёт промежуточный вектор- строку. Это требует создания временных векторов, которые должны заполняться данными,
скопированными из матрицы c. Когда вектор-столбец c(2:i,i) создан, то все элементы находятся друг за другом. Вектор-строка c(i,2:i), напротив, требует пропускать все эле- менты, которые находятся между двумя последовательными элементами в исходной матрице c. Поэтому выделение вектора-строки из матрицы требует большего времени, чем выделе- ние вектора-столбца последовательных элементов: это требует меньше перемещений памяти и лучше использует различные уровни кэша. Следовательно, если данный алгоритм может быть изменён так, чтобы он мог обращаться или обновлять матрицу постолбцово, то следует предпочесть эту ориентацию.
5.4
Оптимизированные библиотеки линейной алгебры
В этом разделе мы представляем библиотеки линейной алгебры, которые используются в Scilab’е. В следующем разделе мы представляем BLAS и LAPACK, которые являются
127

A
11
A
21
A
m1
A
12
A
22
A
m2
A
1n
A
2n
A
mn
Рис. 37: Хранение матрицы чисел типа double размером m × n в Scilab. Два последователь- ных элемента в одном и том же столбце матрицы хранятся последовательно в памяти. Два последовательных элемента в одной строке матрицы разделены в памяти m адресами.
строительными блоками линейной алгебры в Scilab’е. Мы также представляем численные библиотеки ATLAS и Intel MKL. Мы представляем метод установки этих оптимизированных библиотек линейной алгебры в операционных системах Windows и Linux. Мы представляем влияние изменения библиотеки на производительность произведения общих неразреженных матриц.
5.4.1
BLAS, LAPACK, ATLAS и MKL
В этом разделе мы вкратце представляем библиотеки BLAS, LAPACK, ATLAS и MKL
для того, чтобы иметь общее представление об этих библиотеках. Эти библиотеки исполь- зуются Scilab’ом для обеспечения максимальной производительности в каждой системе.
Библиотека BLAS (Basic Linear Algebra Subprograms — основные подпрограммы ли- нейной алгебры) [
2
] является коллекцией процедур Фортрана, предоставляющих низкоуров- невые операции, такие, как добавление векторов, поэлементное и матричное умножения.
Библиотека BLAS делится на три уровня.
• BLAS уровня 1 выполняет скалярные, векторные и векторно-векторные операции.