병렬 계산 - do concurrent 와 OpenMP

이 장에서 배우는 것

앞 장에서 만든 2차원 열 확산 계산은 현재 온도 배열을 읽어 다음 온도 배열을 채운다. 내부 격자점의 계산식은 같고, 서로 다른 격자점에 기록하는 결과도 겹치지 않는다. 이 구조는 반복을 여러 실행 흐름에 나누기 좋은 출발점이다. 그러나 반복문 앞에 병렬 지시문을 붙이는 것만으로 계산의 타당성이 보장되지는 않는다. 먼저 반복 사이의 의존성을 확인하고, 동시에 접근하는 데이터의 읽기와 쓰기를 구분해야 한다.

이 장에서는 같은 열 확산 계산을 순차 반복, do concurrent, OpenMP로 작성한다. 계산량을 늘리거나 실행 시간을 출력하는 대신, 작은 격자를 이용해 세 구현의 결과를 직접 비교한다. 병렬화의 첫 검증 대상은 속도가 아니라 계산 결과다.

• do concurrent가 요구하는 반복 독립성을 설명하고 적합한 반복을 구분한다.
• OpenMP의 parallel do와 reduction으로 격자 계산과 통계 계산을 나눈다.
• 현재 배열과 다음 배열을 분리하고, 공유 변수에 대한 데이터 경쟁을 피한다.
• 순차 결과와 배열 전체를 비교하며, 부동소수점 합산 순서의 영향을 해석한다.

문제 상황

금속판의 격자를 촘촘하게 만들면 한 시간 단계에서 계산할 격자점이 늘어난다. 시간 단계도 많아지므로, 내부 격자점을 순회하는 반복이 프로그램 실행의 큰 부분을 차지할 수 있다. 여러 코어를 가진 컴퓨터에서 이 반복을 나누어 실행하고 싶어지는 상황이다.

앞 장의 명시적 갱신식은 한 격자점의 현재 온도와 상하좌우 온도로 다음 온도를 구한다. 공간 간격이 두 방향에서 같을 때, 계산에 사용할 계수를 r이라고 두면 다음과 같이 쓸 수 있다. 여기서 r은 열 확산 계수와 시간 간격의 곱을 공간 간격의 제곱으로 나눈 값이다.

next(i,j) = current(i,j) + r * (current(i-1,j) + current(i+1,j) &
           + current(i,j-1) + current(i,j+1) - 4.0_real64 * current(i,j))

병렬화가 바꾸어야 하는 것은 격자점 계산을 맡기는 방식이다. 어느 시점의 온도를 읽는지, 경계 온도를 어떻게 유지하는지, 시간 단계를 언제 넘어가는지는 그대로 유지해야 한다. 한 단계의 내부 계산이 모두 끝나기 전에 다음 단계를 시작하면 다른 수치 방법이 된다.

이번 예제는 가로와 세로가 각각 다섯 점인 판을 사용한다. 가장자리 온도는 0으로 고정하고, 중앙 한 점의 초기 온도를 64로 둔다. r은 0.125이며 두 단계만 계산한다. 작은 크기는 성능 평가용이 아니라 각 구현의 계산과 출력을 손으로도 확인하기 위한 선택이다. 이 계수는 앞 장의 명시적 방법에서 사용하는 안정 조건을 만족한다.

do concurrent는 반복의 독립성을 표현한다

do concurrent는 반복마다 해야 할 일이 서로 의존하지 않는다는 뜻을 코드에 담는다. 일반 do의 실행 순서를 이용해 앞 반복의 결과를 다음 반복에서 읽는 코드는 이 구문에 맞지 않는다. 처리기가 반복을 어떤 순서로 수행하더라도 같은 계산이 성립해야 한다.

이 선언이 여러 코어의 사용을 보장하는 것은 아니다. 컴파일러는 순차 실행을 선택할 수도 있고, 벡터 명령 등을 이용할 수도 있다. 이 장의 컴파일 명령만으로 do concurrent가 여러 스레드에서 실행된다고 해석해서는 안 된다. 여기서 확인하는 것은 반복 독립성을 표현한 코드가 순차 구현과 같은 결과를 만드는지다.

