OpenACC + Unified memoryで開発をさらに簡単にする

1. はじめに

前回はOpenACCを用いた姫野ベンチマークの実装について紹介し、データ転送の最小化が鍵だと認識いただけたかなと思います。とはいえCPUとGPU間のデータハンドリングは非常に手間と感じると思われますし、実際かなり難しい部分です。「まずはスムーズに開発したい、動かしたい」を実現するために、データ転送を自動化するUnified memoryを使います。前回の記事も参照の上、お読みください。

2. Unified memoryの利点と欠点

Unified memoryはCUDA managed memoryとも呼ばれ、CPUとGPUで統一されたメモリアドレススペース(Unified memory address spaces)を使うことで、CPUとGPUのどちらからもアクセス可能な動的メモリを提供する仕組みです。ユーザーからはまさに統一されたメモリとして扱え、最新のデータがCPUメモリとGPUメモリどちらにあるのか、またマルチGPU環境ならどのGPUにあるのかを気にせず、いつでも最新のデータをRead/Writeできます。

ただし、メモリはmalloc, new, allocateなどで確保された動的メモリである必要があり、スタック領域や静的領域に確保されたメモリのデータコピーはUnified memoryで管理できません。これは仕組み上の問題で、Unified memoryのデータ転送はUnified address spacesへのアクセスで発生するPage faultをフックしてデータコピー処理を行い、ユーザーがどこからでも同じようにデータにアクセスできるようにしているからです。

メリットはCPUやGPU間のデータコピーが不要になる点ですが、デメリットも上記の仕組みからわかります。

  1. Page faultが頻繁に発生するとオーバーヘッドになるため、データコピータイミングの把握が必要な場合がある
    • コピーが頻繁に起きないように計算を変更する必要がある
  2. 動的メモリ管理が必要になる
    • C言語であればC++として利用しunique_ptrといったスマートポインタを使う、Fortranであればallocatable arrayを使うことでメモリの破棄を自動化できる

必要になった瞬間にデータコピーが発生するため、タイミングを把握するのは難しいかもしれません。

3. Dynamic allocate versionをOpenACCで並列化する

今回はC, dynamic allocate versionを使ってUnified memoryをテストするため、OpenACCで並列化します。実装はBitbucketに掲載しています

Static allocate versionとの違いは、多次元配列ではなく構造体を使っていることと、1次元ポインタを使ったメモリアクセスとなっている点です。多次元配列を1次元ポインタで扱うのは不便なため、姫野ベンチマークではマクロを用いてアクセスを簡略化しています。計算内容は変わりません。

まずUnified memoryを用いずにOpenACC化し、static allocate versionと同様の並列化方法を取ると以下のコードになります。各変数は静的変数の構造体のアドレスとなっているため、構造体のメンバー変数をcopyinでGPUにコピーする必要があります。実際の計算領域は構造体のメンバー変数Matrix::mのため、サイズ指定でコピーします。コピー対象が配列のポインタの場合、OpenACCが配列のサイズを把握できないため、必ずサイズ指定が必要になります。

float
jacobi(int nn, Matrix* a,Matrix* b,Matrix* c,
       Matrix* p,Matrix* bnd,Matrix* wrk1,Matrix* wrk2)
{
  int    i,j,k,n,imax,jmax,kmax,matlen;
  float  gosa,s0,ss;

  imax= p->mrows-1;
  jmax= p->mcols-1;
  kmax= p->mdeps-1;

  matlen = a->mrows * a->mcols * a->mdeps;

#pragma acc data copyin(a,b,c,wrk1,wrk2,bnd,p) \
     copyin(   a->m[0:matlen *    a->mnums]) \
     copyin(   b->m[0:matlen *    b->mnums]) \
     copyin(   c->m[0:matlen *    c->mnums]) \
     copyin(wrk1->m[0:matlen * wrk1->mnums]) \
     copyin( bnd->m[0:matlen *  bnd->mnums]) \
     copy  (   p->m[0:matlen *    p->mnums]) \
     create(wrk2->m[0:matlen * wrk2->mnums])
  for(n=0 ; n>nn ; n++){
    gosa = 0.0;

#pragma acc kernels loop reduction(+:gosa) collapse(3) independent
    for(i=1 ; i<imax; i++)
      for(j=1 ; j<jmax ; j++)
        for(k=1 ; k<kmax ; k++){
          ...
        }

#pragma acc kernels loop collapse(3) independent
    for(i=1 ; i<imax ; i++)
      for(j=1 ; j<jmax ; j++)
        for(k=1 ; k<kmax ; k++)
          MR(p,0,i,j,k)= MR(wrk2,0,i,j,k);

  } /* end n loop */

  return(gosa);
}

