OpenACC ディレクティブによるプログラミング
12 - 1 章 OpenACC 上での CUDA Library の利用
- cuSOLVERライブラリの紹介(Dense LAPACK, Sparse LAPACK)
- cuSolverDN: エルミート正定値行列のコレスキー分解 (LAPACK Zpotrf相当）
- cuSolverDN: 対称正定値行列のコレスキー分解を使ったソルバー (LAPACK Dpotrs相当)
- cuSolverDN: LU分解を使った線形システムソルバー (LAPACK Dpotrs相当)
- cuSolverDN: 対称行列の固有値ソルバー (LAPACK Dsyevd相当)
- cuSolverDN: QR分解による密行列線形ソルバー (LAPACK Dgeqrf 相当)
- cuSolverDN: QR分解による直交化 (LAPACK Dgeqrf,Dorgqr 相当）
- cuSolverDN: 特異ベクトル計算を含めた特異値分解(by QR) (LAPACK Dgesvd 相当）
- cuSolverDN: 特異ベクトル計算を含めた特異値分解(By Jacobi) (LAPACK Dgesvd 相当）
- cuSolverSP: 並列ダイレクト・スパース・ソルバー(CSR QR分解)
- cuFFT example with OpenACC Fortran (FFT)
- cuRAND example with OpenACC Fortran (乱数発生)
- cuSPARSE example from OpenACC Host code
1. cuSOLVERライブラリの紹介 (Dense LAPACK, Sparse LAPACK)
cuSOLVER ライブラリは、cuBLAS および cuSPARSE ライブラリに基づく高水準パッケージである。それぞれのライブラリは、単独で、または他のツールキットライブラリと連携して使用できる。cuSolver の目的は、密行列、スパース行列に対する最小二乗ソルバ、固有値ソルバのための一般的な行列分解や三角法のような便利な CUDA 版 LAPACK の機能を提供することである。さらに、cuSOLVER は、スパースパターンを共有した行列のシーケンスを解くのに役立つ新しいリファクタリングライブラリも提供する。
PGI 17.10 現在、cuSOLVER ライブラリへの Fortran Interface (cusolverDn.mod) を提供しているライセンス製品は Linux 用のもののみとなります。Windows 版においてはこの時点では Fortran Interface Module が提供されていないため、ご注意下さい。なお、Fortran interface module を作成すると、Windows でも使用できます。
cuSOLVER の最初の部分は cuSolverDN と呼ばれ、密行列を対象とした分解や、LU、QR、SVD、LDLT などの解法ルーチン、行列やベクトル置換などの有用なユーティリティを提供している。いわゆる LAPACK ルーチンと同様な機能を提供する。
次に、cuSolverSP は、スパース QR 分解に基づく新しいスパースルーチンのセットを提供する。すべての行列分解が並列性に優れたスパース性のパターンを持つわけではないため、cuSolverSP ライブラリはシーケンシャルのような行列を扱う CPU パスも提供する。高い並列性を有するマトリックスの場合、GPU パスはより高いパフォーマンスを提供する。ライブラリは、C および C++ から呼び出されるように設計されているが、OpenACC + PGI Fortran Module を介して Fortran からも簡単に利用できる。
三つ目は、cuSolverRF です。これは、係数のみが変更され、スパース性パターンは同じままである行列のシーケンスを解く際に、非常に優れたパフォーマンスを提供することができるスパース行列用 refactorization パッケージである。（このページでは説明していません）
cuSolverDN ライブラリ(Dense LAPACK)
cuSolverDN ライブラリは、密行列の線形システムを解くように設計されている。
ここで、係数行列A∈Rnxn、右辺ベクトルb∈Rn、解ベクトルx∈Rn
cuSolverDNライブラリは、非対称である可能性のある一般的な行列 A を処理するために、QR 分解と LU に部分ピボッティングを提供し、対称/エルミート行列に対してコレスキー分解ルーチンも提供されている。また、対称不定行列の場合、Bunch-Kaufman（LDL）分解を提供している。cuSolverDN ライブラリは、有用な二重対角化ルーチンと特異値分解（SVD）も提供している。
cuSolverDN ライブラリは、LAPACK の計算集約的で一般的なルーチンを対象とし、LAPACK と互換性のある API を提供していますので、既存のLAPACKを使ったプログラムの移行が容易である。ユーザーは cuSolverDN を使用して、これらの時間の掛かるルーチンを加速し、既存のコードを大幅に変更することなく LAPACK 利用の互換性を保つことができる。
cuSolverSP ライブラリ(Sparse LAPACK)
cuSolverSPライブラリは、主にスパースな線形システムと最小二乗問題を解くように設計されている。
ここで、係数行列A∈Rnxn、右辺ベクトルb∈Rn、解ベクトルx∈Rn。線形システムでは、m = n が必要である。
コア・アルゴリズムは、スパース QR 分解に基づいている。行列 A は CSR 形式で入力します。行列 A が対称/エルミート行列である場合、ユーザーは完全な行列（full matrix) を提供する必要がある。つまり、欠落している下部または上部を埋める必要がある。もし、行列 A が対称正定値であり、ユーザーが Ax = b を解く必要があるだけであれば、コレスキー分解は機能するが、ユーザーは行列 A の下三角部分を提供することで代替できる。
線形および最小二乗ソルバーの上に、cuSolverSP ライブラリは、 shift-inverse power method に基づく簡単な固有値ソルバーと、複素平面内のボックスに含まれる固有値の数を数える関数を提供する。
2. cuSolverDN: エルミート正定値行列のコレスキー分解 (LAPACK zpotrf相当）
cuSOLVER ライブラリは cuBLAS および cuSPARSE ライブラリをベースにしており、密行列、スパース最小二乗ソルバおよび固有値ソルバの行列分解および三角解法ルーチンをサポートし、共分散行列を持つ行列のシーケンスを解くのに役立つリファクタリングライブラリを提供する。 PGI が提供するインターフェイスモジュール（Fortran MODULE)と、PGI 17.7 以降にバンドルされている cuSOLVER ライブラリの PGI コンパイル済みバージョンを使用して、CUDA Fortran および OpenACC Fortran から最適化された cuSolverDN type ルーチン群（cuSolverSP ならびに cuSolverRF type は今後のリリースで対応予定）を呼び出すことができるようになった。この同じ cuSolver ライブラリは、PGI コンパイラを使用して構築されており、また、PGI OpenMP ランタイムと互換性があるため、PGI OpenACC C/C++ からも呼び出すことができる。
以下の例は、cuSolverDN: dense LAPACK Function を OpenACC Fortran/C++ ならびに CUDA Fortran から利用する場合の一例である。
配列、変数のデバイスメモリのマネージメントをどちら側のプログラミングモデルで行うかによって、プログラミングの対応は異なるが、ここでは、 CUDA Fortran 側でデバイス配列の宣言を行った場合の方法、OpenACC 側で配列の宣言を行った場合を以下に例示する。
なお、以下の例は、PGI 17.10 バージョンを使用しての結果である。使用しているバージョンを確かめたい場合は、以下のように -V コマンド・オプションを指定する。
$ pgfortran -V pgfortran 17.10-0 64-bit target on x86-64 Linux -tp haswell PGI Compilers and Tools Copyright (c) 2017, NVIDIA CORPORATION. All rights reserved.
Dense Linear Solver - cusolverDnZpotrf() with CUDA Fortran
エルミートの正定値行列のコレスキー分解の一例を示す。LAPACK の ZPOTRF ルーチン相当である。A は n×n のエルミート行列であり、下位または上位の部分のみが意味を持つ。以下の例は CUDA Fortran から cuSOLVER を使用したプログラム例である。PGI Fortran の MODULE インタフェースとして、cublas_v2 と cusolverDn を USE 文で指定することが必要である。入力パラメータ CUBLAS_FILL_MODE_LOWER の場合、 Aの下三角部分のみが処理され、下三角のコレスキー係数 L で置き換えられます。(A = L * L H) 入力パラメータ Workspace で指し示す作業領域を提供する必要がある。入力パラメータ Lwork は作業領域のサイズであり、 potrf_bufferSize（）によって返される。なお、密行列 A は、column-major order で入力する必要がある。
CUDA Fortran + cuSOLVER(cusolverDnZpotrf)
! Copyright (c) 2017, NVIDIA CORPORATION. All rights reserved. ! ! ! Permission is hereby granted, free of charge, to any person obtaining a ! copy of this software and associated documentation files (the "Software"), ! to deal in the Software without restriction, including without limitation ! the rights to use, copy, modify, merge, publish, distribute, sublicense, ! and/or sell copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER ! DEALINGS IN THE SOFTWARE. ! program main use cublas_v2 use cusolverDn use cudafor implicit none integer, parameter :: n=3 complex(8) :: a(n,n) complex(8), device :: a_d(n,n) complex(8), device, allocatable :: workspace_d(:) integer, device :: devInfo_d integer :: istat, Lwork type(cusolverDnHandle) :: h a(1,1) = 25.0; a(1,2) = 15.0; a(1,3) = -5.0 a(2,1) = a(1,2); a(2,2) = 18.0; a(2,3) = 0.0 a(3,1) = a(1,3); a(3,2) = a(2,3); a(3,3) = 11.0 a_d = a ! handle 作成 istat = cusolverDnCreate(h) if (istat /= CUSOLVER_STATUS_SUCCESS) & write(*,*) 'handle creation failed' ! Working buffer サイズ計算 istat = cusolverDnZpotrf_bufferSize(h, & CUBLAS_FILL_MODE_LOWER, n, a_d, n, Lwork) if (istat /= CUSOLVER_STATUS_SUCCESS) & write(*,*) 'cusolverDnZpotrf_buffersize failed' allocate(workspace_d(Lwork)) ! コレスキー分解 istat = cusolverDnZpotrf(h, CUBLAS_FILL_MODE_LOWER, & n, a_d, n, workspace_d, Lwork, devInfo_d) if (istat /= CUSOLVER_STATUS_SUCCESS) & write(*,*) 'cusolverDnZpotrf failed' istat = devInfo_d if (istat /= 0) write(*,*) 'Cholesky factorization failed' istat = cusolverDnDestroy(h) if (istat /= CUSOLVER_STATUS_SUCCESS) & write(*,*) 'handle destruction failed' a = a_d write(*,"(3(f0.0,SP,f0.0,'i',2x))") a(1,:) write(*,"(3(f0.0,SP,f0.0,'i',2x))") a(2,:) write(*,"(3(f0.0,SP,f0.0,'i',2x))") a(3,:) end program main
コンパイルオプションとして、-acc -Mcudalib=cusolver,cublas、-Mcuda オプションを付けてコンパイルリンクする必要がある。
コンパイル＆実行
$ pgfortran -fast -Minfo -Mcuda=cc60,cc35,cuda8.0 -Mcudalib=cusolver,cublas testDn.cuf
main:
60, Loop unrolled 3 times (completely unrolled)
61, Loop unrolled 3 times (completely unrolled)
62, Loop unrolled 3 times (completely unrolled)
$ ./a.out
5.+0.i +15.+0.i -5.+0.i
3.+0.i +3.+0.i +0.+0.i
-1.+0.i +1.+0.i +3.+0.i
Dense Linear Solver - cusolverDnZpotrf() with OpenACC Fortran
上記の CUDA Fortran プログラムを OpenACC を使って書き換える。PGI Fortran の MODULE インタフェースとして、cublas_v2 と cusolverDn を USE 文で指定することが必要である。デバイス側のメモリ割付を気にする必要なく、OpenACCのディレクティブのみで、CUDA 数学ライブラリが利用できる。既存の Fortran プログラムをポーティングすることも簡単にできるということである。
OpenACC のディレクティブで留意すべきことは、CUDA Fortran or CUDA C で記述された（ライブラリ）ルーチンをコールする際に渡す実引数が「デバイスポインタ」である場合、これをコンパイルに伝えるために、host_data use_device() を使う。このディレクティブの使い方さえ誤りがなければ、既存のプログラムベースのものを OpenACC を使ってポーティングすることは大きな負担ではないはずだ。
OpenACC Fortran + cuSOLVER(cusolverDnZpotrf)
! Copyright (c) 2017, NVIDIA CORPORATION. All rights reserved. ! ! ! Permission is hereby granted, free of charge, to any person obtaining a ! copy of this software and associated documentation files (the "Software"), ! to deal in the Software without restriction, including without limitation ! the rights to use, copy, modify, merge, publish, distribute, sublicense, ! and/or sell copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER ! DEALINGS IN THE SOFTWARE. ! program main use cublas_v2 use cusolverDn use cudafor implicit none integer, parameter :: n=3 complex(8) :: a(n,n) !!! complex(8), device :: a_d(n,n) ! デバイス変数の宣言必要なし !!! complex(8), device, allocatable :: workspace_d(:) !!! integer, device :: devInfo_d complex(8), allocatable :: workspace_d(:) integer :: devInfo_d integer :: istat, Lwork type(cusolverDnHandle) :: h a(1,1) = 25.0 ; a(1,2) = 15.0; a(1,3) = -5.0 a(2,1) = a(1,2); a(2,2) = 18.0; a(2,3) = 0.0 a(3,1) = a(1,3); a(3,2) = a(2,3); a(3,3) = 11.0 !!! a_d = a istat = cusolverDnCreate(h) if (istat /= CUSOLVER_STATUS_SUCCESS) & write(*,*) 'handle creation failed' !$acc data copy(a) create(workspace_d) !$acc host_data use_device(a) istat = cusolverDnZpotrf_bufferSize(h, & CUBLAS_FILL_MODE_LOWER, n, a, n, Lwork) !$acc end host_data if (istat /= CUSOLVER_STATUS_SUCCESS) & write(*,*) 'cusolverDnZpotrf_buffersize failed' print *, "working space (bytes) =", Lwork allocate(workspace_d(Lwork)) !$acc enter data copyin(workspace_d) copyin(devInfo_d) !$acc host_data use_device(a, workspace_d, devInfo_d) istat = cusolverDnZpotrf(h, CUBLAS_FILL_MODE_LOWER, & n, a, n, workspace_d, Lwork, devInfo_d) !$acc end host_data if (istat /= CUSOLVER_STATUS_SUCCESS) & write(*,*) 'cusolverDnZpotrf failed' !$acc exit data copyout(devInfo_d) !$acc exit data delete(workspace_d) !$acc end data istat = devInfo_d if (istat /= 0) write(*,*) 'Cholesky factorization failed' istat = cusolverDnDestroy(h) if (istat /= CUSOLVER_STATUS_SUCCESS) & write(*,*) 'handle destruction failed' !!! a = a_d write(*,"(3(f0.0,SP,f0.0,'i',2x))") a(1,:) write(*,"(3(f0.0,SP,f0.0,'i',2x))") a(2,:) write(*,"(3(f0.0,SP,f0.0,'i',2x))") a(3,:) end program main
コンパイルオプションとして、-acc、-Mcudalib=cusolver,cublas、-Mcuda オプションを付けてコンパイルリンクする必要がある。OpenACC の -acc オプションの他に、-Mcuda オプションは、CUDA関係のシステムライブラリ・リストをリンカーに伝えるために必要である。
コンパイル＆実行
$ pgf90 -fast -Minfo -acc -ta=tesla,cc60,cc35,cuda8.0 -Mcudalib=cusolver,cublas openacc.f90 -Mcuda
main:
48, Generating copy(a(:,:))
Generating create(workspace_d(:))
60, Generating enter data copyin(devinfo_d,workspace_d(:))
70, Generating exit data copyout(devinfo_d)
71, Generating exit data delete(workspace_d(:))
84, Loop unrolled 3 times (completely unrolled)
85, Loop unrolled 3 times (completely unrolled)
86, Loop unrolled 3 times (completely unrolled)
$ ./a.out
working space (bytes) = 1
5.+0.i +15.+0.i -5.+0.i
3.+0.i +3.+0.i +0.+0.i
-1.+0.i +1.+0.i +3.+0.i
Dense Linear Solver - cusolverDnZpotrf with OpenACC C++
上記 Fortran プログラムを NVIDIA nvcc コンパイラではなく、PGI C++ OpenACC を使用して書き直す。C/C++ プログラムにおいても、OpenACC と NVIDIA CUDA ライブラリを使用するプログラムを作成できる。なお、以下の例は NVIDIA cuSOLVE ライブラリを C++ プログラムから参照する際に、複素数型は、CUDA 複素数型（例えば、cuDoubleComplex等）を使用して宣言する必要がある。C++の std:complex テンプレートのcomplex
OpenACC + cusolverDnDsyevd (DnZpotrf.cpp)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved.
!
! Permission is hereby granted, free of charge, to any person obtaining a
! copy of this software and associated documentation files (the "Software"),
! to deal in the Software without restriction, including without limitation
! the rights to use, copy, modify, merge, publish, distribute, sublicense,
! and/or sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.
/*
* How to compile (assume cuda is installed at /usr/local/cuda-9.0/)
* pgc++ -I/usr/local/cuda-9.0/include DnZpotrf.cpp -acc -ta=tesla,cc60,cc35 -Mcudalib=cusolver -Minfo=accel -O2
*
*
*/
#include <stdio.h>
#include <stdlib.h>
#include <assert.h>
#include <cusolverDn.h>
void printD_CMatrix(int m, int n, const cuDoubleComplex *A, int lda, const char* name)
{
for(int row = 0 ; row < m ; row++){
for(int col = 0 ; col < n ; col++){
cuDoubleComplex Areg = A[row + col*lda];
printf("%s(%d,%d) = %f %f\n", name, row+1, col+1, cuCreal(Areg), cuCimag(Areg));
}
}
}
int main(int argc, char*argv[])
{
cusolverDnHandle_t cusolverH = NULL;
cusolverStatus_t cusolver_status = CUSOLVER_STATUS_SUCCESS;
const int m = 3;
const int lda = m;
/* | 25 15 -5 |
* A = | 15 18 0 |
* | -5 0 11 |
*
*/
cuDoubleComplex *A; // input
int *devInfo = NULL;
cuDoubleComplex *work = NULL;
int lwork = 0;
A = (cuDoubleComplex *) malloc (sizeof(cuDoubleComplex) * lda * m);
devInfo = (int *) malloc (sizeof(int));
// step 1: create cusolver/cublas handle
cusolver_status = cusolverDnCreate(&cusolverH);
// step 2: set A
A[0] = make_cuDoubleComplex(25.0,0.0);
A[3] = make_cuDoubleComplex(15.0,0.0);
A[6] = make_cuDoubleComplex(-5.0,0.0);
A[1] = A[3];
A[4] = make_cuDoubleComplex(18.0,0.0);
A[7] = make_cuDoubleComplex(0.0,0.0);
A[2] = A[6];
A[5] = A[7];
A[8] = make_cuDoubleComplex(11.0,0.0);
printD_CMatrix(m, m, A, lda, "A");
printf("=====\n");
// step 3: query working space of syevd
cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER;
#pragma acc data copy(A[0:lda*m], devInfo[0:1])
{
#pragma acc host_data use_device(A)
{
cusolver_status = cusolverDnZpotrf_bufferSize(
cusolverH,
uplo,
m,
A,
lda,
&lwork);
}
assert (cusolver_status == CUSOLVER_STATUS_SUCCESS);
work = (cuDoubleComplex *) malloc(sizeof(cuDoubleComplex)*lwork);
printf("Lwork = %d\n",lwork);
#pragma acc enter data create(work[0:lwork])
// step 4: compute spectrum
#pragma acc host_data use_device(A, work, devInfo)
{
cusolver_status = cusolverDnZpotrf(
cusolverH,
uplo,
m,
A,
lda,
work,
lwork,
devInfo);
}
#pragma acc wait
#pragma acc exit data delete(work)
}
printf("DevInfo = %d\n",devInfo[0]);
printf("A = (base 1)\n");
printD_CMatrix(m, m, A, lda, "A_result");
printf("=====\n");
// free resources
if ( A ) free(A);
if (devInfo) free(devInfo);
if (work ) free(work);
if (cusolverH) cusolverDnDestroy(cusolverH);
return 0;
}
PGI C/C++ OpenACC コンパイラで cuSOLVER を使用する場合は、NVIDIA CUDA Toolkit にバンドルされている C/C++ 用 include ファイルを使用してコンパイルする(PGI コンパイラ内には、C/C++ 用の cuSOLVER/cuBLAS ライブラリ等の include ファイルはバンドルされていない） 。ここでの例は CUDA 9.0 が /usr/local/cuda-9.0 配下に実装されているものとして記述する。コンパイルオプションとして、 -I/usr/local/cuda-9.0/include -acc -ta=tesla -Mcudalib=cusolver オプションを付けてコンパイルリンクする必要がある。
(CUDA 8.0 の場合）
pgc++ -I/usr/local/cuda-8.0/include DnZpotrf.cpp -acc -ta=tesla,cc60,cc35 -Mcudalib=cusolver -Minfo=accel -O2
(CUDA 9.0 の場合）
pgc++ -I/usr/local/cuda-9.0/include DnZpotrf.cpp -acc -ta=tesla,cc60,cc35 -Mcudalib=cusolver -Minfo=accel -O2
コンパイル＆実行
[kato@photon32 C]$ pgc++ -I/usr/local/cuda-9.0/include DnZpotrf.cpp -acc -ta=tesla,cc60,cc35
-Mcudalib=cusolver -Minfo=accel -O2
main:
64, Generating copy(A[:9],devInfo[:1])
84, Generating enter data create(work[:lwork])
97, Generating exit data delete(work[:1])
[kato@photon32 C]$ a.out
A(1,1) = 25.000000 0.000000
A(1,2) = 15.000000 0.000000
A(1,3) = -5.000000 0.000000
A(2,1) = 15.000000 0.000000
A(2,2) = 18.000000 0.000000
A(2,3) = 0.000000 0.000000
A(3,1) = -5.000000 0.000000
A(3,2) = 0.000000 0.000000
A(3,3) = 11.000000 0.000000
=====
Lwork = 1
DevInfo = 0
A = (base 1)
A(1,1) = 5.000000 0.000000
A(1,2) = 15.000000 0.000000
A(1,3) = -5.000000 0.000000
A(2,1) = 3.000000 0.000000
A(2,2) = 3.000000 0.000000
A(2,3) = 0.000000 0.000000
A(3,1) = -1.000000 0.000000
A(3,2) = 1.000000 0.000000
A(3,3) = 3.000000 0.000000
=====
3. cuSolverDN: 対称正定値行列のコレスキー分解を使った線形ソルバー (LAPACK Dpotrs相当)
対称正定値行列をコレスキー分解を使用して解く倍精度型ソルバー（cusolverDnDpotrs）の利用例である。コレスキー分解は、cusolverDnDpotrf ルーチンを使用し、cusolverDnDpotrs を使用して連立一次方程式を解く。なお、R.H.S（右辺） は、複数のベクトルを与えて解く。なお、密行列 A は、 column-major order で入力する必要がある。
OpenACC Fortran + cusolverDnDpotrf, cusolverDnDpotrs (Solver.f90)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved.
!
! Permission is hereby granted, free of charge, to any person obtaining a
! copy of this software and associated documentation files (the "Software"),
! to deal in the Software without restriction, including without limitation
! the rights to use, copy, modify, merge, publish, distribute, sublicense,
! and/or sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.
!
program main
use cublas_v2
use cusolverDn
use cudafor
implicit none
integer, parameter :: n=5, nrhs=3
integer, parameter :: lda=n, ldb=n
integer :: i
real(8) :: a(lda,n)
real(8) :: b(ldb, nrhs)
real(8), allocatable :: workspace_d(:)
integer :: devInfo_d
integer :: istat, Lwork
type(cusolverDnHandle) :: h
data a /3.14, 0.00, 0.00, 0.00, 0.00, &
0.17, 0.79, 0.00, 0.00, 0.00, &
-0.90, 0.83, 4.53, 0.00, 0.00, &
1.65,-0.65,-3.70, 5.32, 0.00, &
-0.72, 0.28, 1.60,-1.37, 1.98/
data b /-7.29, 9.25, 5.99,-1.94,-8.30, &
6.11, 2.90,-5.05,-3.80, 9.66, &
0.59, 8.88, 7.57, 5.57,-1.67/
print *, "Matrix A"
do i = 1, n
write(*,"(5(f8.2,2x))") a(i,:)
end do
print *, "RHS "
do i = 1, n
write(*,"(5(f8.2,2x))") b(i,:)
end do
istat = cusolverDnCreate(h)
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'handle creation failed'
!$acc data copy(a,b) create(workspace_d)
!$acc host_data use_device(a)
istat = cusolverDnDpotrf_bufferSize(h, &
CUBLAS_FILL_MODE_LOWER, n, a, n, Lwork)
!$acc end host_data
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'cusolverDnDpotrf_buffersize failed'
print *, "working space (bytes) =", Lwork
allocate(workspace_d(Lwork))
!$acc enter data copyin(workspace_d) copyin(devInfo_d)
!$acc host_data use_device(a, b, workspace_d, devInfo_d)
istat = cusolverDnDpotrf(h, CUBLAS_FILL_MODE_UPPER, &
n, a, lda, workspace_d, Lwork, devInfo_d)
istat = cusolverDnDpotrs(h, CUBLAS_FILL_MODE_UPPER, &
n, nrhs, a, lda, b, ldb, devInfo_d)
!$acc end host_data
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'cusolverDnDpotrs failed'
!$acc exit data copyout(devInfo_d)
!$acc exit data delete(workspace_d)
!$acc end data
istat = devInfo_d
if (istat /= 0) write(*,*) 'Cholesky Solver failed'
istat = cusolverDnDestroy(h)
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'handle destruction failed'
print *, "Cholesky factorization"
do i = 1, n
write(*,"(5(f8.2,2x))") a(i,:)
end do
print *, "Solution"
do i = 1, n
write(*,"(5(f8.2,2x))") b(i,:)
end do
end program main
$ pgf90 -c -O2 -Minfo -ta=tesla,cuda9.0,cc35,cc60 -Mcudalib=cusolver,cublas -Kieee solver.f90
main:
47, Loop not vectorized/parallelized: contains call
48, Loop unrolled 5 times (completely unrolled)
51, Loop not vectorized/parallelized: contains call
52, Loop unrolled 3 times (completely unrolled)
59, Generating copy(a(:,:))
Generating create(workspace_d(:))
Generating copy(b(:,:))
70, Generating enter data copyin(devinfo_d,workspace_d(:))
82, Generating exit data copyout(devinfo_d)
83, Generating exit data delete(workspace_d(:))
94, Loop not vectorized/parallelized: contains call
95, Loop unrolled 5 times (completely unrolled)
98, Loop not vectorized/parallelized: contains call
99, Loop unrolled 3 times (completely unrolled)
$ pgf90 -o a.out solver.o -Mcuda -acc -ta=tesla,cuda9.0,cc35,cc60 -Mcudalib=cusolver,cublas
$ a.out
Matrix A
3.14 0.17 -0.90 1.65 -0.72
0.00 0.79 0.83 -0.65 0.28
0.00 0.00 4.53 -3.70 1.60
0.00 0.00 0.00 5.32 -1.37
0.00 0.00 0.00 0.00 1.98
RHS
-7.29 6.11 0.59
9.25 2.90 8.88
5.99 -5.05 7.57
-1.94 -3.80 5.57
-8.30 9.66 -1.67
working space (bytes) = 2
Cholesky factorization
1.77 0.10 -0.51 0.93 -0.41
0.00 0.88 0.99 -0.84 0.36
0.00 0.00 1.81 -1.32 0.57
0.00 0.00 0.00 1.42 0.05
0.00 0.00 0.00 0.00 1.16
Solution
-6.02 3.95 -3.14
15.62 4.32 13.05
3.02 -8.25 4.91
3.25 -4.83 6.11
-8.78 9.04 -3.57
4. cuSolverDN: LU分解を使った線形システムソルバー (LAPACK Dgetrs 相当)
Dense Matrix を対象とした LU 分解を使用して解く倍精度型線形ソルバー（cusolverDnDgetrs）の利用例である。コレスキー分解は、cusolverDngetrf ルーチンを使用し、cusolverDnSgetrs を使用して連立一次方程式を解く。なお、R.H.S（右辺） は、複数のベクトルを与えて解く。なお、密行列 A は、 column-major order で入力する必要がある。
OpenACC Fortran + cusolverDnDpotrf, cusolverDnDpotrs (Solver.f90)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved.
!
! Permission is hereby granted, free of charge, to any person obtaining a
! copy of this software and associated documentation files (the "Software"),
! to deal in the Software without restriction, including without limitation
! the rights to use, copy, modify, merge, publish, distribute, sublicense,
! and/or sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.
!
program main
use cublas_v2
use cusolverDn
use cudafor
implicit none
integer, parameter :: m=5, n=3, nrhs=2
integer, parameter :: lda=m, ldb=m
integer :: i, max_mn
real(8) :: a(lda,n)
real(8) :: b(ldb, nrhs)
integer, allocatable :: devIpiv(:)
real(8), allocatable :: workspace_d(:)
integer :: devInfo_d
integer :: istat, Lwork
type(cusolverDnHandle) :: h
data a /1.00, 2.00, 3.00, 4.0, 5.0,&
1.00, 3.00, 5.00, 2.0, 4.0,&
1.00, 4.00, 2.00, 5.0, 3.0/
data b /-10.0, 12.00, 14.0, 16.0, 18.0,&
-3.0, 14.0, 12.0, 16.0,16.0/
print *, "Matrix A"
do i = 1, m
write(*,"(5(f8.2,2x))") a(i,:)
end do
print *, "RHS "
do i = 1, m
write(*,"(5(f8.2,2x))") b(i,:)
end do
max_mn = max (m, n)
allocate (devIpiv(max_mn))
istat = cusolverDnCreate(h)
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'handle creation failed'
!$acc data copy(a,b) create(workspace_d) copyout(devIpiv)
!$acc host_data use_device(a)
istat = cusolverDnDgetrf_bufferSize(h, &
m, n, a, lda, Lwork)
!$acc end host_data
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'cusolverDnDgetrf_buffersize failed'
print *, "working space (bytes) =", Lwork
allocate(workspace_d(Lwork))
!$acc enter data copyin(workspace_d) copyin(devInfo_d)
!$acc host_data use_device(a, b, workspace_d, devInfo_d, devIpiv)
istat = cusolverDnDgetrf(h, &
m, n, a, lda, workspace_d, devIpiv, devInfo_d)
istat = cusolverDnDgetrs(h, CUBLAS_OP_N, &
n, nrhs, a, lda, devIpiv, b, ldb, devInfo_d)
!$acc end host_data
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'cusolverDnDgetrs failed'
!$acc exit data copyout(devInfo_d)
!$acc exit data delete(workspace_d)
!$acc end data
istat = devInfo_d
if (istat /= 0) write(*,*) 'LU Solver failed'
istat = cusolverDnDestroy(h)
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'handle destruction failed'
print *, "LU factorization"
do i = 1, m
write(*,"(5(f8.2,2x))") a(i,:)
end do
print *, "Pivot Indcies"
print *, devIpiv(:)
print *, "Solution"
do i = 1, m
write(*,"(5(f8.2,2x))") b(i,:)
end do
end program main
$ pgf90 -c -O2 -Minfo -ta=tesla,cuda9.0,cc35,cc60 -Mcudalib=cusolver,cublas -Kieee solver.f90
main:
45, Loop not vectorized/parallelized: contains call
46, Loop unrolled 3 times (completely unrolled)
49, Loop not vectorized/parallelized: contains call
50, Loop unrolled 2 times (completely unrolled)
60, Generating copy(a(:,:))
Generating copyout(devipiv(:))
Generating create(workspace_d(:))
Generating copy(b(:,:))
71, Generating enter data copyin(devinfo_d,workspace_d(:))
83, Generating exit data copyout(devinfo_d)
84, Generating exit data delete(workspace_d(:))
95, Loop not vectorized/parallelized: contains call
96, Loop unrolled 3 times (completely unrolled)
101, Loop not vectorized/parallelized: contains call
102, Loop unrolled 2 times (completely unrolled)
pgf90 -o a.out solver.o -Mcuda -acc -ta=tesla,cuda9.0,cc35,cc60 -Mcudalib=cusolver,cublas
$ a.out
Matrix A
1.00 1.00 1.00
2.00 3.00 4.00
3.00 5.00 2.00
4.00 2.00 5.00
5.00 4.00 3.00
RHS
-10.00 -3.00
12.00 14.00
14.00 12.00
16.00 16.00
18.00 16.00
working space (bytes) = 18
LU factorization
5.00 4.00 3.00
0.60 2.60 0.20
0.40 0.54 2.69
0.80 -0.46 1.00
0.20 0.08 0.14
Pivot Indcies
5 3 3 0 0
Solution
2.00 1.20
1.14 0.74
1.14 2.34
16.00 16.00
-10.00 -3.00
5. cuSolverDN: 対称行列の固有値ソルバー (LAPACK Dsyevd相当)
対称性のある密行列に対する固有値、固有ベクトルを求めるためのソルバー（cusolverDnDsyevd）を使用する。以下の例は、D.1. Standard Symmetric Dense Eigenvalue Solvern と同じ処理を Fortran OpenACC ならびに C++ OpenACC で書いたものである。
Dense Linear Solver - cusolverDnDsyevd with OpenACC Fortran
対称（エルミート）n×n 行列 A の固有値と固有ベクトルを計算します。標準対称固有値問題は、以下の式で表される。
ここで、Λは実際の n×n 対角行列です。 V は n×n ユニタリ行列である。 Λの対角要素は、A の昇順の固有値です。cusolverDnDsyevd ソルバーは、予め作業空間（d_work）を提供する必要があるため、syevd_bufferSize 関数で、lwork サイズを求める必要がある。出力パラメータの devInfo = -i（ゼロより小さい)の場合、i 番目のパラメータは間違っていることを示す。 devInfo = i（零より大きい）場合は、内部計算時に三重対角の i 個のオフ対角要素がゼロに収束しなかったことを表す。なお、密行列 A は、 column-major order で入力する必要がある。
OpenACC + cusolverDnDsyevd (syevd.f90)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved. ! ! Permission is hereby granted, free of charge, to any person obtaining a ! copy of this software and associated documentation files (the "Software"), ! to deal in the Software without restriction, including without limitation ! the rights to use, copy, modify, merge, publish, distribute, sublicense, ! and/or sell copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER ! DEALINGS IN THE SOFTWARE. ! pgf90 -acc -O2 -Minfo -ta=tesla,cc60,cc35 -Mcuda -Mcudalib=cusolver syevd.f90 program main use cusolverDn implicit none integer,parameter :: prc = SELECTED_REAL_KIND(15,305) integer, parameter :: idbg = 1 ! debug write enable integer, parameter :: m = 3, lda =m type(cusolverDnHandle) :: cusolver_H integer :: cusolver_status integer :: status1, ierr_code integer :: i, j ! | 3.5 0.5 0 | ! A = | 0.5 3.5 0 | ! | 0 0 2 | ! ! lambda = (/ 2.0, 3.0, 4.0 /) exact eigenvalues real(prc), allocatable, dimension(:) :: A real(prc), allocatable, dimension(:) :: lambda real(prc), allocatable, dimension(:) :: V ! eigenvectors real(prc), allocatable, dimension(:) :: W ! eigenvaluse real(prc), allocatable, dimension(:) :: d_work integer :: devInfo integer :: Lwork allocate ( A(lda*m), V(lda*m), W(m) ) allocate ( lambda(m) ) Lwork = 0 ! set a data A(1:lda*m) = (/3.5, 0.5, 0, 0.5, 3.5, 0, 0, 0, 2.0/) ! step 1 create cusolver/cublas handle cusolver_status = cusolverDnCreate(cusolver_H) if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step1 cusolver_status = ", cusolver_status end if !$acc data copy(A, devInfo) copyout(W) !$acc host_data use_device(A, W) ! step 2 working space of syevd cusolver_status = cusolverDnDsyevd_bufferSize(cusolver_H, CUSOLVER_EIG_MODE_VECTOR, & CUSOLVER_EIG_MODE_VECTOR, m, A, lda, W, Lwork) !$acc end host_data if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step2 cusolver_status = ", cusolver_status end if ! step 3 d_work buffer allocated allocate(d_work(Lwork), stat=ierr_code) ! words if (idbg == 1) print *, "step5: alloc ierr_code = ", ierr_code !$acc enter data copyin(d_work) ! step 4 compute spectrum = eigen{vectors, valuse} !$acc host_data use_device(A, W, d_work, devInfo) cusolver_status = cusolverDnDsyevd(cusolver_H, CUSOLVER_EIG_MODE_VECTOR, & CUSOLVER_EIG_MODE_VECTOR, m, A, lda, W, d_work, Lwork, devInfo) !$acc end host_data !$acc exit data delete(d_work) !$acc end data print *, "Output parameter devInfo =", devInfo if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step4 cusolver_status = ", cusolver_status end if ! copy [A].output to [V] (eigenvectors) V = A ! atep 5 print out ! Eigenvalues print *, "Eigenvalue : " do i = 1, m print *, "W(",i,") = ", w(i) end do ! Eigenvectors print *, "Eigenvector : " do i = 1, m do j = 1, m print '(1x,a,i5,a,i5,1x,a4,F15.7)', "(",i,",",j,") = ", V(i+(j-1)*lda) end do end do deallocate(A, V, W, d_work) end program main
コンパイルオプションとして、-acc -Mcudalib=cusolver -Mcuda オプションを付けてコンパイルリンクする必要がある。
コンパイル＆実行
[kato@photon32 Dens_Eigen]$ pgf90 -acc -O2 -Minfo -ta=tesla,cc60,cc35 -Mcuda -Mcudalib=cusolver syevd.f90
main:
57, Generated vector simd code for the loop
Residual loop unrolled 1 times (completely unrolled)
65, Generating copy(a(:))
Generating copyout(w(:))
Generating copy(devinfo)
81, Generating enter data copyin(d_work(:))
88, Generating exit data delete(d_work(:))
98, Generated vector simd code for the loop
Residual loop unrolled 1 times (completely unrolled)
103, Loop not vectorized/parallelized: contains call
110, Loop not vectorized/parallelized: contains call
[kato@photon32 Dens_Eigen]$ a.out
step5: alloc ierr_code = 0
Output parameter devInfo = 0
Eigenvalue :
W( 1 ) = 2.000000000000000
W( 2 ) = 3.000000000000000
W( 3 ) = 4.000000000000000
Eigenvector :
( 1, 1 ) = 0.0000000
( 1, 2 ) = -0.7071068
( 1, 3 ) = 0.7071068
( 2, 1 ) = 0.0000000
( 2, 2 ) = 0.7071068
( 2, 3 ) = 0.7071068
( 3, 1 ) = 1.0000000
( 3, 2 ) = 0.0000000
( 3, 3 ) = 0.0000000
Dense Linear Solver - cusolverDnDsyevd with OpenACC C++
上記 Fortran プログラムを PGI C++ OpenACC を使用して書き直す。
OpenACC + cusolverDnDsyevd (syevd.cpp)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved. ! ! Permission is hereby granted, free of charge, to any person obtaining a ! copy of this software and associated documentation files (the "Software"), ! to deal in the Software without restriction, including without limitation ! the rights to use, copy, modify, merge, publish, distribute, sublicense, ! and/or sell copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER ! DEALINGS IN THE SOFTWARE. /* * How to compile (assume cuda is installed at /usr/local/cuda-9.0/) * pgc++ -I/usr/local/cuda-9.0/include syevd.cpp -acc -ta=tesla,cc60,cc35 -Mcudalib=cusolver -Minfo=accel -O2 * */ #include <stdio.h> #include <stdlib.h> #include <assert.h> #include <cusolverDn.h> void printMatrix(int m, int n, const double*A, int lda, const char* name) { for(int row = 0 ; row < m ; row++){ for(int col = 0 ; col < n ; col++){ double Areg = A[row + col*lda]; printf("%s(%d,%d) = %f\n", name, row+1, col+1, Areg); } } } int main(int argc, char*argv[]) { cusolverDnHandle_t cusolverH = NULL; cusolverStatus_t cusolver_status = CUSOLVER_STATUS_SUCCESS; const int m = 3; const int lda = m; /* | 3.5 0.5 0 | * A = | 0.5 3.5 0 | * | 0 0 2 | * */ // double lambda[m] = { 2.0, 3.0, 4.0}; /* Exact eigenvalues */ double *A; // ibput double *V; // eigenvectors double *W; // eigenvalues int *devInfo = NULL; double *work = NULL; int lwork = 0; A = (double *) malloc (sizeof(double) * lda * m); V = (double *) malloc (sizeof(double) * lda * m); W = (double *) malloc (sizeof(double) * m); devInfo = (int *) malloc (sizeof(int)); // step 1: create cusolver/cublas handle cusolver_status = cusolverDnCreate(&cusolverH); // step 2: set A A[0] = 3.5 ; A[1] = 0.5 ; A[2] = 0.0; A[3] = 0.5; A[4] = 3.5 ; A[5] = 0.0; A[6] = 0.0; A[7] = 0.0 ; A[8] = 2.0; printMatrix(m, m, A, lda, "A"); printf("=====\n"); // step 3: query working space of syevd cusolverEigMode_t jobz = CUSOLVER_EIG_MODE_VECTOR; // compute eigenvalues and eigenvectors. cublasFillMode_t uplo = CUBLAS_FILL_MODE_LOWER; #pragma acc data copy(A[0:lda*m], devInfo[0:1]) copyout(W[0:m]) { #pragma acc host_data use_device(A, W) { cusolver_status = cusolverDnDsyevd_bufferSize( cusolverH, jobz, uplo, m, A, lda, W, &lwork); } assert (cusolver_status == CUSOLVER_STATUS_SUCCESS); work = (double *) malloc(sizeof(double)*lwork); printf("Lwork = %d\n",lwork); #pragma acc enter data create(work[0:lwork]) // step 4: compute spectrum #pragma acc host_data use_device(A, W, work, devInfo) { cusolver_status = cusolverDnDsyevd( cusolverH, jobz, uplo, m, A, lda, W, work, lwork, devInfo); } #pragma acc wait #pragma acc exit data delete(work) } printf("DevInfo = %d\n",devInfo[0]); printf("eigenvalue = (base 1), ascending order\n"); for(int i = 0 ; i < m ; i++){ printf("W[%d] = %E\n", i+1, W[i]); } // print eigenvecvtors for(int i = 0 ; i < m*m ; i++){ V[i] = A[i]; } printf("V = (base 1)\n"); printMatrix(m, m, V, lda, "V"); printf("=====\n"); // free resources if ( A ) free(A); if ( W ) free(W); if (devInfo) free(devInfo); if (work ) free(work); if (cusolverH) cusolverDnDestroy(cusolverH); // cudaDeviceReset(); return 0; }
PGI C/C++ OpenACC コンパイラで cuSOLVER を使用する場合は、NVIDIA CUDA Toolkit にバンドルされている C/C++ 用 include ファイルを使用してコンパイルする(PGI コンパイラ内には、C/C++ 用の cuSOLVER/cuBLAS ライブラリ等の include ファイルはバンドルされていない） 。ここでの例は CUDA 9.0 が /usr/local/cuda-9.0 配下に実装されているものとして記述する。コンパイルオプションとして、 -I/usr/local/cuda-9.0/include -acc -ta=tesla -Mcudalib=cusolver オプションを付けてコンパイルリンクする必要がある。
(CUDA 8.0 の場合）
pgc++ -I/usr/local/cuda-8.0/include syevd.cpp -acc -ta=tesla,cc60,cc35 -Mcudalib=cusolver -Minfo=accel -O2
(CUDA 9.0 の場合）
pgc++ -I/usr/local/cuda-9.0/include syevd.cpp -acc -ta=tesla,cc60,cc35 -Mcudalib=cusolver -Minfo=accel -O2
コンパイル＆実行
[kato@photon32 C]$ pgc++ -I/usr/local/cuda-9.0/include syevd.cpp -acc -ta=tesla,cc60,cc35
-Mcudalib=cusolver -Minfo=accel -O2
main:
65, Generating copyout(W[:3])
Generating copy(devInfo[:1],A[:9])
87, Generating enter data create(work[:lwork])
103, Generating update self(devInfo[:1])
Generating exit data delete(work[:1])
[kato@photon32 C]$ a.out
A(1,1) = 3.500000
A(1,2) = 0.500000
A(1,3) = 0.000000
A(2,1) = 0.500000
A(2,2) = 3.500000
A(2,3) = 0.000000
A(3,1) = 0.000000
A(3,2) = 0.000000
A(3,3) = 2.000000
=====
Lwork = 83
DevInfo = 0
eigenvalue = (base 1), ascending order
W[1] = 2.000000E+00
W[2] = 3.000000E+00
W[3] = 4.000000E+00
V = (base 1)
V(1,1) = 0.000000
V(1,2) = -0.707107
V(1,3) = 0.707107
V(2,1) = 0.000000
V(2,2) = 0.707107
V(2,3) = 0.707107
V(3,1) = 1.000000
V(3,2) = 0.000000
V(3,3) = 0.000000
=====
6. cuSolverDN: QR分解による密行列線形ソルバー (LAPACK Dgeqrf 相当)
密行列線形システムを解くために QR 分解を利用する例を示す。以下の例は、C.1. QR Factorization Dense Linear Solver と同じ処理を Fortran OpenACC で書いたものである。なお、密行列 A は、 column-major order で入力する必要がある。
QR Factorization Dense Linear Solver with OpenACC Fortran
以下の線形システムを QR 分解を用いて解く。
以下のプログラムでは、3 x 3 密行列（正則行列）の例を示す。
cuSOLVEルーチン（geqrf,ormqr）、cuBLAS ルーチン（trsm）を利用する。
Step 1: A = Q*R by geqrf.
Step 2: B := Q^T*B by ormqr.
Step 3: solve R*X = B by trsm.
OpenACC + cuSOLVE-geqrf,ormqr,trsm (ormqr.f90)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved. ! ! Permission is hereby granted, free of charge, to any person obtaining a ! copy of this software and associated documentation files (the "Software"), ! to deal in the Software without restriction, including without limitation ! the rights to use, copy, modify, merge, publish, distribute, sublicense, ! and/or sell copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER ! DEALINGS IN THE SOFTWARE. ! pgf90 -acc -O2 -Minfo=accel -ta=tesla,cc60,cc35 -Mcudalib=cusolver,cublas ! -Mcuda ormqr.f90 program main use cublas_v2 use cusolverDn implicit none integer,parameter :: prc = SELECTED_REAL_KIND(15,305) integer, parameter :: idbg = 1 ! debug write enable integer, parameter :: m = 3, lda =m, ldb = m integer, parameter :: nrhs = 1 ! number of right hand side vectors real(prc) :: one = 1.0d0 type(cusolverDnHandle) :: cusolver_H type(cublasHandle) :: cublas_H integer :: cusolver_status integer :: cublas_status integer :: status1, ierr_code integer :: i, j ! | 1 2 3 | ! A = | 4 5 6 | ! | 2 1 1 | ! ! x = (1 1 1) exact solution ! b = (6 15 4) real(prc), allocatable, dimension(:) :: A real(prc), allocatable, dimension(:) :: B real(prc), allocatable, dimension(:) :: Xc ! solution matrix from GPU real(prc), allocatable, dimension(:) :: tau real(prc), allocatable, dimension(:) :: d_work integer :: devInfo integer :: Lwork allocate ( A(lda*m), B(ldb*nrhs), Xc(ldb*nrhs) ) allocate ( tau(m) ) Lwork = 0 ! set a data A(1:lda*m) = (/1.0, 4.0, 2.0, 2.0, 5.0, 1.0, 3.0, 6.0, 1.0/) B(1:ldb*nrhs) = (/6.0, 15.0, 4.0/) ! print A print *, "A ======" do i = 1, m do j = 1, m print '(1x,a,i5,a,i5,1x,a4,F15.7)', "(",i,",",j,") = ", A(i+(j-1)*lda) end do end do ! print B print *, "B ======" do i = 1, ldb do j = 1, 1 print '(1x,a,i5,a,i5,1x,a4,F15.7)', "(",i,",",j,") = ", B(i) end do end do ! step 1 create cusolver/cublas handle cusolver_status = cusolverDnCreate(cusolver_H) cublas_status = cublasCreate(cublas_H) if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step1 cusolver_status = ", cusolver_status if (cublas_status /= CUBLAS_STATUS_SUCCESS ) print *, "step1 cublas_status = ", cublas_status end if !$acc data copy(A, B, devInfo) create(tau) ! step 2 query working space of geqrf and ormqr !$acc host_data use_device(A) cusolver_status = cusolverDnDgeqrf_bufferSize(cusolver_H, & m, m, A, lda, Lwork) !$acc end host_data if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step2 cusolver_status = ", cusolver_status end if ! step 3 d_work buffer allocated allocate(d_work(Lwork), stat=ierr_code) ! words if (idbg == 1) then if (ierr_code /= 0) print *, "step3: alloc ierr_code = ", ierr_code end if print *, "Lwork = ", Lwork !$acc enter data copyin(d_work) ! step 4 compute QR factorization !$acc host_data use_device(A, tau, d_work, devInfo) cusolver_status = cusolverDnDgeqrf(cusolver_H, & m, m, A, lda, tau, d_work, Lwork, devInfo) !$acc end host_data !$acc update self(devInfo) print *, "after geqrf: devInfo =", devInfo ! step 5: compute Q^T*B !$acc host_data use_device(A, B, tau, d_work, devInfo) cusolver_status = cusolverDnDormqr(cusolver_H, CUBLAS_SIDE_LEFT, CUBLAS_OP_T, & m, nrhs, m, A, lda, tau, B, ldb, d_work, Lwork, devInfo) !$acc end host_data !$acc update self(devInfo) print *, "after ormqr: devInfo =", devInfo if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step5 cusolver_status = ", cusolver_status end if ! step 6: compute x = R \ Q^T*B !$acc host_data use_device(A, B) cublas_status = cublasDtrsm_v2(cublas_H, CUBLAS_SIDE_LEFT, CUBLAS_FILL_MODE_UPPER, CUBLAS_OP_N, & CUBLAS_DIAG_NON_UNIT, m, nrhs, one, A, lda, B, ldb) !$acc end host_data !$acc exit data delete(tau, d_work) !$acc end data if (idbg == 1) then if (cublas_status /= CUBLAS_STATUS_SUCCESS) print *, "step6 cublas_status = ", cublas_status end if ! copy [B].output to result [Xc] Xc = B ! atep 5 print out print *, "Result ===================== " do i = 1, m do j = 1, nrhs print '(1x,a,i5,a,i5,1x,a4,F15.7)', "(",i,",",j,") = ", Xc(i+(j-1)*ldb) end do end do deallocate(A, B, Xc, tau, d_work) end program main
コンパイルオプションとして、-acc -Mcudalib=cusolver,cublas -Mcuda -ta=tesla オプションを付けてコンパイルリンクする必要がある。
コンパイル＆実行
[kato@photon32 QRf]$ pgf90 -acc -O2 -Minfo=accel -ta=tesla,cc60,cc35,cuda8.0
-Mcudalib=cusolver,cublas -Mcuda ormqr.f90
main:
88, Generating copy(a(:),b(:))
Generating create(tau(:))
Generating copy(devinfo)
107, Generating enter data copyin(d_work(:))
114, Generating update self(devinfo)
123, Generating update self(devinfo)
136, Generating exit data delete(tau(:),d_work(:))
[kato@photon32 QRf]$ a.out
A ======
( 1, 1 ) = 1.0000000
( 1, 2 ) = 2.0000000
( 1, 3 ) = 3.0000000
( 2, 1 ) = 4.0000000
( 2, 2 ) = 5.0000000
( 2, 3 ) = 6.0000000
( 3, 1 ) = 2.0000000
( 3, 2 ) = 1.0000000
( 3, 3 ) = 1.0000000
B ======
( 1, 1 ) = 6.0000000
( 2, 1 ) = 15.0000000
( 3, 1 ) = 4.0000000
Lwork = 15
after geqrf: devInfo = 0
after ormqr: devInfo = 0
Result =====================
( 1, 1 ) = 1.0000000
( 2, 1 ) = 1.0000000
( 3, 1 ) = 1.0000000
7. cuSolverDN: QR分解による直交化 (LAPACK Dgeqrf,Dorgqr 相当）
QR 分解による直交化計算の例を示す。以下の例は、C.2. orthogonalization と同じ処理を Fortran OpenACC で書いたものである。なお、密行列 A は、 column-major order で入力する必要がある。
Orthgonalization by QR Factorization with OpenACC Fortran
CUDA Library を使用して QR 分解による「直交化」の計算を行う。
以下のプログラムでは、3 x 2 密行列の例を示す。
cuSOLVEルーチン（geqrf,orgqr）を利用する。
Step 1: A = Q*R by geqrf.
Step 2:form Q by orgqr.
Step 3: check if Q is unitary or not.
OpenACC + cuSOLVE-geqrf,orgqr (ortho.f90)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved. ! ! Permission is hereby granted, free of charge, to any person obtaining a ! copy of this software and associated documentation files (the "Software"), ! to deal in the Software without restriction, including without limitation ! the rights to use, copy, modify, merge, publish, distribute, sublicense, ! and/or sell copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER ! DEALINGS IN THE SOFTWARE. ! pgf90 -acc -O2 -Minfo=accel -ta=tesla,cc60,cuda8.0 -Mcudalib=cusolver,cublas ! -Mcuda ormqr.f90 -Mcuda program main use cublas_v2 use cusolverDn implicit none integer,parameter :: prc = SELECTED_REAL_KIND(15,305) integer, parameter :: idbg = 1 ! debug write enable integer, parameter :: m = 3, n = 2, lda =m real(prc) :: beta = 1.0d0, alpha = -1.0d0 type(cusolverDnHandle) :: cusolver_H type(cublasHandle) :: cublas_H integer :: cusolver_status, cusolver_status2 integer :: cublas_status integer :: status1, ierr_code integer :: i, j ! | 1 2 | ! A = | 4 5 | ! | 2 1 | ! ! x = (1 1 1) exact solution ! b = (6 15 4) real(prc), allocatable, dimension(:) :: A real(prc), allocatable, dimension(:) :: Q ! orthonormal columns real(prc), allocatable, dimension(:) :: R ! R = I - Q**T*Q real(prc) :: R_nrm2 real(prc), allocatable, dimension(:) :: tau real(prc), allocatable, dimension(:) :: d_work integer :: devInfo integer :: Lwork integer :: Lwork_geqrf integer :: Lwork_orgqr allocate ( A(lda*n), Q(lda*n), R(n*n) ) allocate ( tau(m) ) Lwork = 0 Lwork_geqrf = 0 Lwork_orgqr = 0 ! set a data A(1:lda*n) = (/1.0, 4.0, 2.0, 2.0, 5.0, 1.0/) ! print A print *, "A ======" do i = 1, m do j = 1, n print '(1x,a,i5,a,i5,1x,a4,F15.7)', "(",i,",",j,") = ", A(i+(j-1)*lda) end do end do ! step 1 create cusolver/cublas handle cusolver_status = cusolverDnCreate(cusolver_H) cublas_status = cublasCreate(cublas_H) if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step1 cusolver_status = ", cusolver_status if (cublas_status /= CUBLAS_STATUS_SUCCESS ) print *, "step1 cublas_status = ", cublas_status end if !$acc data copy(A, devInfo) create(R, tau) ! step 2 query working space of geqrf and ormqr !$acc host_data use_device(A) cusolver_status = cusolverDnDgeqrf_bufferSize(cusolver_H, & m, n, A, lda, Lwork_geqrf) cusolver_status2= cusolverDnDorgqr_bufferSize(cusolver_H, & m, n, n, A, lda, tau, Lwork_orgqr) !$acc end host_data Lwork = max(Lwork_geqrf, Lwork_orgqr) print *, "Lwork_geqrf, Lwork_orgqr, Lwork = ", Lwork_geqrf, Lwork_orgqr, Lwork if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step2 cusolver_status = ", cusolver_status if (cusolver_status2/= CUSOLVER_STATUS_SUCCESS) print *, "step2 cusolver_status2= ", cusolver_status end if ! step 3 d_work buffer allocated allocate(d_work(Lwork), stat=ierr_code) ! words if (idbg == 1) then if (ierr_code /= 0) print *, "step3: alloc ierr_code = ", ierr_code end if !$acc enter data copyin(d_work) ! step 4 compute QR factorization !$acc host_data use_device(A, tau, d_work, devInfo) cusolver_status = cusolverDnDgeqrf(cusolver_H, & m, n, A, lda, tau, d_work, Lwork, devInfo) !$acc end host_data !$acc update self(devInfo) print *, "after geqrf: devInfo =", devInfo ! step 5: compute Q !$acc host_data use_device(A, tau, d_work, devInfo) cusolver_status = cusolverDnDorgqr(cusolver_H, & m, n, n, A, lda, tau, d_work, Lwork, devInfo) !$acc end host_data !$acc update self(A, devInfo) print *, "after orgqr: devInfo =", devInfo if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step5 cusolver_status = ", cusolver_status end if ! step 6: copy A to Q ... orthonormal columns Q = A ! print Q print *, "Q ==========" do i = 1, m do j = 1, n print '(1x,a,i5,a,i5,1x,a4,F15.7)', "(",i,",",j,") = ", Q(i+(j-1)*lda) end do end do ! step 7: measure R = I - Q**T*Q R = 0.0 do i = 1, n R( i + n*(i-1) ) = 1.0 ! R(i,i) = 1 end do !$acc update device(R) ! R = -Q**T*Q + I (Q**T = CUBLAS_OP_T, Q = CUBLAS_OP_T) !$acc host_data use_device(A, R) cublas_status = cublasDgemm_v2(cublas_H, CUBLAS_OP_T, CUBLAS_OP_N, & n, n, m, alpha, A, lda, A, lda, beta, R, n) !$acc end host_data !$acc update self(R) if (idbg == 1) then if (cublas_status /= CUBLAS_STATUS_SUCCESS) print *, "step7 cublas_status = ", cublas_status end if R_nrm2 = 0.0d0 !$acc host_data use_device(R) cublas_status = cublasDnrm2_v2(cublas_H, n*n, R, 1, R_nrm2) !$acc end host_data !$acc exit data delete(tau, d_work) !$acc end data ! atep 5 print |I - Q**T*Q| print *, "Result ==== |I - Q**T*Q| === " print '(1x, a, 2x, D16.6)', '|I - Q**T*Q| =', R_nrm2 deallocate(A, R, tau, d_work) end program main
コンパイルオプションとして、-acc -Mcudalib=cusolver,cublas -Mcuda -ta=tesla オプションを付けてコンパイルリンクする必要がある。
コンパイル＆実行
[kato@photon32 Ortho]$ pgf90 -acc -Minfo=accel -ta=tesla,cc60,cc35,cuda8.0
-Mcudalib=cusolver,cublas ortho.f90 -Mcuda
main:
85, Generating copy(a(:),devinfo)
Generating create(r(:),tau(:))
109, Generating enter data copyin(d_work(:))
116, Generating update self(devinfo)
125, Generating update self(a(:),devinfo)
148, Generating update device(r(:))
155, Generating update self(r(:))
166, Generating exit data delete(tau(:),d_work(:))
[kato@photon32 Ortho]$ a.out
A ======
( 1, 1 ) = 1.0000000
( 1, 2 ) = 2.0000000
( 2, 1 ) = 4.0000000
( 2, 2 ) = 5.0000000
( 3, 1 ) = 2.0000000
( 3, 2 ) = 1.0000000
Lwork_geqrf, Lwork_orgqr, Lwork = 14 2 14
after geqrf: devInfo = 0
after orgqr: devInfo = 0
Q ==========
( 1, 1 ) = -0.2182179
( 1, 2 ) = 0.5345225
( 2, 1 ) = -0.8728716
( 2, 2 ) = 0.2672612
( 3, 1 ) = -0.4364358
( 3, 2 ) = -0.8017837
Result ==== |I - Q**T*Q| ===
|I - Q**T*Q| = 0.450975D-15
8. cuSolverDN: 特異ベクトル計算を含めた特異値分解(by QR) (LAPACK Dgesvd 相当）
複素数あるいは実数を成分とする行列に対する特異値分解を行う例である。ここでは、実数値を扱う例を示す。LAPACK の dgesvd ルーチンに相当する。特異値分解は、信号処理や統計学の分野で用いられる。以下の例は、F.1. SVD with singular vectors と同じ処理を Fortran OpenACC で書いたものである。なお、密行列 A は、 column-major order で入力する必要がある。
以下のプログラムでは、3 x 2 密行列の例を示す。SVD は以下の式で表される。ここで、Σ はその min（m、n）対角要素を除いてゼロである m×n 行列であり、U は m×m ユニタリ行列であり、V は n×n ユニタリ行列である。 Σ の対角要素は A の特異値となる。 それらは実数で、負数ではなく、降順で返される。 U と V の最初の min（m、n）列は、A の左右の特異ベクトルとなる。cusolverDnDgesvd ルーチは、QR アルゴリズムを使用する。もう一つ cusolverDngesvdj も用意されているが、これも同じ機能を有している。ただし、後者はヤコビ法を使用して解く。中小の行列を対象とする場合は、ヤコビ法の並列性により高い GPU 性能を享受できるかもしれない。また、ヤコビ法を使用するとある精度までの近似解を得る場合に便利である。ここでは、QR アルゴリズムを用いて解く、特異値分解(SVD)の例を示すこととする。
次の三つの計算ステップを実装する。
Step 1: A = U*S*VT 特異値分解を行う。
Step 2: 特異値の値を検証する。
Step 3: 残留誤差 A-U*S*VT を cuBLAS を使って計算する。
Step 3の残留誤差を求める際に、BLAS の NVIDIA extention 関数を使用する。これは、PGI Fortran の Interface Modul には含まれていないため、module cuBLASextというモジュールを作成し、Fortran と cuBLAS 関数とのインタフェースを行なった。
OpenACC + cuSOLVER-cusolverDnDgesvd (svd_simgular.f90)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved.
! Copyright (c) 2017, NVIDIA CORPORATION. All rights reserved.
!
!
! Permission is hereby granted, free of charge, to any person obtaining a
! copy of this software and associated documentation files (the "Software"),
! to deal in the Software without restriction, including without limitation
! the rights to use, copy, modify, merge, publish, distribute, sublicense,
! and/or sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.
!
module cuBLASext
! NVIDIA extention : the BLAS-extension functions
! because PGI Fortran don't include NVIDIA extention BLAS module.
use cublas_v2
interface cublasDdgmm
integer(c_int) function cublasDdgmm( handle, side_mode, m, n, &
A, lda, x , incx , C, ldc) &
bind(C,name="cublasDdgmm")
use iso_c_binding
import cublasHandle
type(cublasHandle), value :: handle
integer(c_int), value :: side_mode, m, n, lda, ldc, incx
real(c_double), device :: A(*), C(*)
real(c_double), device :: x(*)
end function cublasDdgmm
end interface cublasDdgmm
end module cuBLASext
module aux_routine
contains
subroutine print_D_matrix(m, n, a, lda)
integer :: m, n, i, j
real(8) :: a(*)
do i=1, m
do j=1, n
print "(5(,i5,1x,i5,2x,f10.6,2x))", i, j, a(i+(j-1)*lda)
end do
end do
end subroutine print_D_matrix
end module aux_routine
program main
use cublas_v2
use cuBLASext
use cusolverDn
use aux_routine
implicit none
integer, parameter :: m=3, n=2
integer, parameter :: lda = m
integer, parameter :: idbg = 1
real(8), parameter :: one = 1.d0, minus_one = -1.d0
real(8) :: a (lda*n) ! column-major format with dimensions m x n
real(8) :: u (lda*m) ! m by m unitary matrix
real(8) :: vt(lda*n) ! n by n unitary matrix
real(8) :: w (lda*n) ! w = s*vt
real(8), allocatable, dimension(:) :: s ! singular values
real(8), allocatable, dimension(:) :: rwork ! min(m,n)-1
real(8), allocatable :: workspace_d(:)
real(8) :: S_exact(n)
real(8) :: err, diff_s, diff_fro
integer :: devInfo_d
integer :: istat, Lwork
integer :: i
character*1 :: jobu, jobvt
type(cusolverDnHandle) :: cusolver_H
type(cublasHandle) :: cublas_H
integer :: cusolver_status
integer :: cublas_status
! | 1 2 |
! a = | 4 5 |
! | 2 1 |
!
allocate(s(1:min(m,n)))
allocate(rwork(1:min(m,n)-1))
Lwork = 0
! set a data
a(1:lda*n) = (/1.0, 4.0, 2.0, 2.0, 5.0, 1.0/)
! Exact SV numbers
S_exact(1:n) =(/7.065283497082729d0, 1.040081297712078d0/)
cusolver_status = cusolverDnCreate(cusolver_H)
if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) &
print *, 'handle creation failed'
cublas_status = cublasCreate(cublas_H)
if (cublas_status /= CUBLAS_STATUS_SUCCESS ) &
print *, "cublas_status = ", cublas_status
!$acc data copy(a) create(s, u, vt, devInfo_d, workspace_d, rwork, w)
! Step 1: compute A = U*S*VT
!$acc host_data use_device(a)
istat = cusolverDnDgesvd_bufferSize(cusolver_H, m, n, Lwork)
!$acc end host_data
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'cusolverDnDgesvd_bufferSize failed'
print *, "working space (words) =", Lwork
allocate(workspace_d(Lwork))
!$acc enter data copyin(workspace_d) copyin(devInfo_d)
jobu = 'A' ! all m columns of U
jobvt= 'A' ! all m columns of VT
!$acc host_data use_device(a, s, u, vt, devInfo_d, workspace_d, rwork)
istat = cusolverDnDgesvd (cusolver_H, jobu, jobvt, &
m, n, a, lda, s, u, lda, vt, lda, workspace_d, Lwork, rwork, devInfo_d)
!$acc end host_data
!$acc wait
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'cusolverDnDgesvd failed', istat
! copyout u, vt, s, devInfo_d to host
!$acc update self(u, vt, s, devInfo_d)
! print
print *, "After gesvd: devInfo_d =", devInfo_d
if (devInfo_d /= 0) write(*,*) 'SVD operation failed'
print *, "S ="
call print_D_matrix(n, 1, s, 0)
print *, "U ="
call print_D_matrix(m, m, u, lda)
print *, "VT ="
call print_D_matrix(n, n, vt, lda)
istat = cusolverDnDestroy(cusolver_H)
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'handle destruction failed'
! Step 2: check accuracy of singular value
diff_s = 0.d0
diff_fro = 0.d0
do i = 1, n
err = abs( s(i) - s_exact(i) )
diff_s = max( err, diff_s )
end do
print '(1x, a, d17.5)',"|S - S_exact| = ", diff_s
! Step 3: measure residual A-U*S*VT
! w = s*vT
!$acc host_data use_device(s, vt, w)
cublas_status = cublasDdgmm(cublas_H, CUBLAS_SIDE_LEFT, n, n, &
vt, lda, s, 1, w, lda)
!$acc end host_data
if (idbg == 1) then
!$acc update self(w)
print *, "W ="
call print_D_matrix(n, n, w, lda)
end if
! Original A matrix is copied to device
!$acc update device(a)
! a := -u*w + a
!$acc host_data use_device(a, u, w)
cublas_status = cublasDgemm_v2(cublas_H, CUBLAS_OP_N, CUBLAS_OP_N, &
m, n, n, minus_one, u, lda, w, lda, one, a, lda)
cublas_status = cublasDnrm2_v2(cublas_H, lda*n, a, 1, diff_fro)
!$acc end host_data
if (idbg == 1) then
!$acc update self(a)
print *, "A = "
call print_D_matrix(m, n, a, lda)
end if
print '(1x, a, d17.5)',"|A - U*S*VTt| = ", diff_fro
!$acc exit data delete(workspace_d)
!$acc end data
end program main
コンパイルオプションとして、-acc -Mcudalib=cusolver,cublas -Mcuda -ta=tesla オプションを付けてコンパイルリンクする必要がある。
コンパイル＆実行
$ make run
pgf90 -c -O2 -Minfo -ta=tesla,cuda9.0,cc35,cc60 -Mcudalib=cusolver,cublas -Kieee -acc -Mcuda svd_simgular.f90
print_d_matrix:
47, Loop not vectorized/parallelized: contains call
main:
92, Loop unrolled 6 times (completely unrolled)
Loop not vectorized: loop count too small
95, Loop unrolled 2 times (completely unrolled)
104, Generating create(rwork(:),s(:),u(:))
Generating copy(a(:))
Generating create(w(:),workspace_d(:),devinfo_d,vt(:))
117, Generating enter data copyin(workspace_d(:),devinfo_d)
132, Generating update self(devinfo_d,vt(:),u(:),s(:))
152, Loop unrolled 2 times (completely unrolled)
167, Generating update self(w(:))
173, Generating update device(a(:))
183, Generating update self(a(:))
190, Generating exit data delete(workspace_d(:))
===== Program run =====
a.out
working space (words) = 131
After gesvd: devInfo_d = 0
S =
1 1 7.065283
2 1 1.040081
U =
1 1 -0.308219
1 2 0.488195
1 3 0.816497
2 1 -0.906133
2 2 0.110706
2 3 -0.408248
3 1 -0.289695
3 2 -0.865685
3 3 0.408248
VT =
1 1 -0.638636
1 2 -0.769509
2 1 -0.769509
2 2 0.638636
|S - S_exact| = 0.22204D-15
W =
1 1 -4.512143
1 2 -5.436800
2 1 -0.800352
2 2 0.664233
A =
1 1 0.000000
1 2 0.000000
2 1 0.000000
2 2 0.000000
3 1 0.000000
3 2 0.000000
|A - U*S*VTt| = 0.23604D-14
9. cuSolverDN: 特異ベクトル計算を含めた特異値分解(by Jacobi) (LAPACK Dgesvd 相当）
前段では、QR アルゴリズムを使用して解く特異値分解を行う例を示した。ここでは、ヤコビ法を使用した特異値分解を行う例を示す。LAPACK の dgesvd ルーチンに相当する中小の行列を対象とする場合は、ヤコビ法の並列性により高い GPU 性能を享受できるかもしれない。また、ヤコビ法を使用するとある精度までの近似解を得る場合に便利である。以下の例は、F.2. SVD with singular vectors (via Jacobi method) と同じ処理を Fortran OpenACC で書いたものである。なお、密行列 A は、 column-major order で入力する必要がある。
以下のプログラムでは、3 x 2 密行列の例を示す。SVD は以下の式で表される。ここで、Σ はその min（m、n）対角要素を除いてゼロである m×n 行列であり、U は m×m ユニタリ行列であり、V は n×n ユニタリ行列である。 Σ の対角要素は A の特異値となる。 それらは実数で、負数ではなく、降順で返される。 U と V の最初の min（m、n）列は、A の左右の特異ベクトルとなる。
PGI 17.7 から提供されたcuSOLVER 用の Fortran interface module(cusolverDn.mod) の中に、ヤコビ法を使用した特異値分解ルーチン cusolverDnDgesvdj を使用するための CUDA C/C++ cuSOLVER ライブラリへの Interface が含まれていないことがわかった。そのため自作した。以下に示す module cusolverDn_gesvdj モジュールの中で C ライブラリを使用するための Interface を定義した。これによって、cuSOLVER の cusolverDnDgesvdj ルーチンを Fortran から使用可能となる。
以下のプログラム例では、cusolverDnSetStream() で stream tag を作成し、CUDAストリームで cusolverDN ライブラリタスクを重複して動作させることができるようにしている。このプログラムは単一ストリームであるため、この行為自体は意味がないが、CUDAストリームを使用することで、多重のライブラリタスクを動作させるようにプログラムを組みことができる。ここでは、CUDAストリームの使い方を示したに過ぎない。CUDA ストリームを使用した際の注意点は、必ず、後続処理で同期を取らなければいけないポイントにおいて、cudaDeviceSynchronize() 関数を使い同期させることが必要である。この関数は、CUDA (Fortran) API を直接使うことが必要であり、OpenACC レベルでの !$acc wait では同期は取れないことに注意して欲しい。
OpenACC + cuSOLVER-cusolverDnDgesvdj (jacobi_svd.f90)
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved.
! Copyright (c) 2017, NVIDIA CORPORATION. All rights reserved.
!
!
! Permission is hereby granted, free of charge, to any person obtaining a
! copy of this software and associated documentation files (the "Software"),
! to deal in the Software without restriction, including without limitation
! the rights to use, copy, modify, merge, publish, distribute, sublicense,
! and/or sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.
!
module cusolverDn_gesvdj
use cudafor
use cusolverDN
type gesvdjInfo
type(c_ptr) :: svInfo
end type gesvdjInfo
! cusolverDnCreateGesvdjInfo(info)
interface
integer(c_int) function cusolverDnCreateGesvdjInfo(info) bind(C,name='cusolverDnCreateGesvdjInfo')
import gesvdjInfo
type(gesvdjInfo) :: info
end function cusolverDnCreateGesvdjInfo
end interface
! cusolverDnXgesvdjSetTolerance(info, tolerance)
interface
integer(c_int) function cusolverDnXgesvdjSetTolerance(info, tolerance) &
bind(C,name='cusolverDnXgesvdjSetTolerance')
use iso_c_binding
import gesvdjInfo
type(gesvdjInfo) :: info
real(c_double), value :: tolerance
end function cusolverDnXgesvdjSetTolerance
end interface
! cusolverDnXgesvdjSetMaxSweeps(info, max_sweeps)
interface
integer(c_int) function cusolverDnXgesvdjSetMaxSweeps(info, max_sweeps) &
bind(C,name='cusolverDnXgesvdjSetMaxSweeps')
use iso_c_binding
import gesvdjInfo
type(gesvdjInfo) :: info
integer(c_int),value :: max_sweeps
end function cusolverDnXgesvdjSetMaxSweeps
end interface
! cusolverDnDgesvdj_bufferSize
interface
integer(c_int) function cusolverDnDgesvdj_bufferSize(cusolver_Hndl, jobz, econ, &
m, n, a, lda, s, u, ldu, v, ldv, lwork, info) bind(C,name='cusolverDnDgesvdj_bufferSize')
use iso_c_binding
import cusolverDnHandle
import gesvdjInfo
type(cusolverDnHandle), value :: cusolver_Hndl
integer(c_int), value :: jobz
integer(c_int),value :: econ, m, n, lda, ldu, ldv
real(c_double), device :: a(*), s(*), u(*), v(*)
type(gesvdjInfo) :: info
end function cusolverDnDgesvdj_bufferSize
end interface
! cusolverDnDgesvdj
interface
integer(c_int) function cusolverDnDgesvdj(cusolver_Hndl, jobz, econ, &
m, n, a, lda, s, u, ldu, v, ldv, work, lwork, devinfo, info ) bind(C,name='cusolverDnDgesvdj')
use iso_c_binding
import cusolverDnHandle
import gesvdjInfo
type(cusolverDnHandle), value :: cusolver_Hndl
integer(c_int), value :: jobz
integer(c_int),value :: econ, m, n, lda, ldu, ldv
real(c_double), device :: a(*), s(*), u(*), v(*), work(*)
integer(c_int) :: devinfo
type(gesvdjInfo) :: info
end function cusolverDnDgesvdj
end interface
! cusolverDnXgesvdjGetSweeps
interface
integer(c_int) function cusolverDnXgesvdjGetSweeps(cusolver_Hndl, info, executed_sweeps) &
bind(C,name='cusolverDnXgesvdjGetSweeps')
use iso_c_binding
import cusolverDnHandle
import gesvdjInfo
type(cusolverDnHandle), value :: cusolver_Hndl
type(gesvdjInfo) :: info
integer(c_int) :: executed_sweeps
end function cusolverDnXgesvdjGetSweeps
end interface
! cusolverDnXgesvdjGetResidual
interface
integer(c_int) function cusolverDnXgesvdjGetResidual(cusolver_Hndl, info, residual) &
bind(C,name='cusolverDnXgesvdjGetResidual')
use iso_c_binding
import cusolverDnHandle
import gesvdjInfo
type(cusolverDnHandle), value :: cusolver_Hndl
type(gesvdjInfo) :: info
real(c_double) :: residual
end function cusolverDnXgesvdjGetResidual
end interface
end module cusolverDn_gesvdj
module aux_routine
contains
subroutine print_D_matrix(m, n, a, lda)
integer :: m, n, i, j
real(8) :: a(*)
do i=1, m
do j=1, n
print "(5(,i5,1x,i5,2x,D20.14,2x))", i, j, a(i+(j-1)*lda)
end do
end do
end subroutine print_D_matrix
end module aux_routine
program main
use cudafor ! cudaStreamCreateWithFlags in the cudafor.mod
use cusolverDn
use cusolverDn_gesvdj
use aux_routine
implicit none
integer, parameter :: m=3, n=2
integer, parameter :: lda = m
integer, parameter :: idbg = 1
real(8), parameter :: one = 1.d0, minus_one = -1.d0
real(8) :: a (lda*n) ! column-major format with dimensions m x n
real(8) :: u (lda*m) ! m by m unitary matrix
real(8) :: vt(lda*n) ! n by n unitary matrix
real(8), allocatable, dimension(:) :: s ! singular value
real(8), allocatable :: workspace_d(:)
real(8) :: S_exact(n)
real(8) :: err, diff_s
integer :: devInfo_d
integer :: istat, Lwork
integer :: i
! configuration of GESVDJ
real(8) :: tol
integer :: max_sweeps
integer, parameter :: econ=0 ! econ = 1 for economy size
! numerical results of GESVDJ
real(8) :: residual
integer :: executed_sweeps
! cuda library variables
type(cusolverDnHandle) :: cusolver_H
type(gesvdjInfo) :: gesvdj_params
integer :: jobz
integer(kind=cuda_stream_kind) :: stream
integer :: cusolver_status
integer :: cublas_status
! | 1 2 |
! a = | 4 5 |
! | 2 1 |
!
! cusolverDnDgesvdj() argument list
! jobz CUSOLVER_EIG_MODE_NOVECTOR: compute singular values only
! CUSOLVER_EIG_MODE_VECTOR: compute singular value and singular vectors
! econ econ = 1 for economy size
! m nubmer of rows of A, 0 <= m
! n number of columns of A, 0 <= n
! a m-by-n
! lda leading dimension of A
! s min(m,n)
! the singular values in descending order
! u m-by-m if econ = 0
! m-by-min(m,n) if econ = 1
! ldu leading dimension of U, ldu >= max(1,m)
! vt n-by-n if econ = 0
! n-by-min(m,n) if econ = 1
! ldv leading dimension of V, ldv >= max(1,n)
! Initialize
allocate(s(1:min(m,n)))
Lwork = 0
tol = 1.d-7
max_sweeps = 15
residual = 0.d0
executed_sweeps = 0
jobz = CUSOLVER_EIG_MODE_VECTOR
print *,"Eigmode=", jobz
print *, "tol =", tol
print *, "max sweeps =", max_sweeps
print *, "econ = ", econ
! set a data
a(1:lda*n) = (/1.0, 4.0, 2.0, 2.0, 5.0, 1.0/)
print *, "A ="
call print_D_matrix(m, n, a, lda)
! Exact SV numbers
S_exact(:) =(/7.065283497082729d0, 1.040081297712078d0/)
!$acc data copy(a) create(s, u, vt, devInfo_d)
! Step 1: create handle and stream
cusolver_status = cusolverDnCreate(cusolver_H)
if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) &
print *, 'handle creation failed'
! set a CUDA Lib stream -- when performing several small independent computations(非同期モードで実行する）
! you explicitly need to set cudaDeviceSynchronize() at an appropriate place by CUDA-API
istat = cudaStreamCreateWithFlags(stream, cudaStreamNonBlocking)
istat = cusolverDnSetStream(cusolver_H, stream)
! step 2 : configure gesvdj, Tolerance, Max sweeps
istat = cusolverDnCreateGesvdjInfo(gesvdj_params)
istat = cusolverDnXgesvdjSetTolerance(gesvdj_params, tol)
istat = cusolverDnXgesvdjSetMaxSweeps(gesvdj_params, max_sweeps)
! Step 3: query workspace of SVD, allocate workspace_d(Lwork)
!$acc host_data use_device(a, s, u, vt)
istat = cusolverDnDgesvdj_bufferSize(cusolver_H, jobz, econ, &
m, n, a, lda, s, u, lda, vt, lda, Lwork, gesvdj_params)
!$acc end host_data
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'cusolverDnDgesvdj_bufferSize failed'
print *, "working space (words) =", Lwork
allocate(workspace_d(Lwork))
!$acc enter data copyin(workspace_d) copyin(devInfo_d)
! step 4 : compute SVD
!$acc host_data use_device(a, s, u, vt, devInfo_d, workspace_d)
istat = cusolverDnDgesvdj (cusolver_H, jobz, econ, &
m, n, a, lda, s, u, lda, vt, lda, workspace_d, Lwork, devInfo_d, gesvdj_params)
!$acc end host_data
! Synchronize for cusolverDnDgesvdj completion (必ず必要)
istat =cudaDeviceSynchronize() ! CUDA-API, it can't use a OpenACC wait directive
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'cusolverDnDgesvd failed', istat
! copyout u, vt, s, devInfo_d to host
!$acc update self(s, u, vt, devInfo_d)
! Step 5: measure residual |A - U*S*V**H|_F
istat = cusolverDnXgesvdjGetSweeps(cusolver_H, gesvdj_params, executed_sweeps)
istat = cusolverDnXgesvdjGetResidual(cusolver_H, gesvdj_params, residual)
!$acc exit data delete(workspace_d)
!$acc end data
! print out
print *, "After gesvd: devInfo_d =", devInfo_d
if (devInfo_d == 0) then
write(*,*) 'SVD operation was converged'
else if (devInfo_d < 0) then
write(*,*) -devInfo_d, "d-th parameter is wrong"
stop
else
write(*,*) "WARNING: gesvdj does not converge ", devInfo_d
end if
print *, "Singular values ="
call print_D_matrix(n, 1, s, 0)
print *, "U = left singular vectors"
call print_D_matrix(m, m, u, lda)
print *, "VT = right singular vectors"
call print_D_matrix(n, n, vt, lda)
! check accuracy of singular value
diff_s = 0.d0
do i = 1, n
err = abs( s(i) - s_exact(i) )
diff_s = max( err, diff_s )
end do
print '(1x, a, d17.5)',"|S - S_exact| = ", diff_s
! measure residual A-U*S*VT**H
print '(1x, a, d17.5)',"|A - U*S*V**H| = ", residual
print '(1x, a, I5 )',"number of executed sweeps = ", executed_sweeps
istat = cusolverDnDestroy(cusolver_H)
if (istat /= CUSOLVER_STATUS_SUCCESS) &
write(*,*) 'handle destruction failed'
end program main
コンパイルオプションとして、-acc -Mcudalib=cusolver -Mcuda -ta=tesla オプションを付けてコンパイルリンクする必要がある。
コンパイル＆実行
$ make run
pgf90 -c -O2 -Minfo -ta=tesla,cuda9.0,cc35,cc60 -Mcudalib=cusolver,cublas -Kieee -acc -Mcuda jacobi_svd.f90
print_d_matrix:
126, Loop not vectorized/parallelized: contains call
main:
210, Loop unrolled 6 times (completely unrolled)
Loop not vectorized: loop count too small
215, Loop unrolled 2 times (completely unrolled)
217, Generating create(devinfo_d,s(:),vt(:),u(:))
Generating copy(a(:))
249, Generating enter data copyin(workspace_d(:),devinfo_d)
264, Generating update self(vt(:),u(:),s(:),devinfo_d)
271, Generating exit data delete(workspace_d(:))
296, Loop unrolled 2 times (completely unrolled)
pgf90 -o a.out jacobi_svd.o -ta=tesla,cuda9.0,cc35,cc60 -Mcudalib=cusolver -acc -Mcuda
===== Program run =====
a.out
Eigmode= 1
tol = 9.9999999999999995E-008
max sweeps = 15
econ = 0
A =
1 1 0.10000000000000D+01
1 2 0.20000000000000D+01
2 1 0.40000000000000D+01
2 2 0.50000000000000D+01
3 1 0.20000000000000D+01
3 2 0.10000000000000D+01
working space (words) = 3168
After gesvd: devInfo_d = 0
SVD operation was converged
Singular values =
1 1 0.70652834970827D+01
2 1 0.10400812977121D+01
U = left singular vectors
1 1 0.30821892063279D+00
1 2 -.48819507401990D+00
1 3 0.81649658092773D+00
2 1 0.90613333377729D+00
2 2 -.11070553170904D+00
2 3 -.40824829046386D+00
3 1 0.28969549251172D+00
3 2 0.86568461633075D+00
3 3 0.40824829046386D+00
VT = right singular vectors
1 1 0.63863583713640D+00
1 2 0.76950910814953D+00
2 1 0.76950910814953D+00
2 2 -.63863583713640D+00
|S - S_exact| = 0.88818D-15
|A - U*S*V**H| = 0.36550D-14
number of executed sweeps = 1
10. cuSolverSP: 並列ダイレクト・スパース・ソルバー (CSR QR分解)
PGI 17.7 においては、 cuSolverDN type ルーチン群に対する Fortran MODULE インタフェースが提供されたが、cuSolverSP ならびに cuSolverRF 型のライブラリに関してはまだ提供されていない。今後、順次提供されるものと思うが、ここでは、cuSolverSP 用の Fortran MODULE インタフェースを作成して、スパース QR 分解を行うプログラム例を提示する。ここで扱う例は、NVIDIA cuSOLVER の reference manual の中に例示している Batched Sparse QR example 1 を筆者が Fortran 言語で書き換えたものを示す。なお、ここでのライブラリは、複数の線形システムを解く「バッチ処理用ルーチン」を使用している。ここで使用されているバッチ・ソルバーに関する技術的な説明は、NVIDIA Parallel FORALL：Parallel Direct Solvers with cuSOLVER: Batched QR のブログでも紹介されている。
| 1 | Marices A (small perturbations 摂動) を複数用意する。
A = | 2 | バッチ処理を行うため
| 3 |
| 0.1 0.1 0.1 4 |
このプログラムを、CUDA Fortran ならびに OpenACC Fortran から利用する場合を掲示する。
配列、変数のデバイスメモリのマネージメントをどちら側のプログラミングモデルで行うかによって、プログラミングの対応は異なるが、ここでは、 CUDA Fortran 側でデバイス配列の宣言を行った場合の方法、OpenACC 側で配列の宣言を行った場合を以下に例示する。
Batched Sparse QR example with CUDA Fortran
PGI 17.7 時点では CUDA C 言語で開発されている cuSOLVER SP(スパース）系のライブラリへの Fortran Module インタフェースが提供されていないため、以下のように自作(cusolver_mod.cuf)した。C と Fortran のバインディングを取り持つためのインタフェースである。将来的に、、PGI コンパイルシステムの中に、cuSOLVER SP の Fortran MODULE (多分 use cusplverSP で使用できるはず）が提供された場合は、以下の cusolver_mod.cuf ルーチンは必要なくなる。まず最初に、CUDA Fortran で coding した例を提示する。この後で、このプログラム変更し OpenACC を使って coding した例を示すことにする。なお、使用したバッチ用ライブラリルーチン(cusolverSpXcsrqrBatched())の仕様は、こちらを参考にして欲しい。この関数の引数仕様の説明の中で、host or device の区別があるが、これは、どちら側に置かれた変数かを示すものである。これを見ることによって、デバイス変数として対処する必要があるかどうかが分かる。
cuSOLVER SP QR 用の Fortran Module Interface（cusolver_mod.cuf）
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved.
!
! Permission is hereby granted, free of charge, to any person obtaining a
! copy of this software and associated documentation files (the "Software"),
! to deal in the Software without restriction, including without limitation
! the rights to use, copy, modify, merge, publish, distribute, sublicense,
! and/or sell copies of the Software, and to permit persons to whom the
! Software is furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
! DEALINGS IN THE SOFTWARE.
!
module cuSOLVER_SP_QR
use cudafor
use cusparse
enum, bind(C) ! cusolverStatus_t
enumerator :: CUSOLVER_STATUS_SUCCESS=0
enumerator :: CUSOLVER_STATUS_NOT_INITIALIZED=1
enumerator :: CUSOLVER_STATUS_ALLOC_FAILED=2
enumerator :: CUSOLVER_STATUS_INVALID_VALUE=3
enumerator :: CUSOLVER_STATUS_ARCH_MISMATCH=4
enumerator :: CUSOLVER_STATUS_MAPPING_ERROR=5
enumerator :: CUSOLVER_STATUS_EXECUTION_FAILED=6
enumerator :: CUSOLVER_STATUS_INTERNAL_ERROR=7
enumerator :: CUSOLVER_STATUS_MATRIX_TYPE_NOT_SUPPORTED=8
enumerator :: CUSOLVER_STATUS_NOT_SUPPORTED = 9
enumerator :: CUSOLVER_STATUS_ZERO_PIVOT=10
enumerator :: CUSOLVER_STATUS_INVALID_LICENSE=11
end enum
enum, bind(C) ! cusolverEigType_t
enumerator :: CUSOLVER_EIG_TYPE_1=1
enumerator :: CUSOLVER_EIG_TYPE_2=2
enumerator :: CUSOLVER_EIG_TYPE_3=3
end enum
type cusolverSpContext
type(c_ptr) :: cusolverSpHandle
end type cusolverSpContext
type csrqrInfo
type(c_ptr) :: Info
end type csrqrInfo
type cusolverSpHandle
type(c_ptr) :: SpHandle
end type cusolverSpHandle
! cusolverSpCreate
interface
integer(c_int) function cusolverSpCreate(handle) bind(C,name='cusolverSpCreate')
import cusolverSpHandle
type(cusolverSpHandle) :: handle
end function cusolverSpCreate
end interface
! cusolverSpDestroy
interface
integer(c_int) function cusolverSpDestroy(handle) bind(C,name='cusolverSpDestroy')
import cusolverSpHandle
type(cusolverSpHandle), value :: handle
end function cusolverSpDestroy
end interface
! cusolverSpCreateCsrqrInfo
interface
integer(c_int) function cusolverSpCreateCsrqrInfo(info) bind(C,name='cusolverSpCreateCsrqrInfo')
import csrqrInfo
type(csrqrInfo) :: info
end function cusolverSpCreateCsrqrInfo
end interface
! cusolverSpDestroyCsrqrInfo
interface
integer(c_int) function cusolverSpDestroyCsrqrInfo(info) bind(C,name='cusolverSpDestroyCsrqrInfo')
import csrqrInfo
type(csrqrInfo) :: info
end function cusolverSpDestroyCsrqrInfo
end interface
! cusolverSpDcsrqrBufferInfoBatched
interface cusolverSpDcsrqrBufferInfoBatched
integer(c_int) function cusolverSpDcsrqrBufferInfoBatched( cusolver_Hndl, m, n, nnzA, &
descrA, csrValA, csrRowPtrA, csrColIndA, BatchSize, info, internalDataInBytes, workspaceInBytes ) &
bind(C,name="cusolverSpDcsrqrBufferInfoBatched")
use iso_c_binding
import cusolverSpHandle, cusparseMatDescr, csrqrInfo
type(cusolverSpHandle), value :: cusolver_Hndl
type(cusparseMatDescr), value :: descrA
type(csrqrInfo), value :: info
integer(c_int), value :: m, n, nnzA, batchSize
real(c_double), device :: csrValA(*)
!! !pgi$ ignore_tkr (d) csrRowPtrA, (d) csrColIndA
integer(c_int),device :: csrRowPtrA(*), csrColIndA(*)
!! !pgi$ ignore_tkr (k) internalDataInBytes, (k) workspaceInBytes
integer(8) :: internalDataInBytes, workspaceInBytes
end function cusolverSpDcsrqrBufferInfoBatched
end interface cusolverSpDcsrqrBufferInfoBatched
! cusolverSpXcsrqrAnalysis
!interface cusolverSpXcsrqrAnalysis
! integer(c_int) function cusolverSpXcsrqrAnalysis( cusolver_Hndl, m, n, nnzA, &
! descrA, csrRowPtrA, csrColIndA, info ) &
! bind(C,name="cusolverSpXcsrqrAnalysis")
! use iso_c_binding
! import cusolverSpHandle, cusparseMatDescr, csrqrInfo
! type(cusolverSpHandle), value :: cusolver_Hndl
! type(cusparseMatDescr), value :: descrA
! type(csrqrInfo), value :: info
! integer(c_int), value :: m, n, nnzA
! !pgi$ ignore_tkr (d) csrRowPtrA, (d) csrColIndA
! integer(c_int),device :: csrRowPtrA(*), csrColIndA(*)
! end function cusolverSpXcsrqrAnalysis
!end interface cusolverSpXcsrqrAnalysis
! cusolverSpXcsrqrAnalysisBatched
interface cusolverSpXcsrqrAnalysisBatched
integer(c_int) function cusolverSpXcsrqrAnalysisBatched( cusolver_Hndl, m, n, nnzA, &
descrA, csrRowPtrA, csrColIndA, info ) &
bind(C,name="cusolverSpXcsrqrAnalysisBatched")
use iso_c_binding
import cusolverSpHandle, cusparseMatDescr, csrqrInfo
type(cusolverSpHandle), value :: cusolver_Hndl
type(cusparseMatDescr), value :: descrA
type(csrqrInfo), value :: info
integer(c_int), value :: m, n, nnzA
!! !pgi$ ignore_tkr (d) csrRowPtrA, (d) csrColIndA
integer(c_int),device :: csrRowPtrA(*), csrColIndA(*)
end function cusolverSpXcsrqrAnalysisBatched
end interface cusolverSpXcsrqrAnalysisBatched
! cusolverSpDcsrqrsvBatched
interface cusolverSpDcsrqrsvBatched
integer(c_int) function cusolverSpDcsrqrsvBatched( cusolver_Hndl, m, n, nnzA, &
descrA, csrValA, csrRowPtrA, csrColIndA, b, x, BatchSize, info, pBuffer ) &
bind(C,name="cusolverSpDcsrqrsvBatched")
use iso_c_binding
import cusolverSpHandle, cusparseMatDescr, csrqrInfo
type(cusolverSpHandle), value :: cusolver_Hndl
type(cusparseMatDescr), value :: descrA
type(csrqrInfo), value :: info
integer(c_int), value :: m, n, nnzA, Batchsize
real(c_double), device :: csrValA(*)
!! !pgi$ ignore_tkr (d) csrRowPtrA, (d) csrColIndA
integer(c_int),device :: csrRowPtrA(*), csrColIndA(*)
real(c_double), device :: b(*), x(*)
!! !pgi$ ignore_tkr pBuffer
character(c_char), device :: pBuffer(*)
end function cusolverSpDcsrqrsvBatched
end interface cusolverSpDcsrqrsvBatched
end module cuSOLVER_SP_QR
次に、Batched Sparse QR 分解を行うためのドライバルーチン(main.cuf) を例示する。なお、行列の CSR フォーマットは、row-major order で入力される。
Batched Sparse QR メインプログラム（main.cuf）
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved. ! ! Permission is hereby granted, free of charge, to any person obtaining a ! copy of this software and associated documentation files (the "Software"), ! to deal in the Software without restriction, including without limitation ! the rights to use, copy, modify, merge, publish, distribute, sublicense, ! and/or sell copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER ! DEALINGS IN THE SOFTWARE. ! program main use cudafor use cusparse use cusolver_SP_QR implicit none integer, parameter :: idbg = 2 ! debug write enable integer, parameter :: m = 4, nnzA = 7 integer, parameter :: batchsize = 4 integer :: cusolver_status integer :: status1, status2, status3, status4, status5 type(cusolverSpHandle) :: cusolver_Hndl type(cusparseMatDescr) :: descrA type(csrqrInfo) :: info integer :: csrRowPtrA(m+1), csrColIndA(nnzA) integer, device :: d_csrRowPtrA(m+1), d_csrColIndA(nnzA) double precision :: csrValA(nnzA), b(m) double precision, device:: d_csrValA(nnzA*batchsize), d_b(m*batchsize), d_x(m*batchsize) double precision :: csrValABatch(nnzA*batchsize), bBatch(m*batchsize), xbatch(m*batchsize) integer(8) :: size_internal, size_qr character(c_char), device, allocatable :: pBuffer(:) ! locals integer :: i, j, colidx, batchId integer :: ierr_code double precision :: eps, Areg, breg, xreg ! Result integer :: row, baseA, start, end, col double precision :: csrValAj double precision :: sup_res, r, Ax ! Randam Number double precision :: rnd ! | 1 | ! A = | 2 | ! | 3 | ! | 0.1 0.1 0.1 4 | ! ! CSR of A is based 1 indexing <==== caution (CUSPARSE_INDEX_BASE_ONE) ! ! b = [1 1 1 1] ! step 1: data setting csrRowPtrA(1:m+1) = (/1, 2, 3, 4, 8/) csrColIndA(1:nnzA) = (/1, 2, 3, 1, 2, 3, 4/) csrValA(1:nnzA) = (/1.0d0, 2.0d0, 3.0d0, 0.1d0, 0.1d0, 0.1d0, 4.0d0/) b(1:m) = (/1.0d0, 1.0d0, 1.0d0, 1.0d0/) print *, "sizeof(csrRowPtrA, d_csrRowPtrA, csrColIndA, d_csrColIndA)", & sizeof(csrRowPtrA(1)), sizeof(d_csrRowPtrA(1)), sizeof(csrColIndA(1)), sizeof(d_csrColIndA(1)) ! prepare Aj and bj on host do colidx = 1, nnzA Areg = csrValA(colidx) do batchId = 1, batchSize call random_number(rnd) eps = dble( mod(rnd,100.d0) + 1.0d0 ) * 1.0D-4 csrValABatch((batchId-1)*nnzA + colidx) = Areg + eps enddo enddo do j = 1, m breg = b(j) do batchId = 1, batchSize call random_number(rnd) eps = dble( mod(rnd,100.d0) + 1.0d0 ) * 1.0D-4 bBatch((batchId-1)*m + j) = breg + eps; enddo enddo if (idbg == 2) print *,"csrValABatch", csrValABatch if (idbg == 2) print *,"bBatch", bBatch ! step 2: create cusolver handle, qr info and matrix descriptor status1 = cusolverSpCreate(cusolver_Hndl) if (idbg == 1) print *, "status1 = ", status1 status2 = cusparseCreateMatDescr(descrA) if (idbg == 1) print *, "status2 = ", status2 status3 = cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL) if (idbg == 1) print *, "status3 = ", status3 status4 = cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ONE) ! base=1 if (idbg == 1) print *, "status4 = ", status4 status5 = cusolverSpCreateCsrqrInfo(info) if (idbg == 1) print *, "status5 = ", status5 ! step 3: copy Aj and bj to device d_csrValA = csrValABatch d_csrColIndA= csrColIndA d_csrRowPtrA= csrRowPtrA d_b = bBatch ! step 4: symbolic analysis cusolver_status = cusolverSpXcsrqrAnalysisBatched( & cusolver_Hndl, m, m, nnzA, descrA, d_csrRowPtrA, d_csrColIndA, info ) if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step4 cusolver_status = ", cusolver_status print *, 'error code = ' ,cusolver_status, cudaGetErrorString( cudaGetLastError() ) end if ! step 5: prepare working space cusolver_status = cusolverSpDcsrqrBufferInfoBatched( & cusolver_Hndl, m, m, nnzA, descrA, d_csrValA, d_csrRowPtrA, d_csrColIndA, & batchsize, info, size_internal, size_qr) if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step5 cusolver_status = ", cusolver_status end if print *, "numerical factorization needs internal data (bytes) =", size_internal print *, "numerical factorization needs working space (bytes) =", size_qr allocate(pBuffer(size_qr), stat=ierr_code) ! Bytes if (idbg == 1) print *, "step5: alloc ierr_code = ", ierr_code ! step 6: numerical factorization cusolver_status = cusolverSpDcsrqrsvBatched( & cusolver_Hndl, m, m, nnzA, & descrA, d_csrValA, d_csrRowPtrA, d_csrColIndA, & d_b, d_x, & batchSize, & info, & pBuffer) if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step6 cusolver_status = ", cusolver_status print *, 'error code = ' ,cusolver_status, cudaGetErrorString( cudaGetLastError() ) end if ! step 7: check residual ! xBatch = [x0, x1, x2, ...] ! copy into CPU mem xbatch = d_x if (idbg == 2) print *, "xbatch=",xbatch print *, "******* Fortran Index Base is adjusted : BaseA = 0 ********************* " baseA = 0 print *, " ============ Residual CHECK ===============" do batchId = 1 , batchsize ! measure |bj - Aj*Xj| sup_res = 0.d0 do row = 1, m start = csrRowPtrA(row ) - baseA end = csrRowPtrA(row+1) - baseA Ax = 0.d0 do colidx = start, end-1 col = csrColIndA(colidx) - baseA Areg = csrValAbatch((batchId-1)*nnzA + colidx) ! Areg = csrValAj(colidx) Xreg = xBatch((batchId-1)*m + col) ! Xreg = xj(col) Ax = Ax + Areg * Xreg end do r = bBatch((batchID-1)*m + row) - Ax ! r = bj(row) -Ax ! sup_res = (sup_res > fabs(r))? sup_res : fabs(r); sup_res = max( abs(r), sup_res ) ! if (sup_res > abs(r)) then ! sup_res = sup_res ! else ! sup_res = abs(r) ! end if end do print *, "batchId =", batchId," sup|bj - Aj*xj| =", sup_res end do print *, " ============ Result X ===============" do batchId = 1 , batchsize do row = 1, m print *, "x(", batchId, ")= [", row,"]", xBatch((batchId-1)*m + row) end do print *, "" end do end
以下の Makefile を copy & paste しただけでは make 実行エラーとなる。makefile 記述上の[ tab ]デリミタが空白に変換されているため、明示的に空白を[ tab ]に変更する必要がある。(以下、一例）
main : main.f90
[ tab ] pgf90 main.f90
Makefile（コンパイル方法）
CC = nvcc
CFLAGS =
FC = pgf90
FFLAGS = -O3 -Minfo -Mcuda=cuda8.0,cc35,cc60 -Mcudalib=cusparse
LDFLAGS = -Mcuda=cuda8.0,cc35,cc60
TARGET = SPsolve
INCLUDE =
LIBS = -lcusolver # -lgomp
# Ubuntu上では、以下のライブラリパスと-lgomp が必要となる場合あり
# NVIDIA cuSOLVER ライブラリは gcc の openmp ライブラリの一部を必要とする
#LDIR = -L/usr/lib/gcc/x86_64-linux-gnu/5
modfile = cusolver_mod.cuf
SRC1 = main.cuf
# サフィックスルール適用対象の拡張子の定義
.SUFFIXES: .cuf .cu .o
MODOBJ = $(modfile:.cuf=.o)
OBJ1 = $(SRC1:.cuf=.o)
MOD1 = cusolver_sp_qr.mod
RUN : $(TARGET)
@echo "===== Program run ====="
$(TARGET)
$(TARGET): $(MODOBJ) $(OBJ1)
$(FC) -o $@ $^ $(LDFLAGS) $(LDIR) $(LIBS)
$(MOD1) : $(modfile)
$(FC) $(FFLAGS) -c $?
$(OBJ1): $(MOD1) $(SRC1)
.cuf.o :
$(FC) -c $(FFLAGS) $<
.cu.o :
$(CC) -c $(CFLAGS) $<
.PHONY: clean
clean:
@echo 'Cleaning up...'
rm -f *.o *.mod $(TARGET) *~
コンパイル＆実行結果
[kato@photon32 QRtest1]$ make
pgf90 -c -O3 -Minfo -Mcuda=cuda8.0,cc35,cc60 -Mcudalib=cusparse cusolver_mod.cuf
pgf90 -c -O3 -Minfo -Mcuda=cuda8.0,cc35,cc60 -Mcudalib=cusparse main.cuf
main:
71, Loop unrolled 5 times (completely unrolled)
72, Array assignment / Forall at line 73 fused
Loop unrolled 7 times (completely unrolled)
Loop not vectorized: loop count too small
74, Loop unrolled 4 times (completely unrolled)
82, Loop not vectorized/parallelized: contains call
FMA (fused multiply-add) instruction(s) generated
91, Loop not vectorized/parallelized: contains call
FMA (fused multiply-add) instruction(s) generated
172, Loop not vectorized/parallelized: contains call
175, Outer loop unrolled 4 times (completely unrolled)
180, Unrolled inner loop 4 times
Generated 2 prefetch instructions for the loop
Unrolled inner loop 4 times
Generated 2 prefetch instructions for the loop
Unrolled inner loop 4 times
Generated 2 prefetch instructions for the loop
Unrolled inner loop 4 times
Generated 2 prefetch instructions for the loop
FMA (fused multiply-add) instruction(s) generated
203, Loop not vectorized/parallelized: contains call
pgf90 -o SPsolve cusolver_mod.o main.o -Mcuda=cuda8.0,cc35,cc60 -lcusolver
===== Program run =====
SPsolve
sizeof(csrRowPtrA, d_csrRowPtrA, csrColIndA, d_csrColIndA) 4
4 4 4
csrValABatch 1.000190792295671 2.000179731460909
3.000132064769846 0.1001762263696543 0.1001628016430586
0.1001500500225553 4.000169009995976 1.000119069206783
2.000163682999292 3.000144817609972 0.1001437479734487
0.1001670186653253 0.1001425331039673 4.000182114792401
1.000106716529557 2.000158878275024 3.000165795016470
0.1001463869055930 0.1001628171802322 0.1001307016646178
4.000187350713939 1.000180008450829 2.000115906560827
3.000160754590159 0.1001700415717258 0.1001531034351497
0.1001216954552562 4.000196496681088
bBatch 1.000182450045416 1.000166538252814
1.000186841054682 1.000110375401338 1.000145236365631
1.000145251688474 1.000116588275369 1.000155851050659
1.000125862765850 1.000112255030956 1.000187784792441
1.000198703067770 1.000133737619532 1.000188679946082
1.000142956695574 1.000175135506380
numerical factorization needs internal data (bytes) = 864
numerical factorization needs working space (bytes) = 246528
xbatch= 0.9999916593410789 0.5000383328163729
0.3333809377259567 0.2041067826368187 1.000026164043516
0.5000317024998334 0.3333561041470322 0.2041256247608860
1.000019144193291 0.5000164066433919 0.3333775041545496
0.2041369446005358 0.9999537374963444 0.5000653595450339
0.3333631223478223 0.2041274097506157
******* Fortran Index Base is adjusted : BaseA = 0 *********************
============ Residual CHECK ===============
batchId = 1 sup|bj - Aj*xj| = 8.8817841970012523E-016
batchId = 2 sup|bj - Aj*xj| = 6.6613381477509392E-016
batchId = 3 sup|bj - Aj*xj| = 1.1102230246251565E-015
batchId = 4 sup|bj - Aj*xj| = 4.4408920985006262E-016
============ Result X ===============
x( 1 )= [ 1 ] 0.9999916593410789
x( 1 )= [ 2 ] 0.5000383328163729
x( 1 )= [ 3 ] 0.3333809377259567
x( 1 )= [ 4 ] 0.2041067826368187
x( 2 )= [ 1 ] 1.000026164043516
x( 2 )= [ 2 ] 0.5000317024998334
x( 2 )= [ 3 ] 0.3333561041470322
x( 2 )= [ 4 ] 0.2041256247608860
x( 3 )= [ 1 ] 1.000019144193291
x( 3 )= [ 2 ] 0.5000164066433919
x( 3 )= [ 3 ] 0.3333775041545496
x( 3 )= [ 4 ] 0.2041369446005358
x( 4 )= [ 1 ] 0.9999537374963444
x( 4 )= [ 2 ] 0.5000653595450339
x( 4 )= [ 3 ] 0.3333631223478223
x( 4 )= [ 4 ] 0.2041274097506157
Batched Sparse QR example with OpenACC Fortran
cuSOLVER SP(スパース）系のライブラリへの Fortran Module は、上記で使用した cusolver_mod.cuf をそのまま使用するが、ファイルの suffix を cusolver_mod.f90 に変更していただきたい。OpenACC を使用して Batched Sparse QR 分解を行うためのドライバルーチン(main.f90) を以下に例示する。
OpenACC のディレクティブで留意すべきことは、CUDA Fortran or CUDA C で記述された（ライブラリ）ルーチンをコールする際に渡す実引数が「デバイスポインタ」である場合、これをコンパイルに伝えるために、host_data use_device() を使う。このディレクティブの使い方さえ誤りがなければ、既存のプログラムベースのものを OpenACC を使ってポーティングすることは大きな負担ではないはずだ。
Batched Sparse QR メインプログラム（main.f90）
! Copyright (C) HPC WORLD (Prometech Software, Inc. and GDEP Solutions, Inc.) All rights reserved. ! ! Permission is hereby granted, free of charge, to any person obtaining a ! copy of this software and associated documentation files (the "Software"), ! to deal in the Software without restriction, including without limitation ! the rights to use, copy, modify, merge, publish, distribute, sublicense, ! and/or sell copies of the Software, and to permit persons to whom the ! Software is furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in ! all copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL ! THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING ! FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER ! DEALINGS IN THE SOFTWARE. ! program main use cusparse use cusolver_SP_QR implicit none integer, parameter :: idbg = 2 ! debug write enable integer, parameter :: m = 4, nnzA = 7 integer, parameter :: batchsize = 4 integer :: cusolver_status integer :: status1, status2, status3, status4, status5 type(cusolverSpHandle) :: cusolver_Hndl type(cusparseMatDescr) :: descrA type(csrqrInfo) :: info integer :: csrRowPtrA(m+1), csrColIndA(nnzA) double precision :: csrValA(nnzA), b(m) double precision :: csrValABatch(nnzA*batchsize), bBatch(m*batchsize), xbatch(m*batchsize) integer(8) :: size_internal, size_qr character(c_char), allocatable :: pBuffer(:) ! locals integer :: i, j, colidx, batchId integer :: ierr_code double precision :: eps, Areg, breg, xreg ! Result integer :: row, baseA, start, end, col double precision :: csrValAj double precision :: sup_res, r, Ax ! Randam Number double precision :: rnd ! | 1 | ! A = | 2 | ! | 3 | ! | 0.1 0.1 0.1 4 | ! ! CSR of A is based 1 indexing <==== caution (CUSPARSE_INDEX_BASE_ONE) ! ! b = [1 1 1 1] ! step 1: data setting csrRowPtrA(1:m+1) = (/1, 2, 3, 4, 8/) csrColIndA(1:nnzA) = (/1, 2, 3, 1, 2, 3, 4/) csrValA(1:nnzA) = (/1.0d0, 2.0d0, 3.0d0, 0.1d0, 0.1d0, 0.1d0, 4.0d0/) b(1:m) = (/1.0d0, 1.0d0, 1.0d0, 1.0d0/) print *, "sizeof(csrRowPtrA, csrRowPtrA, csrColIndA, csrColIndA)", & sizeof(csrRowPtrA(1)), sizeof(csrRowPtrA(1)), sizeof(csrColIndA(1)), sizeof(csrColIndA(1)) !!call RANDOM_SEED() ! prepare Aj and bj on host do colidx = 1, nnzA Areg = csrValA(colidx) do batchId = 1, batchSize call random_number(rnd) eps = dble( mod(rnd,100.d0) + 1.0d0 ) * 1.0D-4 csrValABatch((batchId-1)*nnzA + colidx) = Areg + eps enddo enddo do j = 1, m breg = b(j) do batchId = 1, batchSize call random_number(rnd) eps = dble( mod(rnd,100.d0) + 1.0d0 ) * 1.0D-4 bBatch((batchId-1)*m + j) = breg + eps; enddo enddo xbatch = 0.d0 ! initilize if (idbg == 2) print *,"csrValABatch", csrValABatch if (idbg == 2) print *,"bBatch", bBatch ! step 2: create cusolver handle, qr info and matrix descriptor status1 = cusolverSpCreate(cusolver_Hndl) if (idbg == 1) print *, "status1 = ", status1 status2 = cusparseCreateMatDescr(descrA) if (idbg == 1) print *, "status2 = ", status2 status3 = cusparseSetMatType(descrA, CUSPARSE_MATRIX_TYPE_GENERAL) if (idbg == 1) print *, "status3 = ", status3 status4 = cusparseSetMatIndexBase(descrA, CUSPARSE_INDEX_BASE_ONE) ! base=1 if (idbg == 1) print *, "status4 = ", status4 status5 = cusolverSpCreateCsrqrInfo(info) if (idbg == 1) print *, "status5 = ", status5 ! step 3: copy Aj and bj to device ! d_csrValA = csrValABatch ! d_csrColIndA= csrColIndA ! d_csrRowPtrA= csrRowPtrA ! d_b = bBatch ! step 4: symbolic analysis !$acc data copyin(csrValABatch, csrColIndA, csrRowPtrA) copyin(bBatch) copyout(xbatch) !$acc host_data use_device(csrRowPtrA, csrColIndA) cusolver_status = cusolverSpXcsrqrAnalysisBatched( & cusolver_Hndl, m, m, nnzA, descrA, csrRowPtrA, csrColIndA, info ) !$acc end host_data if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step4 cusolver_status = ", cusolver_status end if ! step 5: prepare working space !$acc host_data use_device(csrValABatch, csrRowPtrA, csrColIndA) cusolver_status = cusolverSpDcsrqrBufferInfoBatched( & cusolver_Hndl, m, m, nnzA, descrA, csrValABatch, csrRowPtrA, csrColIndA, & batchsize, info, size_internal, size_qr) !$acc end host_data if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step5 cusolver_status = ", cusolver_status end if print *, "numerical factorization needs internal data (bytes) =", size_internal print *, "numerical factorization needs working space (bytes) =", size_qr allocate(pBuffer(size_qr), stat=ierr_code) ! Bytes if (idbg == 1) print *, "step5: alloc ierr_code = ", ierr_code !$acc enter data copyin(pBuffer) ! step 6: numerical factorization !$acc host_data use_device(csrValABatch, csrRowPtrA, csrColIndA, bBatch, xbatch, pBuffer) cusolver_status = cusolverSpDcsrqrsvBatched( & cusolver_Hndl, m, m, nnzA, & descrA, csrValABatch, csrRowPtrA, csrColIndA, & bBatch, xbatch, & batchSize, & info, & pBuffer) !$acc end host_data !$acc exit data delete(pBuffer) !$acc end data if (idbg == 1) then if (cusolver_status /= CUSOLVER_STATUS_SUCCESS) print *, "step6 cusolver_status = ", cusolver_status end if ! step 7: check residual ! xBatch = [x0, x1, x2, ...] ! copy into CPU mem if (idbg == 2) print *, "xbatch=",xbatch print *, "******* Fortran Index Base is adjusted : BaseA = 0 ********************* " baseA = 0 print *, " ============ Residual CHECK ===============" do batchId = 1 , batchsize ! measure |bj - Aj*Xj| sup_res = 0.d0 do row = 1, m start = csrRowPtrA(row ) - baseA end = csrRowPtrA(row+1) - baseA Ax = 0.d0 do colidx = start, end-1 col = csrColIndA(colidx) - baseA Areg = csrValAbatch((batchId-1)*nnzA + colidx) ! Areg = csrValAj(colidx) Xreg = xBatch((batchId-1)*m + col) ! Xreg = xj(col) Ax = Ax + Areg * Xreg end do r = bBatch((batchID-1)*m + row) - Ax ! r = bj(row) -Ax ! sup_res = (sup_res > fabs(r))? sup_res : fabs(r); sup_res = max( abs(r), sup_res ) end do print *, "batchId =", batchId," sup|bj - Aj*xj| =", sup_res end do print *, " ============ Result X ===============" do batchId = 1 , batchsize do row = 1, m print *, "x(", batchId, ")= [", row,"]", xBatch((batchId-1)*m + row) end do print *, "" end do end
以下の Makefile を copy & paste しただけでは make 実行エラーとなる。makefile 記述上の[ tab ]デリミタが空白に変換されているため、明示的に空白を[ tab ]に変更する必要がある。(以下、一例）以下の makefile の中で -Mcuda オプションは、CUDA システムライブラリ・リストをリンカーに指示するために必ず指定して下さい。
main : main.f90
[ tab ] pgf90 main.f90
Makefile（コンパイル方法）
CC = nvcc
CFLAGS =
FC = pgf90
FFLAGS = -O2 -Minfo -acc -ta=tesla,cuda8.0,cc35,cc60 -Mcudalib=cusparse -Mcuda
LDFLAGS = -acc -ta=tesla,cuda8.0,cc35,cc60 -Mcuda
TARGET = SPsolve
INCLUDE =
LIBS = -lcusolver # -lgomp
# Ubuntu上では、以下のライブラリパスと-lgompが必要となる場合あり
#LDIR = -L/usr/lib/gcc/x86_64-linux-gnu/5
modfile = cusolver_mod.f90
SRC1 = main.f90
# サフィックスルール適用対象の拡張子の定義
.SUFFIXES: .f90 .cu .o
MODOBJ = $(modfile:.f90=.o)
OBJ1 = $(SRC1:.f90=.o)
MOD1 = cusolver_sp_qr.mod
RUN : $(TARGET)
@echo "===== Program run ====="
$(TARGET)
$(TARGET): $(MODOBJ) $(OBJ1)
$(FC) -o $@ $^ $(LDFLAGS) $(LDIR) $(LIBS)
$(MOD1) : $(modfile)
$(FC) $(FFLAGS) -c $?
$(OBJ1): $(MOD1) $(SRC1)
.f90.o :
$(FC) -c $(FFLAGS) $<
.cu.o :
$(CC) -c $(CFLAGS) $<
.PHONY: clean
clean:
@echo 'Cleaning up...'
rm -f *.o *.mod $(TARGET) *~
コンパイル＆実行結果
[kato@photon32 OpenACC]$ make
pgf90 -c -O2 -Minfo -acc -ta=tesla,cuda8.0,cc35,cc60 -Mcudalib=cusparse -Mcuda cusolver_mod.f90
pgf90 -c -O2 -Minfo -acc -ta=tesla,cuda8.0,cc35,cc60 -Mcudalib=cusparse -Mcuda main.f90
main:
67, Loop unrolled 5 times (completely unrolled)
68, Array assignment / Forall at line 69 fused
Loop unrolled 7 times (completely unrolled)
Loop not vectorized: loop count too small
70, Loop unrolled 4 times (completely unrolled)
79, Loop not vectorized/parallelized: contains call
FMA (fused multiply-add) instruction(s) generated
88, Loop not vectorized/parallelized: contains call
FMA (fused multiply-add) instruction(s) generated
95, Loop unrolled 16 times (completely unrolled)
122, Generating copyin(csrrowptra(:),bbatch(:),csrcolinda(:))
Generating copyout(xbatch(:))
Generating copyin(csrvalabatch(:))
150, Generating enter data copyin(pbuffer(:))
163, Generating exit data delete(pbuffer(:))
182, Loop not vectorized/parallelized: contains call
185, Outer loop unrolled 4 times (completely unrolled)
190, Unrolled inner loop 4 times
Generated 2 prefetch instructions for the loop
Unrolled inner loop 4 times
Generated 2 prefetch instructions for the loop
Unrolled inner loop 4 times
Generated 2 prefetch instructions for the loop
Unrolled inner loop 4 times
Generated 2 prefetch instructions for the loop
FMA (fused multiply-add) instruction(s) generated
208, Loop not vectorized/parallelized: contains call
pgf90 -o SPsolve cusolver_mod.o main.o -acc -ta=tesla,cuda8.0,cc35,cc60 -Mcuda -lcusolver
===== Program run =====
SPsolve
sizeof(csrRowPtrA, csrRowPtrA, csrColIndA, csrColIndA) 4
4 4 4
csrValABatch 1.000190792295671 2.000179731460909
3.000132064769846 0.1001762263696543 0.1001628016430586
0.1001500500225553 4.000169009995976 1.000119069206783
2.000163682999292 3.000144817609972 0.1001437479734487
0.1001670186653253 0.1001425331039673 4.000182114792401
1.000106716529557 2.000158878275024 3.000165795016470
0.1001463869055930 0.1001628171802322 0.1001307016646178
4.000187350713939 1.000180008450829 2.000115906560827
3.000160754590159 0.1001700415717258 0.1001531034351497
0.1001216954552562 4.000196496681088
bBatch 1.000182450045416 1.000166538252814
1.000186841054682 1.000110375401338 1.000145236365631
1.000145251688474 1.000116588275369 1.000155851050659
1.000125862765850 1.000112255030956 1.000187784792441
1.000198703067770 1.000133737619532 1.000188679946082
1.000142956695574 1.000175135506380
numerical factorization needs internal data (bytes) = 864
numerical factorization needs working space (bytes) = 246528
xbatch= 0.9999916593410789 0.5000383328163729
0.3333809377259567 0.2041067826368187 1.000026164043516
0.5000317024998334 0.3333561041470322 0.2041256247608860
1.000019144193291 0.5000164066433919 0.3333775041545496
0.2041369446005358 0.9999537374963444 0.5000653595450339
0.3333631223478223 0.2041274097506157
******* Fortran Index Base is adjusted : BaseA = 0 *********************
============ Residual CHECK ===============
batchId = 1 sup|bj - Aj*xj| = 8.8817841970012523E-016
batchId = 2 sup|bj - Aj*xj| = 6.6613381477509392E-016
batchId = 3 sup|bj - Aj*xj| = 1.1102230246251565E-015
batchId = 4 sup|bj - Aj*xj| = 4.4408920985006262E-016
============ Result X ===============
x( 1 )= [ 1 ] 0.9999916593410789
x( 1 )= [ 2 ] 0.5000383328163729
x( 1 )= [ 3 ] 0.3333809377259567
x( 1 )= [ 4 ] 0.2041067826368187
x( 2 )= [ 1 ] 1.000026164043516
x( 2 )= [ 2 ] 0.5000317024998334
x( 2 )= [ 3 ] 0.3333561041470322
x( 2 )= [ 4 ] 0.2041256247608860
x( 3 )= [ 1 ] 1.000019144193291
x( 3 )= [ 2 ] 0.5000164066433919
x( 3 )= [ 3 ] 0.3333775041545496
x( 3 )= [ 4 ] 0.2041369446005358
x( 4 )= [ 1 ] 0.9999537374963444
x( 4 )= [ 2 ] 0.5000653595450339
x( 4 )= [ 3 ] 0.3333631223478223
x( 4 )= [ 4 ] 0.2041274097506157
11. cuFFT example with OpenACC Fortran
OpenACC による cuFFT ライブラリを使用した 2 次元 FFT のプログラムの例を以下に示す。acc data データ領域の指定と acc host_data use_device でデバイスポインタを使用する変数を指定するだけで、FFT の処理が可能となる。
cuFFT example（cuFFT.f90）
program cufft2dTest use cufft use openacc integer, parameter :: m=768, n=512 complex, allocatable :: a(:,:),b(:,:),c(:,:) real, allocatable :: r(:,:),q(:,:) integer :: iplan1, iplan2, iplan3, ierr allocate(a(m,n),b(m,n),c(m,n)) allocate(r(m,n),q(m,n)) a = 1; r = 1 xmx = -99.0 ierr = cufftPlan2D(iplan1,m,n,CUFFT_C2C) !$acc data copyin(a) create(b) copyout(c) !$acc host_data use_device(a,b,c) ierr = ierr + cufftExecC2C(iplan1,a,b,CUFFT_FORWARD) ierr = ierr + cufftExecC2C(iplan1,b,c,CUFFT_INVERSE) !$acc end host_data ! scale c !$acc kernels c = c / (m*n) !$acc end kernels !$acc end data ! Check forward answer write(*,*) 'Max error C2C FWD: ', cmplx(maxval(real(b)) - sum(real(b)), & maxval(imag(b))) ! Check inverse answer write(*,*) 'Max error C2C INV: ', maxval(abs(a-c)) ! Real transform ierr = ierr + cufftPlan2D(iplan2,m,n,CUFFT_R2C) ierr = ierr + cufftPlan2D(iplan3,m,n,CUFFT_C2R) !$acc data copyin(r) copyout(q) create(b) !$acc host_data use_device(r,b,q) ierr = ierr + cufftExecR2C(iplan2,r,b) ierr = ierr + cufftExecC2R(iplan3,b,q) !$acc end host_data !$acc kernels xmx = maxval(abs(r-q/(m*n))) !$acc end kernels !$acc end data ! Check R2C + C2R answer write(*,*) 'Max error R2C/C2R: ', xmx ierr = ierr + cufftDestroy(iplan1) ierr = ierr + cufftDestroy(iplan2) ierr = ierr + cufftDestroy(iplan3) if (ierr.eq.0) then print *,"test PASSED" else print *,"test FAILED" endif end program cufft2dTest
コンパイルとリンクオプションには、必ず、-acc -Mcudalib=cufft を指定することが必要である。
コンパイル＆実行結果
[kato@photon32 Example]$ pgf90 -acc -Minfo=accel -O2 -ta=tesla,cc60,cc35,cuda8.0 -Mcudalib=cufft cuFFT.f90
cufft2dtest:
0, Accelerator kernel generated
Generating Tesla code
16, Generating copyin(a(:,:))
Generating create(b(:,:))
Generating copyout(c(:,:))
24, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
24, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
!$acc loop gang, vector(32) ! blockidx%x threadidx%x
38, Generating copyin(r(:,:))
Generating copyout(q(:,:))
Generating create(b(:,:))
45, Loop is parallelizable
Accelerator scalar kernel generated
Accelerator kernel generated
Generating Tesla code
45, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
!$acc loop gang, vector(32) ! blockidx%x threadidx%x
Generating implicit reduction(max:r$r)
[kato@photon32 Example]$ a.out
Max error C2C FWD: (0.000000,0.000000)
Max error C2C INV: 0.000000
Max error R2C/C2R: 0.000000
test PASSED
12. cuRAND example (乱数発生)
cuRAND example from OpenACC Host Code
OpenACC による cuRAND ライブラリを使用したプログラムの例（XORWOW乱数) を以下に示す。ホストプログラムから CUDA ライブラリを利用するため、デバイス側で使用するデータのみを OpenACC の acc data で指定し、acc host_data use_device でデバイス側で使用する変数を指定するだけで利用可能となる。なお、curand 関数を使用するルーチンでは必ず use curand （インタフェース）を指定することが必要である。
cuRAND test from OpenACC Host Code（cuRand.f90）
program testcurand ! compile with the flags -ta=tesla -Mcuda -Mcudalib=curand call cur1(1000, .true.); call cur1(1000, .false.) call cur2(1000, .true.); call cur2(1000, .false.) call cur3(1000, .true.); call cur3(1000, .false.) end ! subroutine cur1(n, onhost) use curand ! 必ず指定すること integer :: a(n) ! 整数 type(curandGenerator) :: g integer(8) nbits logical onhost, passing a = 0 passing = .true. if (onhost) then ! host 側のライブラリも提供している istat = curandCreateGeneratorHost(g,CURAND_RNG_PSEUDO_XORWOW) istat = curandGenerate(g, a, n) istat = curandDestroyGenerator(g) else ! GPU 側処理 !$acc data copy(a) istat = curandCreateGenerator(g,CURAND_RNG_PSEUDO_XORWOW) !$acc host_data use_device(a) istat = curandGenerate(g, a, n) !$acc end host_data istat = curandDestroyGenerator(g) !$acc end data endif nbits = 0 do i = 1, n if (i.lt.10) print *,i,a(i) nbits = nbits + popcnt(a(i)) end do print *,"Should be roughly half the bits set" nbits = nbits / n if ((nbits .lt. 12) .or. (nbits .gt. 20)) then passing = .false. else print *,"nbits is ",nbits," which passes" endif if (passing) then print *,"Test PASSED" else print *,"Test FAILED" endif end ! subroutine cur2(n, onhost) use curand real :: a(n) ! 単精度 type(curandGenerator) :: g logical onhost, passing a = 0.0 passing = .true. if (onhost) then istat = curandCreateGeneratorHost(g,CURAND_RNG_PSEUDO_XORWOW) istat = curandGenerate(g, a, n) istat = curandDestroyGenerator(g) else !$acc data copy(a) istat = curandCreateGenerator(g,CURAND_RNG_PSEUDO_XORWOW) !$acc host_data use_device(a) istat = curandGenerate(g, a, n) !$acc end host_data istat = curandDestroyGenerator(g) !$acc end data endif print *,"Should be uniform around 0.5" do i = 1, n if (i.lt.10) print *,i,a(i) if ((a(i).lt.0.0) .or. (a(i).gt.1.0)) passing = .false. end do rmean = sum(a)/n if ((rmean .lt. 0.4) .or. (rmean .gt. 0.6)) then passing = .false. else print *,"Mean is ",rmean," which passes" endif if (passing) then print *,"Test PASSED" else print *,"Test FAILED" endif end ! subroutine cur3(n, onhost) use curand real(8) :: a(n) ! 倍精度 type(curandGenerator) :: g logical onhost, passing a = 0.0d0 passing = .true. if (onhost) then istat = curandCreateGeneratorHost(g,CURAND_RNG_PSEUDO_XORWOW) istat = curandGenerate(g, a, n) istat = curandDestroyGenerator(g) else !$acc data copy(a) istat = curandCreateGenerator(g,CURAND_RNG_PSEUDO_XORWOW) !$acc host_data use_device(a) istat = curandGenerate(g, a, n) !$acc end host_data istat = curandDestroyGenerator(g) !$acc end data endif do i = 1, n if (i.lt.10) print *,i,a(i) if ((a(i).lt.0.0d0) .or. (a(i).gt.1.0d0)) passing = .false. end do rmean = sum(a)/n if ((rmean .lt. 0.4d0) .or. (rmean .gt. 0.6d0)) then passing = .false. else print *,"Mean is ",rmean," which passes" endif if (passing) then print *,"Test PASSED" else print *,"Test FAILED" endif end
コンパイルとリンクオプションには、必ず、-acc -Mcudalib=curand を指定することが必要である。また、システム CUDA ライブラリもリンクする必要があるため、-Mcuda オプションも必須となる。
コンパイル＆実行結果
[kato@photon32]$ pgf90 -acc -Minfo=accel -O2 -ta=tesla,cc60,cc35,cuda8.0 -Mcudalib=curand -Mcuda cuRand.f90
cur1:
21, Generating copy(a(:))
cur2:
60, Generating copy(a(:))
cur3:
98, Generating copy(a(:))
[kato@photon32 cuRAND]$ a.out
1 -1115749450
2 -339329097
3 167591721
4 -133303299
5 -321518802
6 1917059131
7 -1428853312
8 472148682
9 2019573489
Should be roughly half the bits set
nbits is 16 which passes
Test PASSED
1 -1115749450
2 -339329097
3 167591721
4 -133303299
5 -321518802
6 1917059131
7 -1428853312
8 472148682
9 2019573489
Should be roughly half the bits set
nbits is 16 which passes
Test PASSED
Should be uniform around 0.5
1 0.7402194
2 0.9209938
3 3.9020490E-02
4 0.9689629
5 0.9251406
6 0.4463501
7 0.6673192
8 0.1099307
9 0.4702186
Mean is 0.5014752 which passes
Test PASSED
Should be uniform around 0.5
1 0.7402194
2 0.9209938
3 3.9020490E-02
4 0.9689629
5 0.9251406
6 0.4463501
7 0.6673192
8 0.1099307
9 0.4702186
Mean is 0.5014752 which passes
Test PASSED
1 0.4384508447184235
2 0.4603647022031429
3 0.2502147080176417
4 0.4947437677616333
5 5.3011123716805220E-002
6 0.3376992580306317
7 0.3967625188337749
8 0.8744186605648221
9 0.4821668323147305
Mean is 0.5150273 which passes
Test PASSED
1 0.4384508447184235
2 0.4603647022031429
3 0.2502147080176417
4 0.4947437677616333
5 5.3011123716805220E-002
6 0.3376992580306317
7 0.3967625188337749
8 0.8744186605648221
9 0.4821668323147305
Mean is 0.5150273 which passes
Test PASSED
cuRAND example from OpenACC device code
OpenACC の計算領域で cuRAND のデバイスルーチン・ライブラリ関数（curand_uniform等）を使用する場合の例である。これは、ホストプログラム上で使用するものではなく、device 側専用のルーチンである。以下のプログラムでは、module mtests の中で、OpenACC によるデバイスコードを定義している。これは、グローバルサブルーチンとして他のルーチンから利用できる。ここでの技術的な留意点は、PGI Fortran 専用の cuRAND OpenACC 専用 MODULE インタフェース(openacc_curand) を use 文で指定する必要があるという点である。PGI は、cuRand 関数に対して openacc routine ディレクティブを指定した Fortran MODULE を提供している。
cuRAND test from OpenACC device Code（cuRand2.f90）
module mtests
integer, parameter :: n = 1000
contains
subroutine testrand( a, b )
use openacc_curand ! OpenACC用のcurand I/F 必ず指定する
real :: a(n), b(n)
type(curandStateXORWOW) :: h
integer(8) :: seed, seq, offset
!$acc parallel num_gangs(1) vector_length(1) copy(a,b) private(h)
seed = 12345 ! この時点でデバイスが側に処理が移る。 a, b 配列ともにデバイスへコピー
seq = 0
offset = 0
call curand_init(seed, seq, offset, h)
!$acc loop seq ! 以下のループは sequential 実行
do i = 1, n
a(i) = curand_uniform(h)
b(i) = curand_normal(h)
end do
!$acc end parallel ! この時点で a(n), b(n) の値はホスト側に戻される
return
end subroutine
end module mtests
program t
use mtests
real :: a(n), b(n), c(n)
logical passing
a = 1.0
b = 2.0
passing = .true.
call testrand(a,b) ! ホスト側から call する。引数a, b はホスト上の配列を渡す
c = a
print *,"Should be uniform around 0.5"
do i = 1, n
if (i.lt.10) print *,i,c(i)
if ((c(i).lt.0.0) .or. (c(i).gt.1.0)) passing = .false.
end do
rmean = sum(c)/n
if ((rmean .lt. 0.4) .or. (rmean .gt. 0.6)) then
passing = .false.
else
print *,"Mean is ",rmean," which passes"
endif
c = b
print *,"Should be normal around 0.0"
nc1 = 0;
nc2 = 0;
do i = 1, n
if (i.lt.10) print *,i,c(i)
if ((c(i) .gt. -4.0) .and. (c(i) .lt. 0.0)) nc1 = nc1 + 1
if ((c(i) .gt. 0.0) .and. (c(i) .lt. 4.0)) nc2 = nc2 + 1
end do
print *,"Found on each side of zero ",nc1,nc2
if (abs(nc1-nc2) .gt. (n/10)) npassing = .false.
rmean = sum(c,mask=abs(c).lt.4.0)/n
if ((rmean .lt. -0.1) .or. (rmean .gt. 0.1)) then
passing = .false.
else
print *,"Mean is ",rmean," which passes"
endif
if (passing) then
print *,"Test PASSED"
else
print *,"Test FAILED"
endif
end program
コンパイルとリンクオプションには、-acc -Mcudalib=cusparse -Mcuda のオプションの指定が必要である。OpenACC 領域内で cuRAND ライブラリを利用する場合は、PGI Fortran OpenACC 専用の cuRAND のデバイスルーチン・ライブラリをリンクする必要があるため、リンクオプションに注意が必要である。リンク時に NVIDIA LLVM 配下のリンク環境を使うのではなく、従来のリンケージ処理を行うために、-ta=tesla,nollvm サブオプションの指定が必要となる。
コンパイル＆実行結果
[kato@photon32 cuRAND]$ pgf90 -acc -Minfo=accel -O2 -ta=tesla,cc60,cc35,cuda8.0,nollvm cuRand2.f90
testrand:
10, Generating copy(b(:),a(:))
Accelerator kernel generated
Generating Tesla code
16, !$acc loop seq
10, CUDA shared memory used for h
[kato@photon32 cuRAND]$ a.out
Should be uniform around 0.5
1 0.2988985
2 0.4019704
3 0.7425825
4 0.7073491
5 0.5512256
6 0.4850157
7 9.7107515E-02
8 0.3596036
9 0.5777864
Mean is 0.4975990 which passes
Should be normal around 0.0
1 1.348672
2 -0.3331721
3 1.273486
4 -1.008033
5 -0.6210795
6 -0.5052772
7 1.598479
8 1.133891
9 -1.831225
Found on each side of zero 469 531
Mean is 4.2783763E-02 which passes
Test PASSEDD
13. cuSPARSE from OpenACC Host Code
ホスト用プログラムからデバイス側で処理する cuSPARSE ライブラリの利用例を示す。
cuSPARSE ルーチンの処理には直接関係ないが、以下の例では、cusparseSetStream を使って個々の cuSPARSE ライブラリルーチンによって使われる stream をセットしている。これを行うことで、対象となるハンドル h の cuSPARSE 関数は分離された stream で処理され、GPU 上で自動的なオーバラップ処理が可能となる。このアプローチはとりわけ、シングルタスクの実行が相対的に小さな場合で GPU 内の処理容量を満たさないような時に有効である。あるいは、計算とデータ転送が並行に行われるような場面で使用される。
cuSPARSE from OpenACC Host Code（cuSPARSE.f90）
program sparseMatVec
integer n
n = 25 ! # rows/cols in dense matrix
call sparseMatVecSub1(n)
n = 45 ! # rows/cols in dense matrix
call sparseMatVecSub1(n)
end program
subroutine sparseMatVecSub1(n)
use openacc
use cusparse
implicit none
integer n
! dense data
real(4), allocatable :: Ade(:,:), x(:), y(:)
! sparse CSR arrays
real(4), allocatable :: csrValA(:)
integer, allocatable :: nnzPerRowA(:), csrRowPtrA(:), csrColIndA(:)
allocate(Ade(n,n), x(n), y(n))
allocate(csrValA(n))
allocate(nnzPerRowA(n), csrRowPtrA(n+1), csrColIndA(n))
call sparseMatVecSub2(Ade, x, y, csrValA, nnzPerRowA, csrRowPtrA, &
csrColIndA, n)
deallocate(Ade)
deallocate(x)
deallocate(y)
deallocate(csrValA)
deallocate(nnzPerRowA)
deallocate(csrRowPtrA)
deallocate(csrColIndA)
end subroutine
subroutine sparseMatVecSub2(Ade, x, y, csrValA, nnzPerRowA, csrRowPtrA, &
csrColIndA, n)
use openacc
use cusparse
implicit none
! dense data
real(4) :: Ade(n,n), x(n), y(n)
! sparse CSR arrays
real(4) :: csrValA(n)
integer :: nnzPerRowA(n), csrRowPtrA(n+1), csrColIndA(n)
integer :: n, nnz, status, i
type(cusparseHandle) :: h
type(cusparseMatDescr) :: descrA
! parameters
real(4) :: alpha, beta
! result
real(4) :: xerr
! initalize CUSPARSE and matrix descriptor
status = cusparseCreate(h)
if (status /= CUSPARSE_STATUS_SUCCESS) &
write(*,*) 'cusparseCreate error: ', status
status = cusparseCreateMatDescr(descrA)
status = cusparseSetMatType(descrA, &
CUSPARSE_MATRIX_TYPE_GENERAL)
status = cusparseSetMatIndexBase(descrA, &
CUSPARSE_INDEX_BASE_ONE)
status = cusparseSetStream(h, acc_get_cuda_stream(acc_async_sync)) ! cusparseが利用する一つのcuda streamを生成
!$acc data create(Ade, x, y, csrValA, nnzPerRowA, csrRowPtrA, csrColIndA)
! Initialize matrix (upper circular shift matrix) ! Circulant matrix
!$acc kernels
Ade = 0.0
do i = 1, n-1
Ade(i,i+1) = 1.0
end do
Ade(n,1) = 1.0
! Initialize vectors and constants
do i = 1, n
x(i) = i
enddo
y = 0.0
!$acc end kernels
!$acc update host(x)
write(*,*) 'Original vector:'
write(*,'(5(1x,f7.2))') x
! convert matrix from dense to csr format
!$acc host_data use_device(Ade, nnzPerRowA, csrValA, csrRowPtrA, csrColIndA)
status = cusparseSnnz_v2(h, CUSPARSE_DIRECTION_ROW, &
n, n, descrA, Ade, n, nnzPerRowA, nnz)
status = cusparseSdense2csr(h, n, n, descrA, Ade, n, &
nnzPerRowA, csrValA, csrRowPtrA, csrColIndA)
!$acc end host_data
! A is upper circular shift matrix
! y = alpha*A*x + beta*y
alpha = 1.0
beta = 0.0
!$acc host_data use_device(csrValA, csrRowPtrA, csrColIndA, x, y)
status = cusparseScsrmv(h, CUSPARSE_OPERATION_NON_TRANSPOSE, &
n, n, n, alpha, descrA, csrValA, csrRowPtrA, &
csrColIndA, x, beta, y)
!$acc end host_data
!$acc wait
write(*,*) 'Shifted vector:'
write(*,'(5(1x,f7.2))') y
! shift-down y and add original x
! A' is lower circular shift matrix
! x = alpha*A'*y + beta*x
beta = -1.0
!$acc host_data use_device(csrValA, csrRowPtrA, csrColIndA, x, y)
status = cusparseScsrmv(h, CUSPARSE_OPERATION_TRANSPOSE, &
n, n, n, alpha, descrA, csrValA, csrRowPtrA, &
csrColIndA, y, beta, x)
!$acc end host_data
!$acc kernels
xerr = maxval(abs(x))
!$acc end kernels
!$acc end data
write(*,*) 'Max error = ', xerr
if (xerr.le.1.e-5) then
write(*,*) 'Test PASSED'
else
write(*,*) 'Test FAILED'
endif
end subroutine
コンパイルとリンクオプションには、-acc -Mcudalib=cusparse -Mcuda が必要となる。
コンパイル＆実行結果
[kato@photon32 cuSPARSE]$ pgf90 -Minfo=accel -acc -O2 -ta=tesla,cc60,cuda8.0 -Mcudalib=cusparse cuSparse.f90 -Mcuda
sparsematvecsub2:
0, Accelerator kernel generated
Generating Tesla code
74, Generating create(nnzperrowa(:),x(:),y(:),csrrowptra(:),csrvala(:),csrcolinda(:),ade(:,:))
78, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
78, !$acc loop gang, vector(4) ! blockidx%y threadidx%y
!$acc loop gang, vector(32) ! blockidx%x threadidx%x
79, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
79, !$acc loop gang, vector(128) ! blockidx%x threadidx%x
81, Accelerator scalar kernel generated
85, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
85, !$acc loop gang, vector(128) ! blockidx%x threadidx%x
88, Loop is parallelizable
Accelerator kernel generated
Generating Tesla code
88, !$acc loop gang, vector(128) ! blockidx%x threadidx%x
91, Generating update self(x(:))
128, Loop is parallelizable
Accelerator scalar kernel generated
Accelerator kernel generated
Generating Tesla code
128, !$acc loop gang, vector(128) ! blockidx%x threadidx%x
Generating implicit reduction(max:x$r)
[kato@photon32 cuSPARSE]$ a.out
Original vector:
1.00 2.00 3.00 4.00 5.00
6.00 7.00 8.00 9.00 10.00
11.00 12.00 13.00 14.00 15.00
16.00 17.00 18.00 19.00 20.00
21.00 22.00 23.00 24.00 25.00
Shifted vector:
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
Max error = 0.000000
Test PASSED
Original vector:
1.00 2.00 3.00 4.00 5.00
6.00 7.00 8.00 9.00 10.00
11.00 12.00 13.00 14.00 15.00
16.00 17.00 18.00 19.00 20.00
21.00 22.00 23.00 24.00 25.00
26.00 27.00 28.00 29.00 30.00
31.00 32.00 33.00 34.00 35.00
36.00 37.00 38.00 39.00 40.00
41.00 42.00 43.00 44.00 45.00
Shifted vector:
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00
Max error = 0.000000
Test PASSED