열 확산에서는 현재 배열을 읽기 전용으로 취급하고 다음 배열의 각 원소를 한 번씩 기록하면 반복을 독립적으로 만들 수 있다. 서로 다른 반복이 같은 현재 원소를 읽어도 문제가 없다. 문제가 되는 것은 다른 반복이 읽는 원소를 동시에 바꾸거나, 여러 반복이 같은 결과 원소에 기록하는 경우다.

독립성 검토는 배열 이름만 보고 끝내지 않는다. 호출하는 함수가 모듈 변수에 값을 저장하거나, 내부에 유지되는 상태를 바꾸거나, 같은 파일에 출력을 남기는지도 살펴야 한다. do concurrent에서 호출하는 프로시저는 순수 프로시저여야 한다. 완성 코드에서는 한 점의 온도를 구하는 함수를 pure로 선언하고, 전달받은 배열과 인덱스만으로 값을 계산한다.

이 함수 안의 임시값은 각 호출에서 계산되는 값이다. 반복 바깥에 선언한 하나의 임시 변수에 계산값을 넣고 여러 반복이 사용하는 구조를 피할 수 있다. 순수 함수라는 선언만으로 잘못된 배열 갱신 방식이 고쳐지는 것은 아니지만, 계산에 필요한 입력과 결과의 관계를 좁혀서 검토하기 쉽게 한다.

Fortran 2018에서는 do concurrent에 지역성 지정자를 사용하는 방법도 있다. 이 장은 그 문법을 추가하기보다 반복 본문에 공유 임시값을 두지 않는 구조를 택한다. 합계 계산 역시 do concurrent 본문에서 하나의 변수에 누적하지 않는다. do concurrent는 반복을 독립적으로 표현하는 구문이지, 보통의 누적문을 자동으로 안전한 병렬 합산으로 바꾸는 구문이 아니다.

격자점 반복을 독립적으로 만들기 위한 접근 규칙

데이터	반복의 접근	확인할 조건
현재 온도 배열	여러 반복이 읽는다	계산 중 값을 바꾸지 않는다
다음 온도 배열	각 반복이 한 원소를 쓴다	서로 다른 반복의 기록 위치가 겹치지 않는다
임시 계산값	한 점의 계산에만 사용한다	다른 반복과 상태를 공유하지 않는다
시간 단계	앞 단계 결과를 읽는다	단계 사이의 실행 순서를 유지한다

OpenMP로 반복과 합산을 나눈다

OpenMP는 컴파일러 지시문을 사용해 병렬 실행을 지정하는 방식이다. Fortran 자유 형식 소스에서는 !$omp로 시작하는 줄을 사용한다. GNU Fortran에서 -fopenmp를 주면 지시문이 활성화되고, 해당 옵션을 주지 않으면 이 예제의 지시문은 주석으로 취급된다. 따라서 같은 소스로 순차 실행과 병렬 실행을 모두 만들 수 있다.

parallel do는 스레드(thread)들이 이어지는 반복의 일을 나누어 수행하도록 지정한다. 이 예제에서는 이중 반복에 collapse(2)를 붙여 두 반복의 조합을 작업 분배 대상으로 삼는다. 안쪽 반복은 바깥쪽 반복의 첫 실행문이어야 하며, 두 반복 사이에 일반 실행문을 끼우지 않는다. 두 방향의 반복 범위도 서로 독립적인 직사각형 영역으로 둔다.

schedule(static)은 반복 작업을 정적으로 나누는 방식이다. 각 격자점의 계산량이 비슷한 이번 코드에 적합한 선택이다. 이 지정은 결과의 정확성을 대신 보장하지 않는다. 실행 순서에 의존하는 반복은 작업 배분 방식을 바꾸어도 여전히 문제가 있다.

