밀도 매트릭스 - 벡터 곱셈의 가장 효율적인 구현은 무엇입니까?

Question

M이 조밀 한 m × n 매트릭스이고 v이 n 성분 벡터 인 경우, u = Mv은 u[i] = sum(M[i,j] * v[j], 1 <= j <= n)으로 주어진 m 성분 벡터입니다. 이 승산의 하나의 간단한 구현은 한번에 벡터 u 이상의 원소를 축적

allocate m-component vector u of zeroes 
for i = 1:m 
    for j = 1:n 
    u[i] += M[i,j] * v[j] 
    end 
end 

이다. 또 다른 구현은 루프를 교환하는 것입니다 :

allocate m-component vector u of zeroes 
for j = 1:n 
    for i = 1:m 
    u[i] += M[i,j] * v[j] 
    end 
end 

전체 벡터가 함께 구성됩니다.

효율적인 계산을 위해 설계된 C 및 Fortran과 같은 언어에서 이러한 구현 (일반적으로) 중 어느 것을 사용합니까? 내 생각에 C와 같은 언어는 내부적으로 열 우선 순위의 행렬을 저장하고, Fortran과 같은 언어는 열 우선 순위를 사용하여 후자를 사용하므로 내부 루프가 행렬 M의 연속 메모리 사이트에 액세스합니다. 이게 맞습니까? 메모리 위치가 기록되기 때문에

이전 구현

는 후자의 구현에 메모리 위치가 변경 m*n 번 기록되는 동안 만 만 m 독특한 위치가 어느 기록에도 불구하고, m 시간을 변경,보다 효율적으로 보인다. (물론 동일한 논리에 의해, 후자의 구현은 행 - 벡터 - 행렬 곱셈에 더 효율적일 것입니다. 그러나 이것은 훨씬 덜 일반적입니다.) 한편, Fortran은 일반적으로 조밀 한 행렬 벡터에서 더 빠릅니다 C보다 곱하기 때문에 아마도 (a) 자신의 구현을 잘못 추측하거나 (b) 두 구현의 상대적인 효율성을 잘못 판단하고있는 것일 수 있습니다.

Answer 1

확립 된 BLAS 구현을 사용하는 것이 가장 일반적입니다. 그 외에도 살펴볼 흥미로운 간단한 구현에 몇 가지 문제가 있습니다. 예를 들어, C (또는 그 문제에 대한 C++)에 종종 에일리어싱 포인터 최적화 많이 방지하고, 따라서 for example

void matvec(double *M, size_t n, size_t m, double *v, double * u) 
{ 
    for (size_t i = 0; i < m; i++) { 
     for (size_t j = 0; j < n; j++) { 
     u[i] += M[i * n + j] * v[j]; 
     } 
    } 
} 

는

.LBB0_4: # Parent Loop BB0_3 Depth=1 
    vmovsd xmm1, qword ptr [rcx + 8*rax] # xmm1 = mem[0],zero 
    vfmadd132sd xmm1, xmm0, qword ptr [r13 + 8*rax - 24] 
    vmovsd qword ptr [r8 + 8*rbx], xmm1 
    vmovsd xmm0, qword ptr [rcx + 8*rax + 8] # xmm0 = mem[0],zero 
    vfmadd132sd xmm0, xmm1, qword ptr [r13 + 8*rax - 16] 
    vmovsd qword ptr [r8 + 8*rbx], xmm0 
    vmovsd xmm1, qword ptr [rcx + 8*rax + 16] # xmm1 = mem[0],zero 
    vfmadd132sd xmm1, xmm0, qword ptr [r13 + 8*rax - 8] 
    vmovsd qword ptr [r8 + 8*rbx], xmm1 
    vmovsd xmm0, qword ptr [rcx + 8*rax + 24] # xmm0 = mem[0],zero 
    vfmadd132sd xmm0, xmm1, qword ptr [r13 + 8*rax] 
    vmovsd qword ptr [r8 + 8*rbx], xmm0 
    add rax, 4 
    cmp r11, rax 
    jne .LBB0_4 

