Рассмотрим некоторые приёмы повышения быстродействия вычислительных программ на примере алгоритма медианного фильтра.
Медианный фильтр – алгоритм, применяемый к изображению или подобной ему матрице, который значение каждого элемента заменяет на медиану среди его соседей. Сам по себе он может применяться для грубого сглаживания изображения, однако более полезен для устранения битых и горячих пикселей. Идея тут такая, что, если значение исходного элемента очень сильно отличается от значения медианного фильтра в той же точке, то значение элемента заменяется значением медианного фильтра.
В качестве исходного медианного фильтра используем свободно распространяемый исходный код профессора факультета математики Калифорнийского технологического университета Джона Буркардта. Этот код (gray_median_news) входит в модуль image_denoise и представлен на языках C, C++ и Fortran 90. Мы здесь рассмотрим оптимизацию кода на Фортране, однако практически всё то же самое применимо и к сишным версиям, различия чисто синтаксические.
Исходная процедура gray_median_news отлично работает (хотя содержит небольшие ошибки при обработке нижней строки), однако довольно медленна. Посмотрим, что тут можно сделать.
Взглянем на главный цикл обработки данных в фильтре:
Во-первых, тут напрашивается единственная чисто фортрановская оптимизация. Данные в фортрановских массивах хранятся по столбцам, а циклы обходят их по строкам, поэтому для более эффективного регулярного доступа к памяти поменяем внутренний и внешний циклы местами:
Для программ на Си, естественно, надо делать наоборот.
Эта простая перестановка циклов сама по себе дала нам ускорение примерно в 2 раза.
Далее, обратим внимание, что обработка здесь заключается в пересчитывании данных из массива gray в массив gray2 и асинхронна по отношению к порядку элементов. Поэтому хотелось бы распараллелить и заасинхронить эти циклы. Но нам здесь будет мешать массив p. Конечно, можно его сделать приватным для разных веток, но вместо этого давайте сначала разберёмся с подпрограммой i4vec_median.
Что же профессор делает в i4vec_median? А вот что:
Здесь он вызывает другую подпрограмму i4vec_frac, предназначенную для поиска k-го элемента отсортированного массива a длиной n, полагая k равным середине массива. Сама подпрограмма i4vec_frac выглядит так:
Разумеется, учинять сортировку массива произвольного размера – это совсем не то, что бы нам хотелось иметь для решения вполне конкретной задачи нахождения медианы из пяти чисел (то есть такого числа, что два других не меньше его, а ещё два других не больше). Дальше по ходу алгоритма профессора ещё используются медианы тройки чисел на краях матрицы.
Попробуем подумать, как нам найти медиану пятёрки и медиану тройки наиболее эффективным способом.
Первое, что приходит в голову – это тупо расписать все альтернативы упорядоченности чисел вложенными условными операторами. Например, для тройки это может выглядеть так:
К сожалению, для пятёрки таким путём получится монструозная конструкция из 120 веток вложенных условных операторов, которую написать представляется возможным разве что путём автоматической генерации. Но надо ли это делать?
Вспомним, что современные процессоры включают наборы векторных инструкций SSE и подобные им, которые обрабатывают небольшие м��ссивы гораздо эффективнее, чем процессор применял бы последовательные ветвления. Беглый анализ существующей практики показывает, что наиболее эффективными процедурами, использующими SSE, в нашем случае будут такие:
Можно практически убедиться, что функция median3, использующая max и min (компилируемые один-в-один в машинные инструкции вроде vpmaxsd и vpminsd), действительно примерно вдвое быстрее, чем использующая if.
Перепишем наши циклы с использованием процедур median3 и median5:
Обратите внимание, что мы избавились от массива p, так как можем передать параметры просто в 5 аргументах функции median5.
Замена сортировки и поиска k-го элемента на простые минимаксные функции медианы для 3 и 5 элементов дала нам ускорение примерно в 50 раз.
И теперь нам ничего не мешает, действуя чисто формально, распараллелить и заасинхронить циклы обработки:
Аналогичным образом мы можем исправить и менее нагруженные циклы обработки точек на краях матрицы.
Почему мы параллелим внешний цикл и не параллелим внутренний? Потому что параллельности внешнего цикла вполне достаточно при разумном для обычного компьютера предположении о том, что количество столбцов в массиве превышает количество ядер процессоров. А накладные расходы на организацию потоков вычислений лучше нести однократно, для внешнего цикла.
Распараллеливание в данном случае дало нам небольшое ускорение на тестовой настольной машине, примерно в полтора раза. В данном случае мы упёрлись в пропускную способность памяти и не смогли реализовать преимущество от одновременной работы нескольких ядер. На сервере с многоканальной памятью мы получаем более значительный эффект.
Каков же общий итог этих оптимизаций? На тестовых настольных компьютерах автора работу фильтра удалось ускорить в 130–160 раз. Для матрицы в 4 миллиона элементов это соответствует прогрессу от десятых до тысячных долей секунды. Фильтр стал пригоден, например, для обработки кадров в реальном времени.
Разумеется, быстрые алгоритмы фильтрации изображения имеются в специализированных библиотеках, так что эту статью не следует рассматривать как руководство по обработке изображений. Идея статьи в том, чтобы показать некоторые приёмы оптимизации и проиллюстрировать важность периодического подключения к процессу программирования мозга, а не только скачанных из интернета библиотек.
Практические выводы:
Медианный фильтр – алгоритм, применяемый к изображению или подобной ему матрице, который значение каждого элемента заменяет на медиану среди его соседей. Сам по себе он может применяться для грубого сглаживания изображения, однако более полезен для устранения битых и горячих пикселей. Идея тут такая, что, если значение исходного элемента очень сильно отличается от значения медианного фильтра в той же точке, то значение элемента заменяется значением медианного фильтра.
В качестве исходного медианного фильтра используем свободно распространяемый исходный код профессора факультета математики Калифорнийского технологического университета Джона Буркардта. Этот код (gray_median_news) входит в модуль image_denoise и представлен на языках C, C++ и Fortran 90. Мы здесь рассмотрим оптимизацию кода на Фортране, однако практически всё то же самое применимо и к сишным версиям, различия чисто синтаксические.
Исходная процедура gray_median_news отлично работает (хотя содержит небольшие ошибки при обработке нижней строки), однако довольно медленна. Посмотрим, что тут можно сделать.
Взглянем на главный цикл обработки данных в фильтре:
do i = 2, m - 1 do j = 2, n - 1 p(1) = gray(i-1,j ) p(2) = gray(i+1,j ) p(3) = gray(i ,j+1) p(4) = gray(i ,j-1) p(5) = gray(i ,j ) call i4vec_median ( 5, p, gray2(i,j) ) end do end do
Во-первых, тут напрашивается единственная чисто фортрановская оптимизация. Данные в фортрановских массивах хранятся по столбцам, а циклы обходят их по строкам, поэтому для более эффективного регулярного доступа к памяти поменяем внутренний и внешний циклы местами:
do j = 2, n - 1 do i = 2, m - 1 p(1) = gray(i-1,j ) p(2) = gray(i+1,j ) p(3) = gray(i ,j+1) p(4) = gray(i ,j-1) p(5) = gray(i ,j ) call i4vec_median ( 5, p, gray2(i,j) ) end do end do
Для программ на Си, естественно, надо делать наоборот.
Эта простая перестановка циклов сама по себе дала нам ускорение примерно в 2 раза.
Далее, обратим внимание, что обработка здесь заключается в пересчитывании данных из массива gray в массив gray2 и асинхронна по отношению к порядку элементов. Поэтому хотелось бы распараллелить и заасинхронить эти циклы. Но нам здесь будет мешать массив p. Конечно, можно его сделать приватным для разных веток, но вместо этого давайте сначала разберёмся с подпрограммой i4vec_median.
Что же профессор делает в i4vec_median? А вот что:
subroutine i4vec_median ( n, a, median ) implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) k integer ( kind = 4 ) median k = ( n + 1 ) / 2 call i4vec_frac ( n, a, k, median ) return end
Здесь он вызывает другую подпрограмму i4vec_frac, предназначенную для поиска k-го элемента отсортированного массива a длиной n, полагая k равным середине массива. Сама подпрограмма i4vec_frac выглядит так:
Скрытый текст
subroutine i4vec_frac ( n, a, k, frac ) implicit none integer ( kind = 4 ) n integer ( kind = 4 ) a(n) integer ( kind = 4 ) frac integer ( kind = 4 ) i integer ( kind = 4 ) iryt integer ( kind = 4 ) ix integer ( kind = 4 ) j integer ( kind = 4 ) k integer ( kind = 4 ) left integer ( kind = 4 ) t if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4VEC_FRAC - Fatal error!' write ( *, '(a,i8)' ) ' Illegal nonpositive value of N = ', n stop end if if ( k <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4VEC_FRAC - Fatal error!' write ( *, '(a,i8)' ) ' Illegal nonpositive value of K = ', k stop end if if ( n < k ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4VEC_FRAC - Fatal error!' write ( *, '(a,i8)' ) ' Illegal N < K, K = ', k stop end if left = 1 iryt = n do if ( iryt <= left ) then frac = a(k) exit end if ix = a(k) i = left j = iryt do if ( j < i ) then if ( j < k ) then left = i end if if ( k < i ) then iryt = j end if exit end if ! ! Find I so that IX <= A(I). ! do while ( a(i) < ix ) i = i + 1 end do ! ! Find J so that A(J) <= IX. ! do while ( ix < a(j) ) j = j - 1 end do if ( i <= j ) then t = a(i) a(i) = a(j) a(j) = t i = i + 1 j = j - 1 end if end do end do return end
Разумеется, учинять сортировку массива произвольного размера – это совсем не то, что бы нам хотелось иметь для решения вполне конкретной задачи нахождения медианы из пяти чисел (то есть такого числа, что два других не меньше его, а ещё два других не больше). Дальше по ходу алгоритма профессора ещё используются медианы тройки чисел на краях матрицы.
Попробуем подумать, как нам найти медиану пятёрки и медиану тройки наиболее эффективным способом.
Первое, что приходит в голову – это тупо расписать все альтернативы упорядоченности чисел вложенными условными операторами. Например, для тройки это может выглядеть так:
pure integer function median3 (a1, a2, a3) implicit none integer, intent(in) :: a1, a2, a3 if (a1 > a2) then if (a1 > a3) then if (a2 > a3) then median3 = a2 else median3 = a3 end if else median3 = a1 end if else if (a2 > a3) then median3 = a3 else median3 = a2 end if end function median3
К сожалению, для пятёрки таким путём получится монструозная конструкция из 120 веток вложенных условных операторов, которую написать представляется возможным разве что путём автоматической генерации. Но надо ли это делать?
Вспомним, что современные процессоры включают наборы векторных инструкций SSE и подобные им, которые обрабатывают небольшие м��ссивы гораздо эффективнее, чем процессор применял бы последовательные ветвления. Беглый анализ существующей практики показывает, что наиболее эффективными процедурами, использующими SSE, в нашем случае будут такие:
pure integer function median3 (a, b, c) implicit none integer, intent(in) :: a, b, c median3 = max(min(a,b),min(c,max(a,b))) end function median3 pure integer function median5 (a, b, c, d, e) implicit none integer, intent(in) :: a, b, c, d, e integer :: f, g f = max(min(a,b),min(c,d)) g = min(max(a,b),max(c,d)) median5 = median3 (e,f,g) end function median5
Можно практически убедиться, что функция median3, использующая max и min (компилируемые один-в-один в машинные инструкции вроде vpmaxsd и vpminsd), действительно примерно вдвое быстрее, чем использующая if.
Перепишем наши циклы с использованием процедур median3 и median5:
do j = 2, n - 1 do i = 2, m - 1 gray2(i,j) = median5 (gray(i-1,j ), & gray(i+1,j ), & gray(i ,j+1), & gray(i ,j-1), & gray(i ,j )) end do end do
Обратите внимание, что мы избавились от массива p, так как можем передать параметры просто в 5 аргументах функции median5.
Замена сортировки и поиска k-го элемента на простые минимаксные функции медианы для 3 и 5 элементов дала нам ускорение примерно в 50 раз.
И теперь нам ничего не мешает, действуя чисто формально, распараллелить и заасинхронить циклы обработки:
!$omp parallel do private (i) do j = 2, n - 1 do concurrent (i=2:m-1) gray2(i,j) = median5 (gray(i-1,j ), & gray(i+1,j ), & gray(i ,j+1), & gray(i ,j-1), & gray(i ,j )) end do end do !$omp end parallel do
Аналогичным образом мы можем исправить и менее нагруженные циклы обработки точек на краях матрицы.
Почему мы параллелим внешний цикл и не параллелим внутренний? Потому что параллельности внешнего цикла вполне достаточно при разумном для обычного компьютера предположении о том, что количество столбцов в массиве превышает количество ядер процессоров. А накладные расходы на организацию потоков вычислений лучше нести однократно, для внешнего цикла.
Распараллеливание в данном случае дало нам небольшое ускорение на тестовой настольной машине, примерно в полтора раза. В данном случае мы упёрлись в пропускную способность памяти и не смогли реализовать преимущество от одновременной работы нескольких ядер. На сервере с многоканальной памятью мы получаем более значительный эффект.
Каков же общий итог этих оптимизаций? На тестовых настольных компьютерах автора работу фильтра удалось ускорить в 130–160 раз. Для матрицы в 4 миллиона элементов это соответствует прогрессу от десятых до тысячных долей секунды. Фильтр стал пригоден, например, для обработки кадров в реальном времени.
Разумеется, быстрые алгоритмы фильтрации изображения имеются в специализированных библиотеках, так что эту статью не следует рассматривать как руководство по обработке изображений. Идея статьи в том, чтобы показать некоторые приёмы оптимизации и проиллюстрировать важность периодического подключения к процессу программирования мозга, а не только скачанных из интернета библиотек.
Практические выводы:
- если вас заботит быстродействие программы, не используйте общие алгоритмы вместо частных, избегайте ненужных генерализаций задачи;
- используйте встроенные функции max и min, они крайне эффективно транслируются в код современных процессоров;
- старайтесь не использовать лишние ветвления в глубоких циклах, чтобы не сбивать конвейер;
- по возможности распараллеливайте долгие циклы – с OpenMP это достаточно просто;
- избегайте ненужной взаимосвязи по данным между итерациями циклов;
- проверяйте, не имеет ли решаемая вами задача известного неочевидного, но эффективного решения;
- не думайте, что скачанная даже из солидного источника библиотека функций не может быть серьёзно оптимизирована.