default(none)은 데이터의 공유 방식을 명시하도록 돕는다. shared로 지정한 현재 배열은 모든 스레드가 같은 저장 공간을 읽는다. 다음 배열도 공유하지만 각 반복이 쓰는 원소가 다르므로 기록이 충돌하지 않는다. 반복 인덱스는 private로 지정한다. 반복 제어 변수에는 별도의 규칙도 적용되지만, 여기서는 역할을 눈에 보이게 적는다.

합계는 격자점 갱신과 다르게 다루어야 한다. 모든 반복에서 total = total + 값을 실행하면 하나의 변수에 읽기와 쓰기가 겹친다. 환원(reduction)을 사용하면 각 스레드가 자체 누적값을 계산한 뒤 연산에 맞게 결과를 합친다. 완성 코드에서는 합계에 +를, 최고 온도에 max를 사용한다.

환원 대상인 원래 변수는 병렬 영역에 들어가기 전에 초기화한다. 합계의 시작값은 0으로 둔다. 최고 온도의 시작값은 -huge로 두어 음수 온도만 있는 배열에서도 올바른 최대값을 구할 수 있게 한다. 최종 결합에는 원래 변수의 시작값도 관여하므로, 이전 계산의 값을 남겨 두면 의도와 다른 결과가 된다.

parallel do 끝에는 기본적으로 동기화 지점이 있다. 이번 프로그램은 이 지점을 지난 뒤에만 다음 배열을 현재 배열에 복사한다. 시간 단계를 감싸는 바깥 반복은 순차적으로 진행된다. 공간의 독립성을 이용하면서 시간 방향의 의존성은 유지하는 구조다.

같은 결과를 확인하는 기준

데이터 경쟁(data race)은 서로 다른 스레드가 동기화 없이 같은 저장 위치에 접근하고, 그중 하나 이상이 기록하는 상황에서 생긴다. 결과가 가끔 맞는다고 경쟁이 없는 것은 아니다. 격자의 크기나 스레드 수를 바꾸면 잘못된 결과가 드러날 수 있다. 결과 비교와 함께 공유 위치에 대한 접근 구조를 검토해야 한다.

배열 비교에는 두 배열의 차이에 절댓값을 취한 뒤 최대값을 구한다. 중앙값이나 전체 합계만 비교하면 서로 다른 위치에서 생긴 오차가 가려질 수 있다. 순차 배열을 기준으로 do concurrent 배열과 OpenMP 배열을 각각 비교하면, 한 원소라도 다른지 확인할 수 있다.

일반적인 부동소수점 합산에서는 계산 순서가 바뀌면 마지막 자릿수가 달라질 수 있다. 환원은 안전한 합산을 제공하지만 순차 합산과 같은 비트 패턴까지 보장하지는 않는다. 이 차이를 데이터 경쟁과 혼동해서는 안 된다. 실제 계산에서는 값의 규모를 반영한 절대 허용 오차와 상대 허용 오차를 함께 검토한다.

이번 예제는 64, 0.125와 두 단계의 결과가 이진 부동소수점으로 정확히 표현되는 작은 값들이다. 합계도 정확히 표현 가능한 범위에 있다. 따라서 제시한 조건에서는 세 배열의 최대 차이가 0이고, 순차 합계와 환원 합계가 모두 60이다. 코드의 허용 오차 1.0e-12는 이 검증 예제에 맞춘 값이며, 다른 문제에 그대로 적용하는 일반 기준은 아니다.

0으로 고정한 경계는 열이 빠져나가는 조건이다. 초기 내부 온도의 합은 64지만 두 단계 뒤 내부 온도의 합은 60이 된다. 합계 감소 자체가 병렬화 오류를 뜻하지 않는다. 물리적 경계 조건이 허용하는 변화와 구현 오류를 구분해야 한다. 공간 간격과 물성값을 곱하지 않은 온도 합계는 여기서 배열 검증용 통계로 사용한다.

완성 코드