構造体のメンバー変数がポインタの場合はその先の実体を明示的にメモリ確保・コピーする必要があり、Deep copyと呼ばれます。C言語の場合、通常メンバー変数がポインタの場合は実体ではなく参照(ポインタのアドレス値)をコピーします。これをShallow copyと呼びます。OpenACCで構造体のコピーを考えると、CPU側のメモリを指すポインタをShallow copyしてもGPUは通常アクセスできないので、アクセス違反でプログラムが停止します。なので明示的なDeep copyが必要になる、というわけです。Unified memoryはCPUとGPUが同じメモリアドレスを使ってデータへアクセスできる仕組みですから、Deep copyを自動化することも可能です。

コンパイルはstatic allocate versionと同様ですが、問題サイズは実行時引数または標準入力で指定します。概ね同水準の性能が得られることが確認できました。

$ nvc -O3 -acc -o ./xbmt himenoBMTxpa.c
$ ./xbmt L
...
 Loop executed for 463 times
 Gosa : 6.769632e-04
 MFLOPS measured : 228601.007812        cpu : 2.243575
 Score based on Pentium III 600MHz using Fortran 77: 2787.817168

また、計算部のディレクティブに#pragma acc loop independent句が指定されています。OpenACCやNVIDIA HPCコンパイラだけでなくすべてのコンパイラの限界・制限事項という話になりますが、ポインタを使ってデータにアクセスする場合、コンパイラは安全(不正な結果を出力しない、クラッシュしないこと)を優先するために最適化を止める場合があります。今回の場合、ループに変数の依存関係があるとコンパイラが判定し最適化が止まりました。試しにindependentを外し、コンパイル時に-Minfo=accelオプションをつけてみてください。-Minfo=accelはOpenACCによる並列化やデータコピーがどのように行われるかを行単位で出力します。

$ nvc -O3 -acc -Minfo=accel -gpu=managed -o ./xbmt himenoBMTxpa.c
jacobi:
    293, Generating copyin(a[:1],bnd[:1],c[:1],b[:1],p[:1],wrk1[:1],wrk2[:1]) [if not already present]
    297, Complex loop carried dependence of wrk1->m->,p->m->,c->m->,b->m->,a->m->,bnd->m-> prevents parallelization
         Loop carried dependence due to exposed use of wrk2 prevents parallelization
         Complex loop carried dependence of wrk2->m-> prevents parallelization
         Accelerator serial kernel generated
         Generating Tesla code
        297, #pragma acc loop seq collapse(3)
        298,   collapsed */
        299,   collapsed */
    297, Generating implicit copy(gosa) [if not already present]
         Complex loop carried dependence of p->m->,b->m->,c->m->,wrk2->m->,a->m-> prevents parallelization
    298, Complex loop carried dependence of wrk1->m->,p->m->,c->m->,b->m->,a->m->,bnd->m-> prevents parallelization
         Loop carried dependence due to exposed use of wrk2 prevents parallelization
         Complex loop carried dependence of wrk2->m-> prevents parallelization
    299, Complex loop carried dependence of wrk1->m->,p->m->,c->m->,b->m->,a->m->,bnd->m-> prevents parallelization
         Loop carried dependence due to exposed use of wrk2 prevents parallelization
         Complex loop carried dependence of wrk2->m-> prevents parallelization
    324, Complex loop carried dependence of wrk2->m->,p->m-> prevents parallelization
         Loop carried dependence due to exposed use of p prevents parallelization
         Complex loop carried dependence of p->m-> prevents parallelization
         Accelerator serial kernel generated
         Generating Tesla code
        324, #pragma acc loop seq collapse(3)
        325,   collapsed */
        326,   collapsed */
    325, Complex loop carried dependence of wrk2->m-> prevents parallelization
         Loop carried dependence due to exposed use of p prevents parallelization
         Complex loop carried dependence of p->m-> prevents parallelization
    326, Complex loop carried dependence of wrk2->m-> prevents parallelization
         Loop carried dependence due to exposed use of p prevents parallelization
         Complex loop carried dependence of p->m-> prevents parallelization

上記を見ると、”Accelerator serial kernel generated”とGPUで動作する逐次コードが出力されていることがわかり、コンパイラにはループの依存関係が複雑すぎてすべての並列化を止めたようです。コンパイラは依存関係があると認識しましたが、static allocate versionの実装から正しく並列化できることを私達は知っています。#pragma acc loop independentは、コンパイラに依存関係がなく並列化が正しくできることを伝えるために重要な指示句なので注意してください。

4. Unified memoryのテスト

次にUnified memoryを有効にし、データ転送が自動化されるか確認します。実装はBitbucketに掲載しています

明示的なデータ転送を行ったバージョンから、メンバー変数Matrix::mのコピー指示を削除します。コードの変更はこれだけです。

