OpenBLAS project sourceに基づいた BLAS/LAPACKライブラリ

対象 OpenBLAS BLAS LAPACK PGIコンパイラ オプション

 PGI 16.1 以降から OpenBLAS project sourceに基づいた最適化 BLAS/LAPACK を提供しました。これは、予め PGI コンパイラを用いてプリコンパイルされているライブラリで、-lblas -llapack リンク・オプションを付けてリンクできます。また、マルチスレッド・モードでも動作します。この方法を以下に説明します。なお、PGI 15.10 までのバージョンでは ACML ライブラリ内の BLAS/LAPACK をご利用ください。

PGI 製品にバンドルされた BLAS/LAPACK の使用法(PGI 16.1 以降)

 PGIコンパイルシステムから参照可能な OpenBLAS ライブラリは、PGIコンパイラでビルドされたものです。 OpenBLAS project source に基づいた性能最適化が施されております。シングルスレッド、マルチスレッドのどちらにも対応できます。以下のリンクオプションを指定して、当該 BLAS あるいは LAPACK ルーチンをリンクしてください。

PGI 16.1 以降の BLAS/LAPACK ライブラリをリンクする方法

● Linux/Windows/OS X 版シングルスレッド

【Fortran】
  pgf77 -fastsse -Minfo test.f -lblas -llapack
  pgfortran -fastsse -Minfo test.f90 -lblas -llapack
【C/C++】
  pgcc -fastsse -Minfo test.c -lblas -llapack -pgf90libs
  pgc++ -fastsse -Minfo test.cpp -lblas -llapack -pgf90libs (Windows版には C++ コマンドはありません)
  ※C/C++ コンパイラ使用時は-pgf90libsオプションを加えて下さい。

● OpenMPマルチスレッド並列用コンパイル (Windows版の OpenBLAS は、マルチスレッド対応ではありません)

【Fortran】 (Linux / OS X only)
  pgf77 -fastsse -Minfo test.f -lblas -llapack -mp
  pgfortran -fastsse -Minfo test.f90 -lblas -llapack -mp
【C/C++】 (Linux / OS X only)
  pgcc -fastsse -Minfo test.c -lblas -llapack -pgf90libs -mp
  pgc++ -fastsse -Minfo test.cpp -lblas -llapack -pgf90libs -mp
  ※C/C++ コンパイラ使用時は-pgf90libsオプションを加えて下さい。
  

Fortranでの使用例 (matrix multiply BLASを使用)

● ソースプログラム例

      program matmul_time
c
      integer i, j, k
      integer size, ntimes, m, n, p
      parameter (size=1000)
      parameter (m=size,n=size,p=size)
      parameter (ntimes=20)
      real*8 a, b, c, arow
      dimension a(m,n), b(n,p), c(n,p), arow(n)

      integer l
      real*8 walltime, mflops
      integer hz, clock0, clock1, clock11, clock2
      integer t

      real*8 second
      real*8 time0, time1, time11, time2

#if defined (OMPTIME)
      real*8 omp_get_wtime, omp_get_wtick
      print *, "OMP_Wtime resolution: ", omp_get_wtick()
#endif

#if defined (SYSCLK)
      call system_clock(count_rate=hz)
      print *, "system_clock resolution: ", real(1.d0/hz)
#endif
c
      call random_number(a)
      call random_number(b)

#if defined (SYSCLK)
! Elapsed time : Using F90 system_clock : Milli second
      call system_clock(count=clock0)

#elif defined (OMPTIME)
! Elapsed time : OpenMP build-in function: Micro second
      time0 = omp_get_wtime()
#else
! Elapsed time : SECOND() C function    : Micro second
      time0 = second()
#endif

      do l = 1, ntimes

#if defined MATMUL
! F90 Intrinsic function (MATMUL)
         c = matmul(a, b)

#elif defined BLAS
! AMD ACML(LAPACK) Library
         call dgemm('N','N', m, p, n, 1.0d0, a, m, b, n, 0.0d0, c, n)

#else
! Source coding
!$omp parallel default(shared)
!$omp do
         do j = 1, p
            do i = 1, m
               c(i,j) = 0.0
            enddo
         enddo
         do i = 1, m
!$omp do
            do ii = 1, n
               arow(ii) = a(i,ii)
            enddo