모든 프로시저를 한 프로그램의 내부 프로시저로 두었다. 입력 파일과 표준 입력은 필요하지 않다. 출력은 모든 병렬 반복이 끝난 뒤 주 프로그램에서만 수행한다.

main.f90

program heat_parallel
  use iso_fortran_env, only : real64
  implicit none

  integer, parameter :: nx = 5, ny = 5, nsteps = 2
  real(real64), parameter :: r = 0.125_real64
  real(real64), parameter :: tol = 1.0e-12_real64
  real(real64) :: initial(nx,ny)
  real(real64) :: seq(nx,ny), conc(nx,ny), par(nx,ny)
  real(real64) :: next_seq(nx,ny), next_conc(nx,ny), next_par(nx,ny)
  real(real64) :: diff_conc, diff_par
  real(real64) :: total_seq, total_par, peak_par
  integer :: step, j

  initial = 0.0_real64
  initial(3,3) = 64.0_real64
  seq = initial
  conc = initial
  par = initial

  do step = 1, nsteps
    call advance_serial(seq, next_seq)
    call advance_concurrent(conc, next_conc)
    call advance_openmp(par, next_par)
    seq = next_seq
    conc = next_conc
    par = next_par
  end do

  diff_conc = maxval(abs(seq - conc))
  diff_par = maxval(abs(seq - par))
  total_seq = sum(seq)
  call statistics_openmp(par, total_par, peak_par)

  if (diff_conc > tol .or. diff_par > tol) then
    error stop 'array comparison failed'
  end if
  if (abs(total_seq - total_par) > tol) then
    error stop 'sum comparison failed'
  end if
  if (abs(total_seq - 60.0_real64) > tol .or. &
      abs(peak_par - 20.0_real64) > tol) then
    error stop 'reference check failed'
  end if

  write(*,'(A,I2)') 'steps = ', nsteps
  write(*,'(A)') 'sequential temperature:'
  do j = ny, 1, -1
    write(*,'(5F8.2)') seq(:,j)
  end do
  write(*,'(A,F10.6)') 'max difference (do concurrent) = ', diff_conc
  write(*,'(A,F10.6)') 'max difference (OpenMP)        = ', diff_par
  write(*,'(A,F8.2)') 'sum (sequential) = ', total_seq
  write(*,'(A,F8.2)') 'sum (OpenMP)     = ', total_par
  write(*,'(A,F8.2)') 'maximum         = ', peak_par
  write(*,'(A)') 'verification: PASS'

contains

  pure function point_value(src, i, j) result(value)
    real(real64), intent(in) :: src(:,:)
    integer, intent(in) :: i, j
    real(real64) :: value

    value = src(i,j) + r * (src(i-1,j) + src(i+1,j) &
          + src(i,j-1) + src(i,j+1) - 4.0_real64 * src(i,j))
  end function point_value

  subroutine advance_serial(src, dst)
    real(real64), intent(in) :: src(:,:)
    real(real64), intent(out) :: dst(:,:)
    integer :: i, j

    dst = src
    do j = 2, size(src,2) - 1
      do i = 2, size(src,1) - 1
        dst(i,j) = point_value(src, i, j)
      end do
    end do
  end subroutine advance_serial

  subroutine advance_concurrent(src, dst)
    real(real64), intent(in) :: src(:,:)
    real(real64), intent(out) :: dst(:,:)
    integer :: i, j

    dst = src
    do concurrent (j = 2:size(src,2)-1, i = 2:size(src,1)-1)
      dst(i,j) = point_value(src, i, j)
    end do
  end subroutine advance_concurrent

  subroutine advance_openmp(src, dst)
    real(real64), intent(in) :: src(:,:)
    real(real64), intent(out) :: dst(:,:)
    integer :: i, j, nxi, nyi

    nxi = size(src,1)
    nyi = size(src,2)
    dst = src
    !$omp parallel do collapse(2) default(none) &
    !$omp& shared(src, dst, nxi, nyi) private(i, j) schedule(static)
    do j = 2, nyi - 1
      do i = 2, nxi - 1
        dst(i,j) = point_value(src, i, j)
      end do
    end do
    !$omp end parallel do
  end subroutine advance_openmp

  subroutine statistics_openmp(src, total, peak)
    real(real64), intent(in) :: src(:,:)
    real(real64), intent(out) :: total, peak
    integer :: i, j, nxi, nyi

    nxi = size(src,1)
    nyi = size(src,2)
    total = 0.0_real64
    peak = -huge(0.0_real64)
    !$omp parallel do collapse(2) default(none) &
    !$omp& shared(src, nxi, nyi) private(i, j) &
    !$omp& reduction(+:total) reduction(max:peak) schedule(static)
    do j = 1, nyi
      do i = 1, nxi
        total = total + src(i,j)
        peak = max(peak, src(i,j))
      end do
    end do
    !$omp end parallel do
  end subroutine statistics_openmp

