はじめに

Intel MKL (Math Kernel Library) は Intel が開発している数値演算ライブラリで、BLAS や LAPACK と同じインターフェースを持っており、高速化された BLAS や LAPACK としても知られています。

以前は Intel のコンパイラ製品に同梱の形でしか入手できなかったのですが、ライセンスが Intel Simplified Software License となって、取得や利用、再配布が容易になりました。そのため、各 Linux ディストリビューションでも Intel MKL が提供されるようになってきています。

この記事では、Ubuntu 20.04 の標準リポジトリから apt でインストールする、という前提で、Fortran と MKL を利用するための環境構築について紹介します。Ubuntu 20.04 は WSL2 でも提供されているので、Windows のユーザーにとっても参考にできるかと思います。

実際のところは、しばらく Fortran を使用してなかった私が久しぶりに環境構築しようとしたらいろいろと変わっていて戸惑ったので、そのとき調べた情報をまとめてみた、という話で、最近の Fortran の動向にはまったくついていけていないのですが、あまり手をかけないでサクッと Fortran + MKL の環境構築を行う、というところを目標にしています。

GFortran (GNU Fortran) と Intel MKL のインストール

apt install でのインストール

前提が apt からのフリーの環境構築なので、Fortran は GFortran を使います。GFortran は Free Software Foundation により開発されたフリーの Fortran で GCC (GNU Compiler Collection) の一部として配布されています。

GFortran, Intel MKL (以降、単に MKL と記載します) のいずれも、 Ubuntu 標準リポジトリから apt でインストールできるようになっていて、それぞれ gfortran, intel-mkl という名称で提供されています。

したがって、Ubuntu では単に apt install intel-mkl のようにしてインストールすることができます¹。

$ sudo apt update
$ sudo apt install gfortran
$ sudo apt install intel-mkl

intel-mkl のインストール (依存関係で多くのパッケージがインストールされます) 時には、下のように BLAS や LAPACK のデフォルトのライブラリとして MKL を使用するように設定するかどうかを訊かれます。配布するプログラムを作成するのでなければ Yes でよいでしょう。そうでなければライセンスを確認して選択してください。

pkg-config によるコンパイルオプションの取得

一般に、外部のライブラリを使用する際には、コンパイルオプションで使用するライブラリを指定する必要があります。規模の大きいライブラリの場合では、しばしば指定は複雑で理解が難しものになりがちです。

apt で MKL をインストールした場合には、GFortran に指定するオプションを pkg-config で取得することができます。

例えば、MKL の LP64 インターフェースで並列処理を行う場合のオプションを取得するには、ライブラリ名を mkl-static-lp64-iomp として、次のように pkg-condfig --libs mkl-static-lp64-iomp を実行すると必要なオプション文字列を得ることができます。

$ pkg-config --libs mkl-static-lp64-iomp
-Wl,--start-group /usr/lib/x86_64-linux-gnu/libmkl_intel_lp64.a /usr/lib/x86_64-linux-gnu/libmkl_intel_thread.a /usr/lib/x86_64-linux-gnu/libmkl_core.a -Wl,--end-group -liomp5 -lpthread -lm -ldl

pkg-config で指定可能なライブラリ名の一覧は pkg-config --list-all で確認できます。

