前回，LAPACK の関数を cuSOLVER に置き換えて GPU オフロードを実現しました。そしてベンチマークを取って，cuSOLVER が大規模計算向きのツールであることを確認しました。
Fortranで有限要素法解析 04 : cuSOLVER を使う 〜 其の2 高速化 〜
そうであるならば，2次元問題の方が GPU 使用の優位性がより顕著に現れそうです（実際，1次元問題でそんなに多数のメッシュを取るケースはないだろうし，仮に CPU を使って解いたとしてもそのコストは高が知れており 1次元問題でGPUを使うメリットは考えにくい）。
そこで，2次元問題に対して同様の GPU 移植と計算速度の比較を試みてみます。
LSVQR
1次元問題でも使用した cuSOLVER の LSVQR（疎行列用の関数）を使用します。2次元問題では係数行列がブロック三重対角になることが異なりますが，それ以外の基本的なアルゴリズム構造は同じです。従って，1次元問題のときに省略した説明を中心に書いていきたいと思います。
メインプログラムは次のようになります。
温度分布を初期化（initialize）した後， それらを GPU にコピー（fem2d_to_device）してから時間ループのイテレーション（evolve）を実行します。最終結果はプログラムの最後で CPU に戻します（fem2d_to_host）。また，必要に応じて，イテレーションの最中に経過データを出力します。
ここで，temperature に u と u_work の２種類があるのは，時間方向に前進差分スキームを使っているからで，一方（既知時刻のデータ）を代数方程式の右辺に，他方（未知時刻のデータ）を左辺に使用して温度分布を更新していきます。
program heat2d
...
call initialize(nx, ny, temperature)
call fem2d_to_device(nx*ny, nx, ny, temperature%u, temperature%u_work)
do iter = 1, nsteps
call evolve(temperature, kappa, dt)
if (mod(iter,out_interval) == 0) then
!$acc update host(temperature%u)
call plot(temperature%u)
end if
end do
call fem2d_to_host(nx*ny, nx, ny, temperature%u, temperature%u_work)
end program
ホスト - デバイス間のデータコピーのルーチン（fem2d_to_device, fem2d_to_host）はそれぞれ次のようになります。
subroutine fem2d_to_device(n, nx, ny, u, u_work)
integer, intent(in) :: n, nx, ny
real(real64), intent(inout) :: u(nx, ny), u_work(nx, ny)
associate(bvec=>u_work, xvec=>u)
!$acc enter data copyin(xvec) &
!$acc create(csrValA, csrRowPtrA, csrColIndA, bvec)
end associate
end subroutine
subroutine fem2d_to_host(n, nx, ny, u, u_work)
integer, intent(in) :: n, nx, ny
real(real64), intent(inout) :: u(nx, ny), u_work(nx, ny)
associate(bvec=>u_work, xvec=>u)
!$acc exit data copyout(xvec) &
!$acc delete(csrValA, csrRowPtrA, csrColIndA, bvec)
end associate
end subroutine fem2d_to_host
代数方程式を解くプロセスは 1 次元のときと同じです。係数行列の作成（assemble_matrix）, 右辺ベクトルの作成（assemble_rhs）, 求解（solve_cusolver_sp）の順に行います。
計算結果は次のようになりました（中心付近の熱が周辺に拡散していっています）。
ベンチマーク
ステップ数 100 までの計算時間を，LAPACK を使用した CPU コードと cuSOLVER を使用した GPU 移植コードで比較しました（単位は秒）。残念ながら最大規模の計算は CPU で実行できませんでしたが（メモリ不足），計算規模が大きいとき，GPU の方に優位性が移りそうだと見て取れます。
| ソルバ | デバイス | n=64$\times$64 | n=128$\times$128 | n=192$\times$192 | n=256$\times$256 |
|---|---|---|---|---|---|
| GBTRF/GBTRS（帯行列 LAPACK） | CPU | 3.9 | 58.0 | 269.6 | — |
| LSVQR（疎行列 cuSOLVER） | GPU w/ OpenACC | 11.1 | 90.3 | 319.8 | 827.8 |
| LSVQR（疎行列 cuSOLVER） | GPU w/ OpenMP | 11.1 | 89.5 | 321.6 | 830.9 |