end program heat_parallel

줄별 해설

use iso_fortran_env는 계산과 비교에 사용할 real64를 가져온다. 세 구현이 같은 종류의 실수와 같은 계수를 사용해야 비교의 의미가 분명해진다. nx, ny, nsteps는 검증 문제의 크기와 단계 수를 고정한다.

initial = 0.0_real64은 경계를 포함해 모든 원소를 초기화한다. 중앙값만 64로 바꾼 뒤 세 현재 배열에 복사한다. 각 구현은 자기 배열만 갱신하므로, 먼저 실행한 구현의 결과가 다음 구현의 입력으로 들어가지 않는다.

do step = 1, nsteps는 시간 방향의 순서를 지킨다. 세 갱신 프로시저가 반환한 뒤 다음 배열들을 현재 배열들에 복사한다. 이번 코드에서는 배열 복사의 비용보다 갱신 구조를 명확히 보여 주는 데 초점을 둔다. 실제 성능을 평가할 때는 이 복사와 병렬 영역을 만드는 비용도 포함해야 한다.

maxval(abs(seq - conc))와 대응하는 OpenMP 비교문은 격자의 모든 원소를 검사한다. 두 통계만 같다는 이유로 배열이 같다고 판단하지 않는다. 이어서 순차 합계와 환원 합계를 비교하고, 손으로 확인할 수 있는 합계 60과 최고 온도 20도 검사한다.

error stop은 검증 실패를 정상 출력과 구분한다. 출력 자릿수만 보고 비교하면 서식에 의해 작은 차이가 사라질 수 있으므로, 성공 여부는 출력하기 전에 실제 실수값으로 판정한다. 이번 예제는 유한한 작은 값만 사용한다. 입력 범위가 넓은 실제 계산에서는 비정상 실수값에 대한 검사도 별도로 설계해야 한다.

point_value는 세 구현이 공통으로 사용하는 한 점의 계산식이다. src는 intent(in)이므로 함수가 현재 배열을 바꾸지 않는다. r은 상수이며, 함수가 읽는 외부의 변경 가능한 상태는 없다. 함수 호출은 항상 내부 격자점에서 이루어지므로 네 이웃의 인덱스가 배열 범위 안에 있다.

각 갱신 프로시저의 dst = src는 내부 계산 전에 경계를 포함한 값을 복사한다. 이후 반복 범위는 2부터 각 방향의 끝에서 하나 앞까지다. 따라서 경계는 원래 값을 유지하고 내부만 새 값으로 바뀐다. 현재 배열과 다음 배열은 별개의 실제 인수로 전달해야 한다.

advance_serial은 기준 구현이다. 첫 번째 인덱스 i를 안쪽 반복에 두었다. Fortran 배열에서는 첫 번째 인덱스 방향의 원소가 연속적으로 배치되므로, 이 순회 방식은 배열 접근을 이해하기에도 자연스럽다. 이 기준 구현은 병렬 구현과 같은 함수를 호출하되, 작업 분배만 하지 않는다.