!$omp do schedule(runtime)
            do j = 1, p
               do k = 1, n
                  c(i,j) = c(i,j) + arow(k) * b(k,j)
               enddo
            enddo
         enddo
!$omp end parallel
#endif
         call dummy(c)
      enddo

#if defined (SYSCLK)
      call system_clock(count=clock1)
      call system_clock(count=clock11)
#elif defined (OMPTIME)
      time1 = omp_get_wtime()
      time11= omp_get_wtime()
#else
      time1 = second()
      time11= second()
#endif

      do i = 1, ntimes
         call dummy(c)
      enddo

#if defined (SYSCLK)
      call system_clock(count=clock2)
#elif defined (OMPTIME)
      time2 = omp_get_wtime()
#else
      time2 = second()
#endif

#if defined (SYSCLK)
      t = (clock1 - clock0) - (clock2 - clock11)
      walltime = real(t) / real(hz) / real(ntimes)
      mflops = (m*p*(2*n-1)) / (walltime * 1.0d+06)
      print *,"Elapsed time =", walltime , " (sec)"
#else
      walltime = ((time1 - time0) - (time2 - time11)) / real(ntimes)
      mflops = (m*p*(2*n-1)) / (walltime * 1.0d+06)
      print *,"Elapsed time =", walltime ," (sec)"
#endif
c
      print *, "M =",M,", N =",N,", P =",P
      print *, "MFLOPS = ", mflops
      print *, "c(10,10) = ", c(10,10)
c
      end
c
      subroutine dummy(c)
      return
      end

● コンパイル・リンク

BLAS/LAPACKルーチンを使用するため -DBLAS オプションを入れて、マルチスレッド(-mp)

$ pgf90 -fastsse -Minfo -DBLAS -DOMPTIME matmul.F -lblas -llapack -mp
matmul_time:
     29, Loop interchange produces reordered loop nest: 30,29
         2 loops fused
         Generated an alternate version of the loop
         Generated vector sse code for the loop
     30, Loop not fused: function call before adjacent loop
         2 loops fused
     34, Loop interchange produces reordered loop nest: 35,34
     52, Loop not vectorized/parallelized: contains call
    100, Loop not vectorized/parallelized: contains call   

● 実行

Intel(R) Core i5-4690(Haswell Refresh 第4世代)CPU @ 3.50GHz 4コアを使用

$ export OMP_NUM_THREADS=1
$ ./a.out
 OMP_Wtime resolution:    9.9999999999999995E-007
 Elapsed time =   4.1460651159286502E-002  (sec)
 M =         1000 , N =         1000 , P =         1000
 MFLOPS =     48214.38988789874  (48GFLOPS)
 c(10,10) =     244.1990574713752

$ export OMP_NUM_THREADS=2
$ ./a.out
 OMP_Wtime resolution:    9.9999999999999995E-007
 Elapsed time =   2.2754621505737305E-002  (sec)
 M =         1000 , N =         1000 , P =         1000
 MFLOPS =     87850.28568794150  (87GFLOPS)
 c(10,10) =     244.1990574713752
 
$ export OMP_NUM_THREADS=3
$ ./a.out
 OMP_Wtime resolution:    9.9999999999999995E-007
 Elapsed time =   1.7555427551269532E-002  (sec)
 M =         1000 , N =         1000 , P =         1000
 MFLOPS =     113867.9188622462  (113GFLOPS)
 c(10,10) =     244.1990574713752
 

$ export OMP_NUM_THREADS=4
$ ./a.out
 OMP_Wtime resolution:    9.9999999999999995E-007
 Elapsed time =   1.3460004329681398E-002  (sec)
 M =         1000 , N =         1000 , P =         1000
 MFLOPS =     148514.0681264043  (148GFLOPS)
 c(10,10) =     244.1990574713752

● BLASを使わず OpenMP プログラムで実行した場合(性能比較)

Intel(R) Core i5-4690(Haswell Refresh 第4世代)CPU @ 3.50GHz 4コアを使用