$ pkg-config --list-all
udev                   udev - udev
tbb                    Threading Building Blocks - Intel's parallelism library for C++
mkl-static-lp64-seq    mkl - compatible pkg-config
mkl-static-ilp64-seq   mkl - compatible pkg-config
xkeyboard-config       XKeyboardConfig - X Keyboard configuration data
mkl-dynamic-lp64-iomp  mkl-dynamic-lp64-iomp - compatible pkg-config
blas64                 mkl-sdl-ilp64 - Debian Specific pkg-config for MKL SDL
blas                   mkl-sdl-lp64 - Debian Specific pkg-config for MKL SDL
blas-netlib            BLAS - FORTRAN reference implementation of BLAS Basic Linear Algebra Subprograms
systemd                systemd - systemd System and Service Manager
iso-codes              iso-codes - ISO country, language, script and currency codes and translations
mkl-static-ilp64-iomp  mkl-static-ilp64-iomp - compatible pkg-config
mkl-dynamic-ilp64-seq  mkl-dynamic-ilp64-seq - compatible pkg-config
libdmmp                libdmmp - Device mapper multipath management library
mkl-dynamic-ilp64-iomp mkl-dynamic-ilp64-iomp - compatible pkg-config
mkl-sdl-ilp64          mkl-sdl-ilp64 - Debian Specific pkg-config for MKL SDL
libxcrypt              libxcrypt - Extended crypt library for DES, MD5, Blowfish and others
mkl-static-lp64-iomp   mkl - compatible pkg-config
mkl-sdl-lp64           mkl-sdl-lp64 - Debian Specific pkg-config for MKL SDL
lapack64               mkl-sdl-ilp64 - Debian Specific pkg-config for MKL SDL
shared-mime-info       shared-mime-info - Freedesktop common MIME database
mkl-dynamic-lp64-seq   mkl-dynamic-lp64-seq - compatible pkg-config
lapack                 mkl-sdl-lp64 - Debian Specific pkg-config for MKL SDL
bash-completion        bash-completion - programmable completion for the bash shell
lapack-netlib          LAPACK - FORTRAN reference implementation of LAPACK Linear Algebra PACKage
libcrypt               libxcrypt - Extended crypt library for DES, MD5, Blowfish and others

システム全体に登録されている一覧なので、MKL と無関係なものも含まれていますが、mkl を接頭辞として static/dynamic, ilp64/lp64, iomp/seq の組み合わせを選択するようになっています。

static/dynamic は、ライブラリを静的リンク (static) とするか動的リンク (dynamic) とするかを選択します。

ilp64/lp64 は、配列インデックスのサイズに関わる指定です。大規模な配列を使用する場合は ILP64 (ilp64) を指定します。従来と互換性があるのは LP64 (lp64) です。詳細は公式ドキュメント「ILP64 インターフェイスと LP64 インターフェイスの使用」を参照してください。

iomp/seq は、並列計算を行う (iomp) か、それとも逐次計算を行う (seq) かを選択します。

つまり、先の例 mkl-static-lp64-iomp は、静的リンクでLP64インターフェースを使用して並列計算を行う、という指定になります。

サンプルコードによる MKL の動作確認

MKL の動作確認を行うために、LAPACK のサンプルプログラムを用意する必要があります。自分で書けばよいという話もありますが、LAPACK ののサンプルプログラムを白紙から書き起こすのは慣れないと苦痛です (少なくとも私は苦痛です)。

幸いにして Intel による LAPACK 用のサンプルコードが Intel® oneAPI Math Kernel Library LAPACK Examples で提供されているので、このコードを使って動作確認を行うことができます。

例えば、倍精度で連立一次方程式の解を求めるライブラリ DGESV の場合であれば、DGESV Example から、下記のような FORTRAN77 用のサンプルコードを取得することができます。

DGESV のサンプルコード dgesv_ex.f (クリックするとコードが表示されます)