advance_concurrent의 두 제어 변수는 내부 격자점의 조합을 나타낸다. 본문은 그 조합의 한 원소에만 기록한다. 특정 조합이 먼저 실행된다는 가정이 없으며, 현재 배열을 바꾸지 않으므로 이웃 값을 읽는 순서도 결과에 영향을 주지 않는다.

advance_openmp에서는 배열 크기를 nxi와 nyi에 먼저 저장한다. 이어지는 지시문에서 배열과 반복 범위를 shared로, 인덱스를 private로 명시한다. 긴 지시문은 첫 줄 끝의 &와 다음 줄의 !$omp&로 연결한다. 이 표기는 일반 Fortran 실행문과 OpenMP 지시문의 연속 규칙을 함께 고려한 것이다.

statistics_openmp는 갱신과 달리 경계를 포함한 전체 배열을 순회한다. 경계가 0이므로 이번 합계에는 영향을 주지 않지만, 검사 대상을 전체 배열로 통일한다. total과 peak는 shared 목록에 넣지 않고 각각 환원 절에 둔다. 각 반복의 누적 연산은 스레드별 환원 저장 공간을 대상으로 수행된다.

실행 결과

먼저 OpenMP를 활성화하지 않고 컴파일한다. 지시문은 주석이 되며, OpenMP용으로 작성한 프로시저도 일반 반복으로 실행된다. GNU Fortran 16을 사용하는 환경에서 아래 명령으로 빌드한다.

gfortran -std=f2018 -Wall main.f90 -o heat_serial
./heat_serial

OpenMP를 활성화할 때는 컴파일과 링크를 함께 수행하는 명령에 -fopenmp를 추가한다. macOS와 Linux의 셸에서는 환경 변수로 스레드 수를 지정할 수 있다. 이 프로그램은 스레드 수나 실행 시간을 출력하지 않는다.

gfortran -std=f2018 -Wall -fopenmp main.f90 -o heat_openmp
OMP_NUM_THREADS=4 ./heat_openmp

두 실행의 예상 출력은 아래와 같다. 배열은 j가 큰 행부터 출력하여 판의 위쪽에서 아래쪽으로 읽는 형태로 나타낸다. F8.2와 F10.6이 각 실수의 표시 폭과 소수 자릿수를 고정한다.

steps =  2
sequential temperature:
    0.00    0.00    0.00    0.00    0.00
    0.00    2.00    8.00    2.00    0.00
    0.00    8.00   20.00    8.00    0.00
    0.00    2.00    8.00    2.00    0.00
    0.00    0.00    0.00    0.00    0.00
max difference (do concurrent) =   0.000000
max difference (OpenMP)        =   0.000000
sum (sequential) =    60.00
sum (OpenMP)     =    60.00
maximum         =    20.00
verification: PASS

첫 단계에서 중앙은 32, 중앙에 붙은 네 점은 각각 8이 된다. 두 번째 단계에서 중앙은 20, 네 이웃은 각각 8, 네 대각선 위치는 각각 2가 된다. 이 값들을 더하면 60이다. 출력이 이러한 기준값과 일치하는지 확인하면 세 구현이 같은 잘못된 식을 공유하는 경우를 발견하는 데도 도움이 된다.

스레드 수를 1, 2, 4로 바꾸어 같은 출력을 확인할 수 있다. 다만 작은 예제의 성공은 모든 입력과 크기에 대한 증명이 아니다. 격자를 키우거나 경계 조건을 바꿀 때도 배열 비교와 접근 규칙 검토를 이어가야 한다. 특히 이 격자는 작업량이 작아 병렬 실행 준비 비용이 계산 비용보다 클 수 있다.

실무에서 자주 틀리는 것

현재 배열에 바로 다음 온도를 기록한다

아래 반복은 한 점을 갱신한 뒤 다음 점이 바뀐 값을 읽을 수 있다. 순차 실행에서도 앞 장의 명시적 갱신과 다른 계산이 되며, 병렬 실행에서는 읽기와 쓰기가 충돌한다.