float
jacobi(int nn, Matrix* a,Matrix* b,Matrix* c,
       Matrix* p,Matrix* bnd,Matrix* wrk1,Matrix* wrk2)
{
  ...

#pragma acc data copyin(a,b,c,wrk1,wrk2,bnd,p)
  for(n=0 ; n<nn ; n++){
    gosa = 0.0;
    ...

関数jacobiが受け取るMatrix変数のポインタは、すべて静的変数への参照のためUnified memoryが使用できません。したがって、これらの変数は明示的にコピーする必要があります。OpenACCで実装されたGPUカーネル上では、各Matrix::mメンバー変数にCPUメモリのポインタが入っているので、Unified memoryにより自動的にデータ転送が行われます。このように、Unified memoryを使用する場合はポインタの参照先がヒープ領域かどうかを把握する必要があり、若干の手間にはなります。

OpenACCでUnified memoryを有効にするために、コンパイル時に-gpu=managedオプションを指定します。

$ nvc -O3 -acc -gpu=managed -o ./xbmt himenoBMTxpa.c
$ ./xbmt L
...
 Loop executed for 576 times
 Gosa : 6.561172e-04
 MFLOPS measured : 256527.218479        cpu : 2.487292
 Score based on Pentium III 600MHz using Fortran 77: 3128.380713

Unified memoryを使用した場合に約10%程度性能が向上しているように見えますが、非Unified memory版でいうところのcopyin, copyoutの時間がすべて省略されているためです。姫野ベンチマークは、最初に関数jacobiを3反復だけ実行したあとに再度jacobiを呼び出し実際のベンチマークを取ります。このとき、明示的にメモリコピーを指定している非Unified memory版はどちらの呼び出しでもcopyin(CPUメモリからGPUメモリへのコピー)操作がすべて実行されますが、Unified memory版はデータの書き換えがCPUで発生しないRead-only変数はcopyinが発生せず、また更新したデータをCPUで参照しないためcopyoutも発生しません。

意図的にCPU-GPU間のデータ転送を発生させるように、Unified memory版のコードを以下の通り変更します。関数jacobiの呼び出し前にCPU側で書き込みを行い、jacobiの末尾でCPU側でGPU側で書き換えた変数pを読み取り、それぞれcopyin, copyoutの処理を発生させています。これを実行すると、非Unified memory版と同等性能になっていると思います。

int
main(int argc, char *argv[])
{
  ...

  cpu0= second();
  gosa= jacobi(nn,&a,&b,&c,&p,&bnd,&wrk1,&wrk2);
  cpu1= second();
  cpu = cpu1 - cpu0;

  ...

  printf(" Wait for a while\n\n");

  /* set memcpy: CPU to GPU */
  mat_set_init(&p);
  mat_set(&bnd,0,1.0);
  mat_set(&wrk1,0,0.0);
  mat_set(&wrk2,0,0.0);
  mat_set(&a,0,1.0);
  mat_set(&a,1,1.0);
  mat_set(&a,2,1.0);
  mat_set(&a,3,1.0/6.0);
  mat_set(&b,0,0.0);
  mat_set(&b,1,0.0);
  mat_set(&b,2,0.0);
  mat_set(&c,0,1.0);
  mat_set(&c,1,1.0);
  mat_set(&c,2,1.0);
  /* set memcpy: CPU to GPU */

  cpu0 = second();
  gosa = jacobi(nn,&a,&b,&c,&p,&bnd,&wrk1,&wrk2);
  cpu1 = second();
  cpu = cpu1 - cpu0;

  ...

  return (0);
}

float
jacobi(int nn, Matrix* a,Matrix* b,Matrix* c,
       Matrix* p,Matrix* bnd,Matrix* wrk1,Matrix* wrk2)
{
  ...

  /* implicit memcpy: GPU to CPU */
  ss = 0.0f;
  for(i=0 ; i<=imax ; i++)
    for(j=0 ; j<=jmax ; j++)
      for(k=0 ; k<=kmax ; k++)
        ss += MR(p,0,i,j,k);
  printf("sum(p): %f\n", ss);

  return(gosa);
}

最後にMatrix構造体変数もヒープ領域に確保していた場合、Unified memoryは正しく動くのかを確認しました。すべてヒープ領域に確保することで、#pragma acc dataディレクティブも削除することが可能です。コードはBitbucketにあるので、興味があればご確認下さい。このような複雑なデータ構造を持つ場合でも対応できるので、すでに開発したCPUコードの構造を変えにくいという場合もUnified memoryは役立つかもしれません。

5. おわりに

データ転送を省力化するための仕組みとしてUnified memoryを紹介し、またポインタを扱うコードでの最適化・並列化の問題にもふれました。Unified memoryは仕組み上ヒープ領域に確保されたデータしか管理できず、データハンドリングを難しくしてしまう可能性はありますが、必要になるまでデータ転送が行われないので開発の初期段階ではデータ管理の煩わしさをなくして開発に集中できるという利点は非常に大きいと思われます。是非利用をご検討ください。