*  Copyright (C) 2009-2015 Intel Corporation. All Rights Reserved.
*  The information and material ("Material") provided below is owned by Intel
*  Corporation or its suppliers or licensors, and title to such Material remains
*  with Intel Corporation or its suppliers or licensors. The Material contains
*  proprietary information of Intel or its suppliers and licensors. The Material
*  is protected by worldwide copyright laws and treaty provisions. No part of
*  the Material may be copied, reproduced, published, uploaded, posted,
*  transmitted, or distributed in any way without Intel's prior express written
*  permission. No license under any patent, copyright or other intellectual
*  property rights in the Material is granted to or conferred upon you, either
*  expressly, by implication, inducement, estoppel or otherwise. Any license
*  under such intellectual property rights must be express and approved by Intel
*  in writing.
*  =============================================================================
*
*  DGESV Example.
*  ==============
*
*  The program computes the solution to the system of linear
*  equations with a square matrix A and multiple
*  right-hand sides B, where A is the coefficient matrix:
*
*    6.80  -6.05  -0.45   8.32  -9.67
*   -2.11  -3.30   2.58   2.71  -5.14
*    5.66   5.36  -2.70   4.35  -7.26
*    5.97  -4.44   0.27  -7.17   6.08
*    8.23   1.08   9.04   2.14  -6.87
*
*  and B is the right-hand side matrix:
*
*    4.02  -1.56   9.81
*    6.19   4.00  -4.09
*   -8.22  -8.67  -4.57
*   -7.57   1.75  -8.61
*   -3.03   2.86   8.99
*
*  Description.
*  ============
*
*  The routine solves for X the system of linear equations A*X = B,
*  where A is an n-by-n matrix, the columns of matrix B are individual
*  right-hand sides, and the columns of X are the corresponding
*  solutions.
*
*  The LU decomposition with partial pivoting and row interchanges is
*  used to factor A as A = P*L*U, where P is a permutation matrix, L
*  is unit lower triangular, and U is upper triangular. The factored
*  form of A is then used to solve the system of equations A*X = B.
*
*  Example Program Results.
*  ========================
*
* DGESV Example Program Results
*
* Solution
*  -0.80  -0.39   0.96
*  -0.70  -0.55   0.22
*   0.59   0.84   1.90
*   1.32  -0.10   5.36
*   0.57   0.11   4.04
*
* Details of LU factorization
*   8.23   1.08   9.04   2.14  -6.87
*   0.83  -6.94  -7.92   6.55  -3.99
*   0.69  -0.67 -14.18   7.24  -5.19
*   0.73   0.75   0.02 -13.82  14.19
*  -0.26   0.44  -0.59  -0.34  -3.43
*
* Pivot indices
*      5      5      3      4      5
*  =============================================================================
*
*     .. Parameters ..
      INTEGER          N, NRHS
      PARAMETER        ( N = 5, NRHS = 3 )
      INTEGER          LDA, LDB
      PARAMETER        ( LDA = N, LDB = N )
*
*     .. Local Scalars ..
      INTEGER          INFO
*
*     .. Local Arrays ..
      INTEGER          IPIV( N )
      DOUBLE PRECISION A( LDA, N ), B( LDB, NRHS )
      DATA             A/
     $  6.80,-2.11, 5.66, 5.97, 8.23,
     $ -6.05,-3.30, 5.36,-4.44, 1.08,
     $ -0.45, 2.58,-2.70, 0.27, 9.04,
     $  8.32, 2.71, 4.35,-7.17, 2.14,
     $ -9.67,-5.14,-7.26, 6.08,-6.87
     $                  /
      DATA             B/
     $  4.02, 6.19,-8.22,-7.57,-3.03,
     $ -1.56, 4.00,-8.67, 1.75, 2.86,
     $  9.81,-4.09,-4.57,-8.61, 8.99
     $                  /
*
*     .. External Subroutines ..
      EXTERNAL         DGESV
      EXTERNAL         PRINT_MATRIX, PRINT_INT_VECTOR
*
*     .. Executable Statements ..
      WRITE(*,*)'DGESV Example Program Results'
*
*     Solve the equations A*X = B.
*
      CALL DGESV( N, NRHS, A, LDA, IPIV, B, LDB, INFO )
*
*     Check for the exact singularity.
*
      IF( INFO.GT.0 ) THEN
         WRITE(*,*)'The diagonal element of the triangular factor of A,'
         WRITE(*,*)'U(',INFO,',',INFO,') is zero, so that'
         WRITE(*,*)'A is singular; the solution could not be computed.'
         STOP
      END IF
*
*     Print solution.
*
      CALL PRINT_MATRIX( 'Solution', N, NRHS, B, LDB )
*
*     Print details of LU factorization.
*
      CALL PRINT_MATRIX( 'Details of LU factorization', N, N, A, LDA )
*
*     Print pivot indices.
*
      CALL PRINT_INT_VECTOR( 'Pivot indices', N, IPIV )
      STOP
      END
*
*     End of DGESV Example.
*
*  =============================================================================
*
*     Auxiliary routine: printing a matrix.
*
      SUBROUTINE PRINT_MATRIX( DESC, M, N, A, LDA )
      CHARACTER*(*)    DESC
      INTEGER          M, N, LDA
      DOUBLE PRECISION A( LDA, * )