! 잘못된 구조
do concurrent (j = 2:ny-1, i = 2:nx-1)
  current(i,j) = point_value(current, i, j)
end do

읽는 배열과 기록하는 배열을 분리하고, 반복이 끝난 뒤에 현재 상태를 바꾼다. 두 배열은 저장 공간도 구분되어야 한다.

! 고친 구조
next = current
do concurrent (j = 2:ny-1, i = 2:nx-1)
  next(i,j) = point_value(current, i, j)
end do
current = next

공유 합계를 여러 반복에서 갱신한다

total을 shared로 둔 채 더하면 두 스레드가 같은 이전 값을 읽고 각각 결과를 기록할 수 있다. 일부 덧셈이 최종값에 반영되지 않는 상황이다. 같은 코드를 do concurrent에 넣는 것도 반복 간 누적 의존성 때문에 적절하지 않다.

! 잘못된 구조
total = 0.0_real64
!$omp parallel do default(none) shared(values, n, total) private(i)
do i = 1, n
  total = total + values(i)
end do
!$omp end parallel do

누적 목적에 맞는 환원을 지정한다. OpenMP를 끄면 이 반복은 같은 초기값에서 시작하는 순차 합산이 된다.

! 고친 구조
total = 0.0_real64
!$omp parallel do default(none) shared(values, n) private(i) reduction(+:total)
do i = 1, n
  total = total + values(i)
end do
!$omp end parallel do

중간 계산값을 shared로 둔다

기록할 배열 원소가 서로 달라도 임시 변수 temp를 공유하면 다른 스레드의 값이 끼어들 수 있다. 다음 코드는 인덱스의 독립성만으로 안전성을 판단할 수 없는 사례다.

! 잘못된 구조
!$omp parallel do default(none) shared(src, dst, n, temp) private(i)
do i = 1, n
  temp = 2.0_real64 * src(i)
  dst(i) = temp
end do
!$omp end parallel do

반복마다 별도로 필요한 임시 변수는 private로 둔다. 또는 완성 코드처럼 순수 함수의 결과를 배열에 바로 대입하여 공유 임시값을 없앨 수 있다.

! 고친 구조
!$omp parallel do default(none) shared(src, dst, n) private(i, temp)
do i = 1, n
  temp = 2.0_real64 * src(i)
  dst(i) = temp
end do
!$omp end parallel do

환원 합계의 비트 일치를 일반 조건으로 삼는다

아래 검사는 일반적인 실수 데이터에서 지나치게 엄격하다. 합산 순서가 바뀌어 생긴 작은 반올림 차이까지 실패로 판정한다. 실수의 직접적인 불일치 비교는 컴파일 옵션에 따라 경고의 대상이 되기도 한다.

! 일반적인 합계 검증에 부적절한 구조
if (total_seq /= total_par) error stop 'sum mismatch'

값의 규모에 맞는 허용 오차를 사용한다. 아래의 atol과 rtol은 문제의 단위, 크기, 요구 정확도에 따라 정해야 한다. 허용 오차를 크게 잡아 데이터 경쟁을 숨기는 방식으로 사용해서는 안 된다.

! 고친 구조
limit = atol + rtol * max(abs(total_seq), abs(total_par))
if (abs(total_seq - total_par) > limit) then
  error stop 'sum mismatch'
end if

배열 갱신과 합계 환원은 따로 검증한다. 각 점의 계산이 같은데 합계의 끝자리만 다른 경우와, 배열 자체가 크게 달라지는 경우는 원인을 다르게 조사해야 한다.

한눈에 보기

순차 반복과 두 병렬화 표현의 역할 및 검증 기준

