разминка
В качестве разминки рассмотрим два более простых примера.
## example 1
x <- 1:11
for (i in 1:10) x[i] <- x[i + 1]
x
# [1] 2 3 4 5 6 7 8 9 10 11 11
x <- 1:11
x[1:10] <- x[2:11]
x
# [1] 2 3 4 5 6 7 8 9 10 11 11
## example 2
x <- 1:11
for (i in 1:10) x[i + 1] <- x[i]
x
# [1] 1 1 1 1 1 1 1 1 1 1 1
x <- 1:11
x[2:11] <- x[1:10]
x
# [1] 1 1 2 3 4 5 6 7 8 9 10
«Векторизация» успешна в 1-м примере, но не во 2-м.Почему?
Вот разумный анализ.«Векторизация» начинается с развертывания цикла, а затем выполняет несколько инструкций параллельно.Возможность векторизации цикла зависит от зависимости данных, переносимой циклом.
Развертывание цикла в примере 1 дает
x[1] <- x[2]
x[2] <- x[3]
x[3] <- x[4]
x[4] <- x[5]
x[5] <- x[6]
x[6] <- x[7]
x[7] <- x[8]
x[8] <- x[9]
x[9] <- x[10]
x[10] <- x[11]
Выполнение этих инструкций один за другим и выполнение их одновременнодать одинаковый результат.Таким образом, этот цикл можно «векторизовать».
Цикл в примере 2 - это
x[2] <- x[1]
x[3] <- x[2]
x[4] <- x[3]
x[5] <- x[4]
x[6] <- x[5]
x[7] <- x[6]
x[8] <- x[7]
x[9] <- x[8]
x[10] <- x[9]
x[11] <- x[10]
К сожалению, выполнение этих инструкций одна за другой и выполнение их одновременно не даст одинакового результата.Например, когда они выполняются один за другим, x[2] изменяется в 1-й инструкции, затем это измененное значение передается в x[3] во 2-й инструкции.Так что x[3] будет иметь то же значение, что и x[1].Однако при параллельном выполнении x[3] равняется x[2].В результате этот цикл не может быть «векторизован».
В теории «векторизации»
- Пример 1 имеет зависимость «запись после чтения» в данных:
x[i] изменяется после чтения; - Пример 2 имеет зависимость «чтение после записи» в данных:
x[i] читается после изменения.
Aцикл с зависимостью данных «запись после чтения» может быть «векторизован», а цикл с зависимостью данных «чтение после записи» не может.
в глубину
Возможно, многие люди были смущены сейчас.«Векторизация» - это «параллельная обработка»?
Да.В 1960-х годах, когда люди задавались вопросом, какой компьютер параллельной обработки предназначен для высокопроизводительных вычислений, Флинн разделил идеи дизайна на 4 типа.Категория «SIMD» (одна инструкция, несколько данных) называется «векторизацией», а компьютер с возможностью «SIMD» называется «векторным процессором» или «процессором массива».
В 1960-х годах не быломного языков программирования.Люди писали ассемблер (затем FORTRAN, когда был изобретен компилятор) для программирования регистров процессора напрямую«SIMD» компьютер может загружать несколько данных в векторный регистр с помощью одной инструкции и выполнять одинаковую арифметику с этими данными одновременно.Таким образом, обработка данных действительно параллельна.Рассмотрим наш пример 1 еще раз.Предположим, что векторный регистр может содержать два векторных элемента, тогда цикл можно выполнить с 5 итерациями, используя векторную обработку, а не 10 итераций, как при скалярной обработке.
reg <- x[2:3] ## load vector register
x[1:2] <- reg ## store vector register
-------------
reg <- x[4:5] ## load vector register
x[3:4] <- reg ## store vector register
-------------
reg <- x[6:7] ## load vector register
x[5:6] <- reg ## store vector register
-------------
reg <- x[8:9] ## load vector register
x[7:8] <- reg ## store vector register
-------------
reg <- x[10:11] ## load vector register
x[9:10] <- reg ## store vector register
Сегодня существует много языков программирования, таких как R .«Векторизация» больше не означает «SIMD». R - это не язык, на котором мы можем программировать регистры процессора.«Векторизация» в R - это просто аналогия «SIMD».В предыдущем вопросе и ответах: Означает ли термин "векторизация" разные вещи в разных контекстах? Я пытался объяснить это.Следующая карта иллюстрирует, как проводится эта аналогия:
single (assembly) instruction -> single R instruction
CPU vector registers -> temporary vectors
parallel processing in registers -> C/C++/FORTRAN loops with temporary vectors
Итак, R "векторизация" цикла в примере 1 выглядит примерно так:
## the C-level loop is implemented by function "["
tmp <- x[2:11] ## load data into a temporary vector
x[1:10] <- tmp ## fill temporary vector into x
Большую часть времени мы простоделать
x[1:10] <- x[2:10]
без явного присвоения временного вектора переменной.На созданный временный блок памяти не указывает ни одна переменная R, и поэтому он подлежит сборке мусора.
полная картинка
В приведенном выше разделе «векторизация» не вводитсяс самым простым примером.Очень часто «векторизация» вводится с чем-то вроде
a[1] <- b[1] + c[1]
a[2] <- b[2] + c[2]
a[3] <- b[3] + c[3]
a[4] <- b[4] + c[4]
, где a, b и c не являются псевдонимами в памяти, то есть блоки памяти, хранящие векторы a,b и c не перекрываются.Это идеальный случай, так как отсутствие псевдонимов памяти не подразумевает зависимости от данных.
Помимо «зависимости от данных», существует также «зависимость управления», то есть работа с «если ... еще ...»"в" векторизации ".Однако по причине времени и пространства я не буду подробно останавливаться на этом вопросе.
вернуться к примеру в вопросе
Теперь пришло время исследовать цикл в вопросе.
set.seed(0)
x <- round(runif(10), 2)
sig <- sample.int(10)
# [1] 1 2 9 5 3 4 8 6 7 10
for (i in seq_along(sig)) x[i] <- x[sig[i]]
Развертывание цикла дает
x[1] <- x[1]
x[2] <- x[2]
x[3] <- x[9] ## 3rd instruction
x[4] <- x[5]
x[5] <- x[3] ## 5th instruction
x[6] <- x[4]
x[7] <- x[8]
x[8] <- x[6]
x[9] <- x[7]
x[10] <- x[10]
Между 3-й и 5-й инструкциями существует зависимость данных «чтение после записи», поэтому цикл нельзя «векторизовать» (см. Замечание 1 ).
Что же тогдаx[] <- x[sig] делает?Давайте сначала явно выпишем временный вектор:
tmp <- x[sig]
x[] <- tmp
Поскольку "[" вызывается дважды, за этим "векторизованным" кодом на самом деле находятся два цикла уровня C:
tmp[1] <- x[1]
tmp[2] <- x[2]
tmp[3] <- x[9]
tmp[4] <- x[5]
tmp[5] <- x[3]
tmp[6] <- x[4]
tmp[7] <- x[8]
tmp[8] <- x[6]
tmp[9] <- x[7]
tmp[10] <- x[10]
x[1] <- tmp[1]
x[2] <- tmp[2]
x[3] <- tmp[3]
x[4] <- tmp[4]
x[5] <- tmp[5]
x[6] <- tmp[6]
x[7] <- tmp[7]
x[8] <- tmp[8]
x[9] <- tmp[9]
x[10] <- tmp[10]
Таким образом, x[] <- x[sig] эквивалентно
for (i in 1:10) tmp[i] <- x[sig[i]]
for (i in 1:10) x[i] <- tmp[i]
rm(tmp); gc()
, что совсем не соответствует первоначальному циклу, указанному в вопросе.
Замечание 1
При реализации циклав Rcpp рассматривается как "векторизация", тогда пусть это будет.Но нет никакой возможности в дальнейшем «векторизовать» цикл C / C ++ с помощью «SIMD».
Замечание 2
Эти вопросы и ответы мотивированы этими вопросамиA .OP первоначально представил цикл
for (i in 1:num) {
for (j in 1:num) {
mat[i, j] <- mat[i, mat[j, "rm"]]
}
}
Заманчиво «векторизовать» его как
mat[1:num, 1:num] <- mat[1:num, mat[1:num, "rm"]]
, но это потенциально неправильно.Позже OP изменил цикл на
for (i in 1:num) {
for (j in 1:num) {
mat[i, j] <- mat[i, 1 + num + mat[j, "rm"]]
}
}
, что устраняет проблему псевдонимов памяти, поскольку заменяемые столбцы - это первые num столбцы, а столбцы, которые нужно искать, - после первого num.столбцы.
Замечание 3
Я получил несколько комментариев относительно того, выполняет ли цикл в вопросе "11-местную" модификацию x.Да, это.Мы можем использовать tracemem:
set.seed(0)
x <- round(runif(10), 2)
sig <- sample.int(10)
tracemem(x)
#[1] "<0x28f7340>"
for (i in seq_along(sig)) x[i] <- x[sig[i]]
tracemem(x)
#[1] "<0x28f7340>"
Мой сеанс R выделил блок памяти, указанный адресом <0x28f7340> для x, и вы можете увидеть другое значение при запуске кода.Однако вывод tracemem не изменится после цикла, что означает, что копия x не создается.Таким образом, цикл действительно выполняет модификацию «на месте» без использования дополнительной памяти.
Однако цикл не выполняет перестановку «на месте».Перестановка на месте - более сложная операция.Меняются не только элементы x вдоль цикла, но и элементы sig также должны быть поменяны местами (и, в конце концов, sig будет 1:10).