*
      INTEGER          I, J
*
      WRITE(*,*)
      WRITE(*,*) DESC
      DO I = 1, M
         WRITE(*,9998) ( A( I, J ), J = 1, N )
      END DO
*
 9998 FORMAT( 11(:,1X,F6.2) )
      RETURN
      END
*
*     Auxiliary routine: printing a vector of integers.
*
      SUBROUTINE PRINT_INT_VECTOR( DESC, N, A )
      CHARACTER*(*)    DESC
      INTEGER          N
      INTEGER          A( N )
*
      INTEGER          I
*
      WRITE(*,*)
      WRITE(*,*) DESC
      WRITE(*,9999) ( A( I ), I = 1, N )
*
 9999 FORMAT( 11(:,1X,I6) )
      RETURN
      END

残念ながら Intel の当該ページは HTML で記述されているので、テキストをコピーして適当なエディタにペーストし、dgesv_ex.f というファイル名で保存するようにしてください。

サンプルプログラム dgesv_ex.f が用意できたら、下記のように gfortran を実行すると MKL とリンクされた実行ファイル dgesv_ex が作成されます。バッククオートで括った部分は pkg-config --libs mkl-dynamic-lp64-iomp の結果を gfortran のオプション文字列として渡すための指定です。ライブラリ名に dynamic の文字列が含まれていることからもわかるように、ここでは動的リンクとなるように指定しました。

$ gfortran -o dgesv_ex dgesv_ex.f `pkg-config --libs mkl-dynamic-lp64-iomp`

実行ファイル dgesv_ex を実行すると、下記のような出力が得られるはずです。ソースコード dgesv_ex.f のコメントの部分に想定される出力が記述されているので見比べてください。同じ結果となっているでしょうか。

$ ./dgesv_ex
 DGESV Example Program Results

 Solution
  -0.80  -0.39   0.96
  -0.70  -0.55   0.22
   0.59   0.84   1.90
   1.32  -0.10   5.36
   0.57   0.11   4.04

 Details of LU factorization
   8.23   1.08   9.04   2.14  -6.87
   0.83  -6.94  -7.92   6.55  -3.99
   0.69  -0.67 -14.18   7.24  -5.19
   0.73   0.75   0.02 -13.82  14.19
  -0.26   0.44  -0.59  -0.34  -3.43

 Pivot indices
      5      5      3      4      5

出力結果が正しいだけでは MKL が使われているかどうか確認できません。Netlib の方の LAPACK ライブラリが使われたかもしれないからです。動的リンクされた実行ファイルであれば、コマンド ldd を使用してどのライブラリが使用されているかを簡単に確認することができます。

ldd dgesv_ex を実行すると、下記のように MKL のライブラリ libmkl_* がリンクされていることを確認できます。

$ ldd dgesv_ex
        linux-vdso.so.1 (0x00007ffc9b5fa000)
        libmkl_intel_lp64.so => /lib/x86_64-linux-gnu/libmkl_intel_lp64.so (0x00007efd5916e000)
        libmkl_intel_thread.so => /lib/x86_64-linux-gnu/libmkl_intel_thread.so (0x00007efd56c02000)
        libmkl_core.so => /lib/x86_64-linux-gnu/libmkl_core.so (0x00007efd528e2000)
        libomp.so.5 => /lib/x86_64-linux-gnu/libomp.so.5 (0x00007efd527f1000)
        libgfortran.so.5 => /lib/x86_64-linux-gnu/libgfortran.so.5 (0x00007efd52529000)
        libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007efd52337000)
        libdl.so.2 => /lib/x86_64-linux-gnu/libdl.so.2 (0x00007efd5232f000)
        libpthread.so.0 => /lib/x86_64-linux-gnu/libpthread.so.0 (0x00007efd5230c000)
        /lib64/ld-linux-x86-64.so.2 (0x00007efd59cea000)
        libquadmath.so.0 => /lib/x86_64-linux-gnu/libquadmath.so.0 (0x00007efd522c2000)
        libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007efd52173000)
        libgcc_s.so.1 => /lib/x86_64-linux-gnu/libgcc_s.so.1 (0x00007efd52158000)