방식	핵심 역할	주의점	검증
일반 do	기준 결과를 계산한다	제자리 갱신은 계산법을 바꿀 수 있다	손으로 구한 기준값과 비교한다
do concurrent	반복 독립성을 표현한다	여러 코어 실행을 보장하지 않는다	순차 배열과 최대 차이를 구한다
parallel do	스레드에 반복을 나눈다	공유 위치의 쓰기가 겹치지 않아야 한다	스레드 수를 바꾸어 비교한다
reduction	부분 결과를 안전하게 결합한다	실수 합산 순서가 달라질 수 있다	값의 규모에 맞는 오차 기준을 쓴다

작업을 나누기 전에 현재 상태, 다음 상태, 임시값, 누적값을 구분한다. 격자점 계산은 서로 다른 다음 원소에 기록하고, 합계는 환원으로 결합한다. 시간 단계의 완료 지점은 유지한다. 이 세 가지를 코드 구조로 드러내면 병렬화 검토와 결과 검증을 함께 진행하기 쉬워진다.

지시문의 세부 규칙을 확인할 때는 OpenMP 공식 명세를, GNU Fortran의 옵션과 지원 상태를 확인할 때는 GNU Fortran 공식 문서를 참고할 수 있다. 실제 성능은 컴파일러와 실행 환경에 따라 달라지므로, 이 장의 결과 일치 검사를 성능 향상의 근거로 사용하지 않는다.

연습 문제

1. 완성 코드의 nsteps를 1로 바꾼다. 중앙 온도, 전체 합계, 최고 온도의 예상값을 구하고 기준값 검사도 함께 수정한다.
2. OpenMP를 켠 실행에서 스레드 수를 1, 2, 4로 지정한다. 같은 출력이 나와야 하는 이유를 배열 갱신과 환원으로 나누어 설명한다.
3. 모든 갱신 프로시저에서 dst = src를 삭제하면 어떤 문제가 생기는지 설명한다. 고정된 0 경계만 사용하는 조건에서 대체할 수 있는 초기화문도 제시한다.
4. statistics_openmp에 온도 제곱합 squares를 추가한다. 필요한 선언, 초기화, 환원 절, 반복 본문의 변경을 쓰고 두 단계 뒤 예상값을 구한다.

정답과 해설

1. 한 단계 뒤 중앙은 32이고 네 이웃은 각각 8이다. 나머지는 0이므로 전체 합계는 64, 최고 온도는 32다. reference check failed를 검사하는 조건의 60을 64로, 20을 32로 바꾼다. 배열 비교와 순차 합계 대 환원 합계 비교는 그대로 유지한다.

2. 배열 갱신은 같은 현재 배열을 읽고 서로 다른 다음 원소에 기록한다. 각 원소에 적용하는 식도 같으므로 스레드 수에 따라 결과가 달라질 이유가 없다. 환원은 부분 합계의 결합 순서가 달라질 수 있지만, 이번 데이터의 모든 합산이 정확히 표현 가능한 작은 값이므로 합계도 같다. 일반 실수 데이터에서는 환원 합계의 끝자리가 달라질 수 있다.

3. dst는 intent(out)이므로 프로시저 진입 시 이전 값이 정의된 상태라고 기대할 수 없다. 내부 원소만 쓰면 경계가 미정의 상태로 남고, 다음 단계에서 이 경계를 읽게 된다. 이번 조건에서는 반복 전에 dst = 0.0_real64를 실행해 대체할 수 있다. 경계 온도가 0이 아니거나 시간에 따라 바뀌는 경우에는 해당 경계 조건을 직접 설정해야 한다.

4. squares를 real(real64)의 출력 인수로 추가하고 호출부에도 같은 종류의 변수를 둔다. 반복 전에 0으로 초기화하며, 환원 절을 reduction(+:total, squares)로 바꾼다. 기존 최대값 환원은 유지한다. 반복 본문에 아래 누적문을 넣는다.

   squares = squares + src(i,j) * src(i,j)

   예상값은 20의 제곱, 8의 제곱 네 개, 2의 제곱 네 개를 더한 672다. 이 경우에도 제곱합만으로 배열 일치를 판단하지 않고 기존 배열 비교를 함께 사용한다.