$ pgf90 -fastsse -Minfo -DOMPTIME matmul.F -mp -Mvect=simd:256 (-DBLASを除く)
matmul_time:
     44, Loop not vectorized/parallelized: too deeply nested
     56, Parallel region activated
     58, Parallel loop activated with static block schedule
     59, Memory zero idiom, loop replaced by call to __c_mzero8
     63, Barrier
         Loop not vectorized/parallelized: too deeply nested
     65, Parallel loop activated with static block schedule
         Loop not vectorized: may not be beneficial
         Unrolled inner loop 4 times
         Generated 3 prefetches in scalar loop
     69, Barrier
     70, Parallel loop activated with  runtime schedule schedule
         Loop not vectorized/parallelized: contains call
         Generated 1 prefetches in scalar loop
     71, Generated vector sse code for the loop
         Generated 2 prefetch instructions for the loop
     75, Barrier
         Parallel region terminated
     92, Loop not vectorized/parallelized: contains call
$ export OMP_NUM_THREADS=1
$ ./a.out
 OMP_Wtime resolution:    9.9999999999999995E-007
 Elapsed time =   0.3339718699455261       (sec)
 M =         1000 , N =         1000 , P =         1000
 MFLOPS =     5985.534052092637 (5.9GFLOPS)
 c(10,10) =     244.1990574713752
 
$ export OMP_NUM_THREADS=2
$ ./a.out
 OMP_Wtime resolution:    9.9999999999999995E-007
 Elapsed time =   0.1942952513694763       (sec)
 M =         1000 , N =         1000 , P =         1000
 MFLOPS =     10288.46554874702 (10.2GFLOPS)
 c(10,10) =     244.1990574713752
 
$ export OMP_NUM_THREADS=4
$ ./a.out
 OMP_Wtime resolution:    9.9999999999999995E-007
 Elapsed time =   0.1415630519390106       (sec)
 M =         1000 , N =         1000 , P =         1000
 MFLOPS =     14120.91624629021 (14.1GFLOPS)
 c(10,10) =     244.1990574713752

C++ での使用例 (LAPACKを使用)、Linux and OS X only

● ソースプログラム例

// LAPACK dgetrs test code
#include <iostream>
#include <vector>
using namespace std;

extern "C" void dgetrf_(int *M, int *N,  double *A, int *LDA, int *IPIV, int *INFO );
extern "C" void dgetrs_(char *TRANS, int *N, int *NRHS, double *A,
                      int *LDA, int *IPIV, double *B, int *LDB, int *INFO );

int main()
{
    char trans = 'N';
    int dim = 3;

    int info= 3;
    int LDA = dim;
    int LDB = dim;
    int n = 3;
    int nrhs = 1;

    // c++11 list initialization
    vector a {
    1, 2, 3,
    2, 3, 4,
    3, 4, 1
    };

    vector<double> b {
    -4,
    -1,
    -2
    };

    int ipiv[3];

    std::cout << "compute the LU factorization..." << endl << endl;
      dgetrf_(&n, &n, & *a.begin(), &LDA, ipiv, &info);
      dgetrs_(&trans, &n, &nrhs, & *a.begin(), &LDA, ipiv, &*b.begin(), &LDB, &info);

    std::cout << "Result is :";
    std::cout << "{" <<  " ";
        int i;
            for (i=0;i<n;i++)
            {
                std::cout << b[i] << " ";
            }
            std::cout << "}" << endl << endl;

    std::cout << "Info flag = " << info << std::endl;

    return(0);
}

● コンパイル・リンク&実行 (マルチスレッド用 -mp オプション付き)

$ pgc++ -w -fastsse main.cpp -std=c++11 -mp -lblas -llapack -pgf90libs 
$ a.out
compute the LU factorization...

Result is :{ 11 -9 1 }

Info flag = 0 

OpenMP 並列スレッド数指定の方法

 上記のコンパイルの方法において、-mp オプションを付加することにより、スレッド実行モードのモジュールが生成されます。BLAS/LAPACK ライブラリの中で OpenMP 並列対応となっているサブルーチンは、このライブラリ・モジュールの中で自動的に並列実行されます(全てのライブラリ・ルーチンが並列対応とはなっていませんのでご注意ください)。
並列実行での環境変数 OMP_NUM_THREADS を予め、以下のように設定してから実行してください(4スレッドの例)。

 $ export OMP_NUM_THREADS=4 
 $ 実行
------------------------------------------------------------------------------------
スタックサイズの不足等のエラーがある場合、スタックサイズの指定を行ってください