Intel MKL インストールのまとめ

以上で、apt install による GFortran と MKL のインストール、そしてその動作確認についての説明は終わりです。ライブラリの指定の方法などで若干わかりにくいところもありますが、導入そのものは以前に比べると非常に容易にできるようになったと思います。

記事の後半は Lapack95 についての話です。使えるようにはなっていないので、記録として残していますが読む必要はありません。

MKL 同梱の Lapack95 インターフェースについて

Lapack95 について

本記事での Lapack95 は、Intel MKL で提供されている LAPACK に対する Fortran95 用のインターフェースを指します。Netlib にも同様の目的で LAPACK95 というものが提供されていますが、両者に互換性はないようです。

結論から言うと、apt でインストールされる MKL の Lapack95 を GFortran で使用することはできませんでした。MKL をビルドしなおせば可能²なようですが、記事の趣旨から外れるので見送りました。

以下、失敗の記録なので折りたたんでいます。

モジュールファイル lapack95.mod の作成

Lapack95 を使用するためには lapack95.mod などのファイルが必要になります。apt で intel-mkl をインストールした場合、lapack95.mod は /usr/include/mkl/intel64/lp64, /usr/include/mkl/intel64/ilp64 などのディレクトリに設置されています。

しかし、Lapack95 を使用するコード lapacktest.f90 (use lapack95 を含むコード) を GFortran でコンパイルしようすると下記のようなエラーとなりました。lapack95.mod の読み込み自体に失敗しているようです。

$ gfortran -I/usr/include/mkl/intel64/lp64 lapacktest.f90 
f951: Fatal Error: Reading module ‘lapack95’ at line 1 column 2: Unexpected EOF
compilation terminated.

ここで使用した lapack95.mod は intel-mkl 同梱の lapack95.modです。つまり Intel Fortran 用として作成されたもので、GFortran とは互換性がなかったと考えられます。

Intel のドキュメント Fortran 95 Interface Conventions for LAPACK Routines によれば、MKL 同梱の lapack.f90 から lapack95.mod を作成すればよさそうです。lapack.f90 は /usr/include/mkl にありました。

$ gfortran -c /usr/include/mkl/lapack.f90
$ ls *.mod
f95_precision.mod  lapack95.mod

カレントディレクトリに lapack95.mod と f95_precision.mod が作成されました。他にも lapack.o が作成されましたが、こちらは使用しないはずです。

ライブラリ libmkl_lapack95_lp64.a のリンク

GFortran 用の lapack95.mod が用意できたので、再度 lapacktest.f90 のコンパイルを試みます。

lapack95.mod は読み込めたようですが、dsyev_f95 への参照が未定義であるとのエラーが出ました。