연타 5 (내부 루프 발췌)하여 설정되어 있습니다

정말보기가 아플뿐 아니라 실행하는 데 더 많은 상처를줍니다. 컴파일러는 u이 M 및/또는 v으로 별칭을 지정할 수 있으므로 u에 저장하는 것이 큰 의구심으로 취급됩니다 (컴파일러에서 앨리어싱 테스트를 삽입 할 수 있었고 " 좋은 경우 경로). Fortran에서 기본적으로 프로 시저 인수는 별칭이 될 수 없으므로이 문제는 없었습니다. 이것은 특별한 트릭없이 임의로 입력 된 코드가 C보다 Fortran에서 더 빠른 전형적인 이유입니다 - 나머지는 내 대답이 아닙니다. C 코드를 좀 더 빠르게 만들려고합니다. 결국 나는 열 - 주요 M)로 돌아갑니다. C에서 앨리어싱 문제는 be fixed이 restrict과 같을 수 있지만 방해하지 않는 것이 유일한 것입니다 (u[i]에 합산하는 대신 명시 적 누적기를 사용하면 트릭을 수행하지만 마술 키워드는 사용하지 않아도됩니다).)

void matvec(double *M, size_t n, size_t m, double *v, double * restrict u) 
{ 
    for (size_t i = 0; i < m; i++) { 
     for (size_t j = 0; j < n; j++) { 
     u[i] += M[i * n + j] * v[j]; 
     } 
    } 
} 

지금이 발생 : 그 좋은, 그래서 그것은, 더 이상 스칼라 아니에요

.LBB0_8: # Parent Loop BB0_3 Depth=1 
    vmovupd ymm5, ymmword ptr [rcx + 8*rbx] 
    vmovupd ymm6, ymmword ptr [rcx + 8*rbx + 32] 
    vmovupd ymm7, ymmword ptr [rcx + 8*rbx + 64] 
    vmovupd ymm8, ymmword ptr [rcx + 8*rbx + 96] 
    vfmadd132pd ymm5, ymm1, ymmword ptr [rax + 8*rbx - 224] 
    vfmadd132pd ymm6, ymm2, ymmword ptr [rax + 8*rbx - 192] 
    vfmadd132pd ymm7, ymm3, ymmword ptr [rax + 8*rbx - 160] 
    vfmadd132pd ymm8, ymm4, ymmword ptr [rax + 8*rbx - 128] 
    vmovupd ymm1, ymmword ptr [rcx + 8*rbx + 128] 
    vmovupd ymm2, ymmword ptr [rcx + 8*rbx + 160] 
    vmovupd ymm3, ymmword ptr [rcx + 8*rbx + 192] 
    vmovupd ymm4, ymmword ptr [rcx + 8*rbx + 224] 
    vfmadd132pd ymm1, ymm5, ymmword ptr [rax + 8*rbx - 96] 
    vfmadd132pd ymm2, ymm6, ymmword ptr [rax + 8*rbx - 64] 
    vfmadd132pd ymm3, ymm7, ymmword ptr [rax + 8*rbx - 32] 
    vfmadd132pd ymm4, ymm8, ymmword ptr [rax + 8*rbx] 
    add rbx, 32 
    add rbp, 2 
    jne .LBB0_8 

. 그러나 이상적이지 않습니다.여기에 8 개의 FMA가 있지만 4 개의 종속 FMA 쌍으로 배열됩니다. 전체 루프를 통해 취해진 FMA의 독립적 종속성 체인은 실제로 4 개뿐입니다. FMA는 일반적으로 Skillake에서 대기 시간이 4이고 처리량이 2/사이클이므로 긴 대기 시간과 적절한 처리량을 갖기 때문에 모든 계산 처리량을 활용하려면 8 개의 독립적 인 FMA 체인이 필요합니다. Haswell은 더욱 심각합니다. FMA는 대기 시간이 5이고 이미 처리량이 2/cycle이므로 10 개의 독립적 인 체인이 필요했습니다. 다른 문제는 실제로 모든 FMA를 실제로 공급하기 어렵다는 것입니다. 위의 구조는 실제로 시도하지도 않습니다. FMA 당 2 개의로드를 사용하는 반면로드는 실제로 FMA와 동일한 처리량을 갖기 때문에 비율은 1 : 1이어야합니다.

부하 개선 : FMA 비율은 우리가이도 캐싱 고려하지 않습니다 (M의 여러 행에 대해 v에서 부하를 다시 사용 할 수있는 외부 루프를 줄이기 수행 할 수 있지만, 그것을 위해 도움이 , 너무). 바깥 고리를 푸는 것은 FMA의 더 독립적 인 사슬을 가지고 있다는 목표로도 작용합니다. 컴파일러는 일반적으로 안쪽 루프를 제외한 모든 것을 언로 드하는 것을 좋아하지 않으므로 약간의 시간이 걸립니다. manual work. "꼬리"반복 생략 (또는 : m이 4의 배수라고 가정).

void matvec(double *M, size_t n, size_t m, double *v, double * restrict u) 
{ 
    size_t i; 
    for (i = 0; i + 3 < m; i += 4) { 
     for (size_t j = 0; j < n; j++) { 
     size_t it = i; 
     u[it] += M[it * n + j] * v[j]; 
     it++; 
     u[it] += M[it * n + j] * v[j]; 
     it++; 
     u[it] += M[it * n + j] * v[j]; 
     it++; 
     u[it] += M[it * n + j] * v[j]; 
     } 
    } 
} 

불행하게도, 연타는 여전히 "틀린"그 순진한 직렬 풀다 인 상태, 내부 루프의 잘못을 풀다하기로 결정합니다.

void matvec(double *M, size_t n, size_t m, double *v, double *u) 
{ 
    size_t i; 
    for (i = 0; i + 3 < m; i += 4) { 
     double t0 = 0, t1 = 0, t2 = 0, t3 = 0; 
     for (size_t j = 0; j < n; j++) { 
     size_t it = i; 
     t0 += M[it * n + j] * v[j]; 
     it++; 
     t1 += M[it * n + j] * v[j]; 
     it++; 
     t2 += M[it * n + j] * v[j]; 
     it++; 
     t3 += M[it * n + j] * v[j]; 
     } 
     u[i] += t0; 
     u[i + 1] += t1; 
     u[i + 2] += t2; 
     u[i + 3] += t3; 
    } 
} 

이제 연타이 수행합니다 :

우리가 게으른 그만 마지막으로 일부 explicit accumulators을 할 경우

가

.LBB0_8: # Parent Loop BB0_3 Depth=1 
    vmovupd ymm5, ymmword ptr [rcx + 8*rdx] 
    vmovupd ymm6, ymmword ptr [rcx + 8*rdx + 32] 
    vfmadd231pd ymm4, ymm5, ymmword ptr [r12 + 8*rdx - 32] 
    vfmadd231pd ymm3, ymm5, ymmword ptr [r13 + 8*rdx - 32] 
    vfmadd231pd ymm2, ymm5, ymmword ptr [rax + 8*rdx - 32] 
    vfmadd231pd ymm1, ymm5, ymmword ptr [rbx + 8*rdx - 32] 
    vfmadd231pd ymm4, ymm6, ymmword ptr [r12 + 8*rdx] 
    vfmadd231pd ymm3, ymm6, ymmword ptr [r13 + 8*rdx] 
    vfmadd231pd ymm2, ymm6, ymmword ptr [rax + 8*rdx] 
    vfmadd231pd ymm1, ymm6, ymmword ptr [rbx + 8*rdx] 
    add rdx, 8 
    add rdi, 2 
    jne .LBB0_8 

이 문제는 사라집니다 : 여전히 4 독립적 인 체인이 있기 때문에 많은 점은 오래가 아니다

.LBB0_6: # Parent Loop BB0_3 Depth=1 
    vmovupd ymm8, ymmword ptr [r10 - 32] 
    vmovupd ymm9, ymmword ptr [r10] 
    vfmadd231pd ymm6, ymm8, ymmword ptr [rdi] 
    vfmadd231pd ymm7, ymm9, ymmword ptr [rdi + 32] 
    lea rax, [rdi + r14] 
    vfmadd231pd ymm4, ymm8, ymmword ptr [rdi + 8*rsi] 
    vfmadd231pd ymm5, ymm9, ymmword ptr [rdi + 8*rsi + 32] 
    vfmadd231pd ymm1, ymm8, ymmword ptr [rax + 8*rsi] 
    vfmadd231pd ymm3, ymm9, ymmword ptr [rax + 8*rsi + 32] 
    lea rax, [rax + r14] 
    vfmadd231pd ymm0, ymm8, ymmword ptr [rax + 8*rsi] 
    vfmadd231pd ymm2, ymm9, ymmword ptr [rax + 8*rsi + 32] 
    add rdi, 64 
    add r10, 64 
    add rbp, -8 
    jne .LBB0_6 

괜찮습니다. 하중 : FMA 비율이 10 : 8이고 Haswell 용 축전지가 너무 적어 일부 개선이 여전히 가능합니다. 다른 흥미로운 언롤링 조합은 (외부 x 내부) 4x3 (12 축 압기, 3 임시, 5/4로드 : FMA), 5x2 (10,2,6/5), 7x2 (14,2,8/7), 15x1 (15, 1, 16/15). 그렇게하면 외부 루프가 더 잘 풀리는 것처럼 보이지만 너무 많은 다른 스트림 ("스트리밍로드"의 의미에서 "스트리밍"이 아닐지라도)이 자동 프리 페치에는 좋지 않으며 실제로 스트리밍하면 채우기 버퍼의 수를 초과합니다 (실제 세부 사항은 부족합니다). 수동 프리 페치도 옵션입니다. 실제로 좋은 MVM 절차를 얻으려면 이러한 많은 작업을 수행하면서 훨씬 많은 작업이 필요합니다.

스토어를 u에 저장하면 내부 루프 외부에서 restrict이 더 이상 필요하지 않습니다. 가장 인상적인 것은, 이것이 생각하기에는 SIMD 내장 함수가 필요 없다는 것입니다. Clang은 무서운 잠재적 앨리어싱이 없다면 꽤 좋습니다. GCC와 ICC는 시도하지만 충분하게 전개하지는 않지만 더 많은 수작업 풀림이 트릭을 수행합니다.

루프 타일링도 옵션이지만 MVM입니다. 타일링은 MMM에 매우 필요하지만 MMM에는 MVM에없는 데이터 재사용이 거의 무제한입니다. 벡터 만 재사용되며 매트릭스는 한 번만 스트리밍됩니다. 거대한 행렬을 스트리밍하는 메모리 처리량은 캐시에 적합하지 않은 벡터보다 큰 문제가됩니다.

주요 열 M을 사용하면 루프 전달 종속성이 크게 달라집니다. 메모리를 통해 의존성이 있지만 시간이 많이 걸립니다. 하중 : FMA 비율은 여전히 ​​감소해야합니다. 따라서 외부 루프의 일부 풀리기가 여전히 필요하지만 전반적으로 다루기가 더 쉽습니다. 대부분 추가 ​​기능을 사용하도록 재정렬 할 수 있지만 FMA는 어쨌든 높은 처리량을 제공합니다 (HSW에서 추가보다 높습니다).그것은 짜증나는 수평 합계를 필요로하지 않지만 어쨌든 내부 루프 외부에서 발생했습니다. 그 대신 내부 루프에 상점이 있습니다. 이를 시도하지 않고서도 이러한 접근 방식간에 고유 한 차이가 크지는 않을 것으로 예상됩니다. 두 방법 모두 컴퓨팅 처리량 (캐시 가능한 크기의 경우)의 80 ~ 90 %까지 조정할 수 있어야합니다. "성가신 추가로드"는 본질적으로 어느 쪽이든 100 % 가까이 임의로 접근하지 못하게합니다.