$ gfortran lapacktest.f90 `pkg-config --libs mkl-dynamic-lp64-iomp`
/usr/bin/ld: /tmp/ccaCroPn.o: in function `test_':
lapacktest.f90:(.text+0xbc8): undefined reference to `dsyev_f95_'
/usr/bin/ld: lapacktest.f90:(.text+0xc64): undefined reference to `dsyev_f95_'
/usr/bin/ld: lapacktest.f90:(.text+0x1b73): undefined reference to `zheev_f95_'
/usr/bin/ld: lapacktest.f90:(.text+0x1c0f): undefined reference to `zheev_f95_'
collect2: error: ld returned 1 exit status

dsyev_f95 は LAPACK のルーチン DSYEV へのラッパーであることが lapack.f90 のコードから推測できます。どうやら MKL のコアライブラリには dsyev_f95 が含まれていないようです。

コマンド nm を使用すると動的ライブラリ内で定義されているシンボルを確認することができます。dsyev_f95 は libmkl_lapack95_lp64.a, libmkl_lapack95_ilp64.a で定義されていました。

オプション -lmkl_lapack95_lp95 を追加してコンパイルを実行すると、dsyev_f95 は解決したようですが、ライブラリが PIE オブジェクトではないというエラーが発生しました。

$ gfortran lapacktest.f90 -lmkl_lapack95_lp64 `pkg-config --libs mkl-dynamic-lp64-iomp`
/usr/bin/ld: /usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu/libmkl_lapack95_lp64.a(dsyev.o): relocation R_X86_64_32 against `__STRLITPACK_2' can not be used when making a PIE object; recompile with -fPIE
/usr/bin/ld: /usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu/libmkl_lapack95_lp64.a(zheev.o): relocation R_X86_64_32 against `__STRLITPACK_2' can not be used when making a PIE object; recompile with -fPIE
collect2: error: ld returned 1 exit status

PIE とは Position Independent Executable (位置独立実行形式) を表します³。エラーはライブラリに含まれるオブジェクトの形式の違いによるもので、手っ取り早く回避するにはオプション -no-pic を指定します。

未定義関数 for_check_mult_overflow64, for_allocate

オプション -no-pie を追加してコンパイルを実行したところ、今度は for_check_mult_overflow64, for_allocate が未定義であるとのエラーが発生しました。

$ gfortran lapacktest.f90 -lmkl_lapack95_lp64 `pkg-config --libs mkl-dynamic-lp64-iomp`
 -no-pie
/usr/bin/ld: /usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu/libmkl_lapack95_lp64.a(dsyev.o): in function `dsyev_f95_':
dsyev.f90:(.text+0x83f): undefined reference to `for_check_mult_overflow64'
/usr/bin/ld: dsyev.f90:(.text+0x875): undefined reference to `for_allocate'
/usr/bin/ld: dsyev.f90:(.text+0x8cd): undefined reference to `for_check_mult_overflow64'
/usr/bin/ld: dsyev.f90:(.text+0x1318): undefined reference to `for_dealloc_allocatable'
/usr/bin/ld: /usr/lib/gcc/x86_64-linux-gnu/9/../../../x86_64-linux-gnu/libmkl_lapack95_lp64.a(zheev.o): in function `zheev_f95_':
zheev.f90:(.text+0xad): undefined reference to `for_check_mult_overflow64'
/usr/bin/ld: zheev.f90:(.text+0xe3): undefined reference to `for_allocate'
/usr/bin/ld: zheev.f90:(.text+0x15f): undefined reference to `for_check_mult_overflow64'
/usr/bin/ld: zheev.f90:(.text+0x972): undefined reference to `for_dealloc_allocatable'
/usr/bin/ld: zheev.f90:(.text+0xa92): undefined reference to `for_dealloc_allocatable'
/usr/bin/ld: zheev.f90:(.text+0xada): undefined reference to `for_check_mult_overflow64'
/usr/bin/ld: zheev.f90:(.text+0xb13): undefined reference to `for_allocate'
/usr/bin/ld: zheev.f90:(.text+0xb7a): undefined reference to `for_check_mult_overflow64'
collect2: error: ld returned 1 exit status

先程と同様に for_check_mult_overflow64, for_allocate が定義されているライブラリをリンクすれば動くのではないかと考えられますが、MKL でインストールされているライブラリの中には見つけることはできませんでした。もしかすると Intel Fortran の中で定義されているのかもしれません。

Lapack95 インターフェース使用についてのまとめ

apt でインストールした MKL で Lapack95 インターフェースを利用しようと試みましたが、未定義シンボルが解決できず、実行ファイルを作成するところまで進めることはできませんでした。GFortran で MKL のビルドを行えば Lapack95 を利用できる可能性がありますが、手軽さは大きく損なわれてしまうので、別の機会に検討したいと思います。

Qiita: 「Intel MKLの各Linuxへのインストール方法」では他の Linux ディストリビューションの場合や、Intel 公式の apt リポジトリを使う方法が紹介されています。 ↩
Qiita: 「MKLのLapack95インターフェースを使ってFortranで行列を対角化しよう」では Lapack95 のビルドを行っているようです。残念なことに参照している公式ドキュメントのリンクが切れてしまっていて一部が確認できません。 ↩
Qiita: 「makeで "can not be used when making a PIE object" が出たら」にわかりやすい説明があります。 ↩