VSCodeで最小限のLaTeX環境を整備する

Visual Studio Code（以下、VSCode）1.5でLaTeX環境を整備したのですが、参考にしたページの設定が心理的な参入障壁になりそうなほど頑張っていたので、備忘録をかねてWindows 10でplatex/pbibtexを使うだらっとした最小限の設定を紹介します。

1. LaTeXのインストールと設定

4GBあるTeX Live 2020のISOイメージをダウンロードしてきてインストールします。Windows 10であれば、ダウンロードしたファイルを右クリックしてマウントすると、ドライブが割り当てられるので、ドライブを開いてinstall-tl-windows.batを実行します。

イマドキはlatexmkコマンドを使ってTeXのコンパイルをするので、その設定になる%USERPROFILE%\.latexmkrcを作ります。notepadに、

$latex = 'platex -synctex=1 %O %S';
$bibtex = 'pbibtex %O %B';
$dvipdf = 'dvipdfmx %O -o %D %S';
$makeindex = 'mendex %O -o %D %S';
$max_repeat = 10;

と書いて、ファイル種類を*.*、ファイル名を%USERPROFILE%\.latexmkrcで、名前をつけて保存してください。これでlatexmk -pdfdvi example.texで、example.pdfが作れるようになります。
分野ごとの参考文献の表示スタイル（e.g. jecon.bst）などもこの段階で追加しておきましょう*1。

2. VSCodeのインストールと設定

VSCodeの公式サイトからダウンロードしてきて、インストールします。続けて、LaTeX WorkshopとJapanese Language Pack for Visual Studio Codeの2つの拡張機能を、ショートカットCTRL + SHIFT + Xを押すか、以下のアイコンをクリックしてVS Code Marketplaceのペーンを呼び出し、検索してインストールします *2。

これで.texファイルと.bibファイルの編集でシンタックスハイライトと言うか色分け表示してくれるようになり、コード補間が効くようになります。

次に、VSCodeで.texファイルを編集中にウィンドウ内右上の三角形のRunボタンを押したときの挙動を設定します。notepadかVSCodeで%APPDATA%\Code\User\settings.jsonを開いて、

    // for LaTeX
    "latex-workshop.latex.recipes": [
        {
            "name": "latexmk",
            "tools": [
                "latexmk"
            ]
        },
    ],
    "latex-workshop.latex.tools": [
        {
            "name": "latexmk",
            "command": "latexmk",
            "args": [
                "-pdfdvi",
                "%DOC%",
            ],
        },
    ],

をファイル内先頭の{と最後の}の間に挟みます。なお、このファイルはフォントサイズやタブ幅などVSCodeのエディターの諸々の設定を書いておける場所で作りこむことが可能です*3が、上のようなコマンド変数以外は、VSCodeでCTRL + ,を押してダイアログ的な画面を使った方がよいかも知れません。

*1: jecon.bstの場合は、jecon.bstを落としてきて、C:\texlive\2020\texmf-dist\pbibtex\bstに放り込むだけになると思います。

*2: Managing Extensions in Visual Studio Code

*3: vscodeがいい感じになる設定を継ぎ足していく（俺流vscode秘伝のタレ） - Qiita

OpenBLASをリンクしたWindows版R 4.1 PR

ループでスカラ演算を繰り返すよりもベクトル演算をした方が速いRですが、標準構成ではNetlib BLASと言う線形代数ライブラリBLASのリファレンス実装を使っていて、行列演算、ベクトル演算も速いとは言えないものとなっています。チューニングされたOpenBLASに差し替えることで、高速化をしてみましょう。

Ubuntu Linuxでは、

sudo apt install libopenblas0-pthread

で、OpenBLASがインストールされ、aptでインストールしたRはOpenBLASを使うようになります。簡単ですね。しかし、Windowsだと時間のかかる作業が必要です。
基本的にはチートペーパーを見ながら、Rの開発版をコンパイルできる環境を揃えて、ビルドに使う設定ファイルを書き換えれば終わりなのですが、微妙に書いてある通りではコンパイルできませんでした。

まず、Step 4: Create blas patchのところですが、現行のMakefile.winがホワイトスペースに空白ではなくタブを使っているせいか、示されているblas.diffではパッチがあたりません。blas.diffをつくりなおす必要がありました。作り直したblas.diffをr-base-master.zipを展開してできるフォルダーC:\r-base-masterの中に入れました。

次に、Step 5: Adjust existing filesのところですが、示されているPKGBUILDの末尾の</jeroen@berkeley.edu>は要らないので消しました。見比べるとmakedependsに"${MINGW_PACKAGE_PREFIX}-openblas"、sourceにblas.diff、sha256sumsに'SKIP'、# Add your patches hereにpatch -Np1 -i "${srcdir}/blas.diff"をそれぞれ追加しているだけのようなので、チートペーパーのPKGBUILDをコピペする必要は無いかもです。（これはチートペーパーに記述の問題はありませんが）MkRules.local.inのQPDFへのパスは、

QPDF = C:/qpdf-10.0.1

と、インストール先のQPDFのホームにしました。
export PATH="$PATH:/c/progra~1/MiKTeX 2.9/miktex/bin/x64"の行は、

export PATH="$PATH:/c/progra~1/MikTeX/miktex/bin/x64"

と、（自明ですが）MikTeXのインストール先にあわせました。

設定は以上で終わりですが、別のチートペーパーを参考にビルドする前にMikTeXでフォントを作成しておかないとFont ts1-zi4r at 540 not foundとエラーが出て止まります*1。

Step 6: Build Rはマイナーなところですが、rtools40のmsys.exeを使えと書いてありますが、手元の環境ではC:\rtools40\msys2.exeでした。
なお、コンパイル途中でQPDFがインストール確認をしてくるので、少なくとも1回は了解ボタンを押すまでは、放置しておいてもコンパイルは終了しません。また、マイナーな計算結果差のようですが、微妙に1つのテストでエラーが出ました。これは4.1 PRのせいかも知れません*2。

コンパイルが終わってR-devel-win.exeができたらすぐにインストールしてみたのですが、Windows版ではsessionInfo()をしても使っているBLASの種類を教えてくれないです。以下のように行列演算をさせ、4コアCPUで35倍以上の経過時間差を出して、差し替えた感を実感しました。

set.seed(1013)
n <- 5000
M1 <- matrix(runif(n^2, min=1, max=10), n, n)
M2 <- matrix(runif(n^2, min=1, max=10), n, n)
system.time({ M3 <- M1 %*% M2 })

ところで先日書いてみたマクロ経済学の数値解析は、ほとんど速くなりませんでした。速くなりませんでした。

*1:チートペーパーのコマンドをうってinitexmf: security risk: running with elevated privilegesとエラーが出る人は、付属ユーティリティMiKTeX Consoleで操作した方がよいかもです。起動してadministration modeを選択し、メニューのTasksのRefresh file name database、Refresh font map filesを順番に実行します。

*2:export rsource_url='https://cran.r-project.org/src/base-prerelease/R-latest.tar.gz'; ./full-build.shとコンパイルを実行すると、開発版ではなくリリース候補版のソースコードを指定できます。

Rの拡張でOpenMPを使ってみる

プロセッサのマルチコア化が進んだ現代なので、時間のかかる計算では並列処理をするコードを書くのが望ましいです。スレッドやセマフォの制御を直接プログラミングすると骨が折れるのですが、シングルコア向け逐次処理コードを並列処理に変換してくれるOpenMPという仕組みがあります。Cで書いたRの拡張で使ってみましょう。

1. Makevarsの編集

gccでOpenMPを使うためには、引数に-fopenmpフラグを加える必要があります。R CMDでこの引数をつけるためには、LinuxやMacOS Xでは$HOME/.R/Makevarsに、MS-Windowsでは%HOME%*1/.R/Makevars.winと言うテキストファイルに、

CFLAGS= -fopenmp
CXXFLAGS= -fopenmp
# PKG_CXXFLAGSか、PKG_CXXFLAGSに-std=c++11と書いてCXX11FLAGSに-fopenmpをつける（C++14/17でも同様）ほうがよいときもあります。

と言うようにかいておきます。なお、-O3オプションなども同様につきます。

2. Cのソースコードの編集

ループされていて一定以上の処理時間がかかるコードの例として、2つの一行を除けば同じ（無意味な）関数をOMPexample.cとして保存します。

#include<R.h>
#include<Rinternals.h>

SEXP example(SEXP m){
  int i, j, n;
  double s1 = 0, s2 = 0;
  int nor, noc;
  SEXP rv, dim;

  /* 第一引数が行列か確認 */
  if(!isMatrix(m))
    error("A matrix is required for the first argument.");

  /* 行列の行数、列数を得る */
  dim = getAttrib(m, R_DimSymbol);
  nor = INTEGER(dim)[0]; /* 行数 */
  noc = INTEGER(dim)[1]; /* 列数 */

  for(j=0;j<noc;j++){
    for(i=0;i<nor;i++){
      s1 += REAL(m)[i+j*nor]*(i+1)/(j+1); /* 行列と言っても一次元配列 */
      s2 -= REAL(m)[i+j*nor]*(i+1)/(j+1);
    }
  }

  rv = PROTECT(allocVector(REALSXP, 2));
  REAL(rv)[0] = s1;
  REAL(rv)[1] = s2;
  UNPROTECT(1);

  return rv;
}

SEXP OMPexample(SEXP m){
  int i, j, n;
  double s1 = 0, s2 = 0;
  int nor, noc;
  SEXP rv, dim;

  /* 第一引数が行列か確認 */
  if(!isMatrix(m))
    error("A matrix is required for the first argument.");

  /* 行列の行数、列数を得る */
  dim = getAttrib(m, R_DimSymbol);
  nor = INTEGER(dim)[0]; /* 行数 */
  noc = INTEGER(dim)[1]; /* 列数 */

  #pragma omp parallel for private(i, j) reduction(+:s1) reduction(-:s2)
  for(j=0;j<noc;j++){
    for(i=0;i<nor;i++){
      s1 += REAL(m)[i+j*nor]*(i+1)/(j+1); /* 行列と言っても一次元配列 */
      s2 -= REAL(m)[i+j*nor]*(i+1)/(j+1);
    }
  }

  rv = PROTECT(allocVector(REALSXP, 2));
  REAL(rv)[0] = s1;
  REAL(rv)[1] = s2;
  UNPROTECT(1);

  return rv;
}

大学の情報基盤センターなるところが書いたネット上にあるOpenMPの解説を読めばすぐ分かりますが、プリプロセッサ命令

  #pragma omp parallel for private(i, j) reduction(+:s1) reduction(-:s2)

がOpenMPの記述になります。private(i, j)は並列処理で変数iとjは、それぞれ異なる値をシェアしない変数として扱うことを意味し、reduction(+:s1)は、それぞれの並列処理で加算されていくそれぞれの変数s1は、並列処理の終了時に一つに合算されることを意味します。reduction(-:s2)は、加算ではなく減算処理をあわせます。
なお、OpenMPが有効でない場合は、このプリプロセッサは無視されて逐次処理されるコードが出力されます。

3. Cのソースコードをコンパイルする

OMPexample.cを用意したら、以下のようにすると、OMPexample.so（もしくはOMPexample.dll）と言うRの拡張が出来上がります。

R CMD SHLIB OMPexample.c

Rtoolsが入っていないと動かないので注意してください。

4. Rから拡張を呼び出して使う

OpenMPを使ったOMPexample関数と、使っていないexample関数の速度比較をRでしてみましょう。

dll <- dyn.load(paste("OMPexample", .Platform$dynlib.ext, sep = ""))
example <- function(m) .Call("example", m)
OMPexample <- function(m) .Call("OMPexample", m)

set.seed(1002)
n <- 5000 # 一瞬で処理が終わる人は増やすと時間をかけられます
m <- matrix(rnorm(n*n), n, n)

print(system.time({
  cat("No OpenMP: ", example(m), "\n")
}))

print(system.time({
  cat("OpenMP: ", OMPexample(m), "\n")
}))

搭載されているプロセッサ数に応じて速くなります。並列化と言っても、手軽に試せますね。

A. Rcpp/Eigen

Rcppの行列はスレッドセーフではないのでOpenMPの並行処理の中で使えませんが、Eigenは最初にEigen::initParallel();をして行列をスレッドセーフにしておくだけで、OpenMPの並列処理の中で行列を使うことができます。つまり、Rcppの行列からEigenの行列を作成して処理を行うことで、RcppからOpenMPを使うことができます。RcppPararellを使うべきと思うかも知れませんが、RcppPararellはWorkerクラスを定義しないといけないので記述が煩雑な一方、Eigen/OpenMPと比較して速度的なアドバンテージはほぼ無いようです。

*1:OSで用意される環境変数ではないので、Rで Sys.getenv("HOME") をして確認してください。

疎な巨大対称行列の一般化固有値問題のC++/Eigenでの効率的な解き方

「水素原子に対するSchrödinger方程式を有限要素法で数値的に解いてみる（C++のソースコード付き）」で紹介されているC++/Eigenのコードを、（ページ著者が）Juliaに移植してみたら12倍くらい速くなったと言うツイートが流れてきて、気になってC++の方のソースコードを見てみたのですが、

• 巨大な疎行列であること
• 必要な固有値が1つであること*1

を活かしておらず、この2点を利用すればもっと時間を短縮できるものでした*2。実際、そう大きくない変更で、手元の環境では25分の1の時間で計算が終わるようになりました。

追記(2021/02/15 01:42):ブログ主が、有限要素法による水素原子のエネルギー固有値が、小さい方から見てどこまで正確かを検証する部分を加筆しました（図1, 図2）。これによって行列の固有値がすべて必要になるので、以下の方法は適切な計算方法にならないです。

なお、すべての固有値を得る場合にJulia比較でC++/Eigenが遅いのは、どうもC++/Eigenが自動的に並列化してマルチコアを活かしてこないところにあるので、その辺を何らかの方法*3で解消すると高速化できます。

EigenはC++用の線形代数ライブラリなのですが、名前がEigenなのに疎行列（Sparse Matrix）用の固有値問題のソルバーが無かったりします。上述のコードでは、GeneralizedSelfAdjointEigenSolverと言う計算する行列式に、疎行列、対称、正定値と言った条件を使わない上に、行列の固有値をすべて用いようとするソルバーが用いられています。汎用ソルバーとしても遅いというような指摘もあるようですが、入力データの特性を活かすのは高速化のイロハです。

この手の数値演算アルゴリズムを実装しようとすると数学の理解やデバッグの手間隙に大変な労力を費やすことになりますが、幸い、手軽に使えるテンプレート・ライブラリでSpectraと言うのがあります。具体的には、一般化固有値問題の計算対象の2つの行列をEigen::SparseMatrixで持ち、SymGEigsSolverと言うソルバーと解きます。なお、使ってみると癖があるなと思う部分があって、特にコンストラクタの第4引数を小さく取りすぎると固有値が収束しなかったとエラーを出して来るところには注意が必要です。

差分は以下になります。コンパイル時にダウンロードして展開したSpectraへインクルード・パスを通さないとコンパイル・エラーになるので気をつけてください。

--- a/src/hydrogen_fem/hydrogen_fem.cpp
+++ b/src/hydrogen_fem/hydrogen_fem.cpp
@@ -12,19 +12,24 @@
#include <memory> // for std::unique_ptr
#include <boost/assert.hpp> // for BOOST_ASSERT
#include <Eigen/Eigenvalues> // for Eigen::GeneralizedSelfAdjointEigenSolver

+#include <Eigen/Sparse>
+#include <Spectra/SymGEigsSolver.h>
+#include <Spectra/MatOp/SparseGenMatProd.h>
+#include <Spectra/MatOp/SparseCholesky.h>
+#include <iostream>
namespace hydrogen_fem {
// #region コンストラクタ

Hydrogen_FEM::Hydrogen_FEM()
- : hg_(Eigen::MatrixXd::Zero(NODE_TOTAL, NODE_TOTAL)),
+ : hg_(Eigen::SparseMatrix<double>(NODE_TOTAL, NODE_TOTAL)),
length_(ELE_TOTAL),
mat_A_ele_(boost::extents[ELE_TOTAL][2][2]),
mat_B_ele_(boost::extents[ELE_TOTAL][2][2]),
nod_num_seg_(boost::extents[ELE_TOTAL][2]),
node_r_ele_(boost::extents[ELE_TOTAL][2]),
node_r_glo_(NODE_TOTAL),
- ug_(Eigen::MatrixXd::Zero(NODE_TOTAL, NODE_TOTAL))
+ ug_(Eigen::SparseMatrix<double>(NODE_TOTAL, NODE_TOTAL))
{
}

@@ -42,20 +47,33 @@ namespace hydrogen_fem {

// 全体行列を生成
make_global_matrix();

// 一般化固有値問題を解く
- Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> es(hg_, ug_);

+  Spectra::SparseSymMatProd<double> Aop(hg_);
+  Spectra::SparseCholesky<double> Bop(ug_);
+
+  // 第3引数の1は要求する固有値の数，第4引数は収束速度（とメモリー利用量と行列操作の複雑さ）を定める値
+  Spectra::SymGEigsSolver<double, Spectra::SMALLEST_ALGE, Spectra::SparseSymMatProd<double>, Spectra::SparseCholesky<double>, Spectra::GEIGS_CHOLESKY> geigs(&Aop, &Bop, 1, NODE_TOTAL<10 ? NODE_TOTAL : 10 + (NODE_TOTAL - 10)/100);
+
+  geigs.init();
+  int nconv = geigs.compute();
+  if(Spectra::SUCCESSFUL != geigs.info()){
+    std::cerr << "geigs.info(): " << geigs.info() << std::endl;
+    exit(-1);
+  }

// エネルギー固有値Eを取得
- auto const e = es.eigenvalues()[0];
+ Eigen::VectorXcd evalues = geigs.eigenvalues();
+ auto const e = evalues[0].real(); // 計算結果は複素数なことに注意．なお、複数の固有値を求めた場合、なぜか小さい順に固有値が並んでいるわけではない.

// 固有ベクトル（波動関数）を取得
- phi_ = es.eigenvectors().col(0);
+ phi_ = geigs.eigenvectors(1).col(0).real(); // なぜか添字が1からなのと、計算結果は複素数なのに注意

// 固有ベクトル（波動関数）を規格化
normalize();

return e;
}

void Hydrogen_FEM::save_result() const
@@ -190,8 +208,8 @@ namespace hydrogen_fem {
for (auto e = 0; e < ELE_TOTAL; e++) {
for (auto i = 0; i < 2; i++) {
for (auto j = 0; j < 2; j++) {
- hg_(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_A_ele_[e][i][j];
- ug_(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_B_ele_[e][i][j];
+ hg_.coeffRef(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_A_ele_[e][i][j];
+ ug_.coeffRef(nod_num_seg_[e][i], nod_num_seg_[e][j]) += mat_B_ele_[e][i][j];
}
}
}

--- a/src/hydrogen_fem/hydrogen_fem.h
+++ b/src/hydrogen_fem/hydrogen_fem.h
@@ -15,6 +15,9 @@
#include <Eigen/Core> // for Eigen::MatrixXd
#include <boost/multi_array.hpp> // for boost::multi_array

+#include <Eigen/Core>
+#include <Eigen/SparseCore>
+
namespace hydrogen_fem {
//! A class.
/*!
@@ -144,7 +147,7 @@ namespace hydrogen_fem {
/*!
左辺の全体行列
*/
- Eigen::MatrixXd hg_;
+ Eigen::SparseMatrix<double> hg_;

//! A private member variable.
/*!
@@ -192,7 +195,7 @@ namespace hydrogen_fem {
/*!
右辺の全体行列
*/
- Eigen::MatrixXd ug_;
+ Eigen::SparseMatrix<double> ug_;

// #region 禁止されたコンストラクタ・メンバ関数

*1:すべての固有値を取りたい場合もあるわけですが、大学のレジュメなどを拝見する限りは、最小値か最大値のどちらか一つを計算すれば間に合うケースが多いようです。

*2:最初からすべて分かったように書いていますが、ページ著者と試行錯誤して、この結論にたどり着いています。

*3:OpenBlASを使う方法が紹介されていました。

FASTGLを使ってC++で積分

数値解析における数値積分法の代表としてガウス求積と言う技法があり、それに必要なガウス点と重みの計算でJuliaがC++より桁違いに速いと言う話が流れていたのですが、Juliaで使われているFastGaussQuadrature.jlが使っているアルゴリズムのC++版を使ったところ、案の定、だいたい同じの速度になりました。

Fast computation of Gauss-Legendre Quadrature Points and Weights（FASTGL）のアルゴリズム開発者のページからfastgl.cppとfastgl.hppをダウンロードしてきて、積分に使うだけです。サブディレクトリfastglにこの2つのファイルを配置したところで、以下のソースコードのファイルを作成してコンパイルしたら簡単に実行できます。-O3 -lmオプションはいると思いますが。

#include <iostream>
#include <chrono>
#include <iomanip>
#include <functional>
#include "fastgl/fastgl.hpp"

/*
  FASTGLを使って積分をする関数
  args: ラムダ関数f(x), xの最小値, xの最大値, 節点の数
*/
double integral(std::function<double(double)> f, double s, double e, int NumberOfKnots){
  double Res = 0.0;
  double xm = 0.5 * (s + e);
  double xr = 0.5 * (e - s);

  for ( int k = 1 ; k <= NumberOfKnots ; ++k ){
    fastgl::QuadPair p = fastgl::GLPair ( NumberOfKnots, k );
    Res += p.weight * f( xm + xr * p.x());
  }

  return Res*xr;
}

int main(){
  std::chrono::system_clock::time_point stime, etime;
  stime = std::chrono::system_clock::now();

  // boost::math::quadratureあたりと異なり、節点が100万個でも大した時間にならない
  double Res = integral([](const double& t) { return 1.0 / 2 / std::sqrt(t); }, 1, 4, 1000000);

  etime = std::chrono::system_clock::now();

  std::cout << std::setprecision(16);
  std::cout << "Results: " << Res << std::endl;
  std::cout << "Elapsted Time: " << std::chrono::duration_cast<std::chrono::milliseconds>(etime - stime).count() << " ms" << std::endl;

  return 0;
}

モンテカルロ法で円周率を計算するベンチマーク

モンテカルロ法で円周率を計算するベンチマークで、Juliaと比較してgccが悪い言われているのですが、本当なのか検証してみましょう。

CとJuliaでこの分岐が単純なマイクロベンチマークを実行すると、乱数生成の速度が結果を大きく左右しえます。後方互換性の問題で、Cの標準の乱数生成関数は線形合同法になっています。Juliaの方は標準がdouble precision SIMD oriented Fast Mersenne Twister pseudorandom number generator (dSFMT)で、明示的にMersenneTwisterを使った場合もdSFMTを用いているようです。Cの方もdSFMTを用いてベンチマークを取ることにします。

Cのソースコード

dSFMTがよくチューニングされているので、私には工夫の余地は無く、オーソドックスに書いています。なお、生成した実行ファイルをreadelfで調べたところ、dSFMTがヘッダーファイルでinline static宣言している関数はインライン展開されていました。10回計測して、一番短い経過時間を採用しています。

#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <limits.h>
#include "dSFMT.h"
#include <float.h>
#include <math.h>

/*
dSFMT-src-2.2.3.tgz を展開したディレクトリで、ソースコードのファイル名をpi_mc_dsfmt.cとして、

clang -o pi_mc pi_mc_dsfmt.c dSFMT-src-2.2.3/dSFMT.c -IdSFMT-src-2.2.3/ -DDSFMT_MEXP=521 -msse2 -DHAVE_SSE2 -O3 -lm

とコンパイルしました。gccとclangはオプションが同じなので、gccとの比較は容易に行えます。
*/

double pi_mc(long n_iter) {
  double  x, y;
  dsfmt_t rng;

  dsfmt_init_gen_rand(&rng, time(NULL));

  int n_in = 0;
  for (long i = 1; i <= n_iter; i++) {

    x = dsfmt_genrand_open_open(&rng);
    y = dsfmt_genrand_open_open(&rng);

    if (x * x + y * y < 1.0)
      n_in++;
  }

  return (double)4.0*n_in/n_iter;
}

int main(){
  struct timespec t1, t2;
  double  dt, sum_dt, min_dt;
  double  r = 0;
  size_t  n;
  int  j, nos = 10;

  printf("log10(n)\telapsed time (seconds)\n");

  for(n=1; n<=1000000*1000; n*=10){
    sum_dt = 0;
    min_dt = DBL_MAX;
    for(j=0; j<nos; j++){
      clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t1); // clock()を使って計測しても、ほぼ同じ数字になります

      r = pi_mc(n);

      clock_gettime(CLOCK_THREAD_CPUTIME_ID, &t2);

//      fprintf(stderr, "estimated pi: %f\n", r);

      dt = (double)(t2.tv_sec - t1.tv_sec) + (double)(t2.tv_nsec - t1.tv_nsec)/1000000000;
      
      sum_dt += dt;
      if(min_dt > dt)
        min_dt = dt;
    }

    printf("%.0f\t%.09f\n", log10(n), min_dt);
  }

  return 0;  
}

Juliaのスクリプト

Cのとだいたい同じですが、グローバル変数へのアクセスが無いように一番外側ループ部分も関数内にしまいました。

using Random
# using Random: default_rng

function pi_mc(n_iter)
n_in = 0
# rng = default_rng()
rng = MersenneTwister(time_ns())
for i in 1:n_iter
x, y = rand(rng), rand(rng)
n_in += ifelse(x^2 + y^2 ≤ 1, 1, 0)
end
4n_in/n_iter
end

function loop(s, e)
  n = s
  while(n <= e)
    sum_dt = 0;
    min_dt = 3600;
    for i in 1:10
      t1 = time_ns()
      r = pi_mc(n)
      t2 = time_ns()
      dt = (t2 - t1)/1000000000
      sum_dt = sum_dt + dt
      if(min_dt > dt)
        min_dt = dt
      end
#      println("estimated pi: ", r)
    end
    println(log10(n), "\t", min_dt)
    n = n*10
  end
end

println("log10(n)", "\t", "elapsed time (seconds)")

loop(1, 10^9)

SIMD拡張命令セットでgcc用のsum関数を書いてみる

とあるJuliaを使った教材に、JuliaやNumPyのsumはSIMD拡張命令セットを使うから速いけれども、C言語でsumを書くとそうではないから遅いような記述があって、CからSIMDを叩けば良いだろうと思ったので、私にしては珍しくC99のコードで、単浮動小数点のsumを書いてみました。gccでしかコンパイルできません。

謳い文句によるとgccは自動的にSIMD拡張命令セットを使うことになっていますが、実際のところはsumのような単純なものでも使いません。こんなにSIMD向きの単純計算でも、128bitsを32bits×4で使っても2倍ちょっとしか速くならないわけで、倍精度だとAVX-256以降でないと意味がないからでしょうか。

なお、C++のテンプレート・ライブラリEigenは行列演算などにSIMD拡張命令セットを使ってくれるので、こういう苦しい記述をしないで浮動小数点演算で速度を出したい人はC++/Eigenを使うべきです。CPUのサポートしている命令セットごとにコードを書き換えないといけないって、手間隙な上に混乱してきますからね。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <emmintrin.h>
#include <time.h>
#include <cpuid.h>

float  x87_sum(float *X, int n){
float s = 0.0;
for (size_t i = 0; i < n; ++i) {
s += X[i];
}
return s;
}

float  simd_sum(float *X, int n){
  static int ws = sizeof(__m128)/sizeof(float);
  __m128  r128, s128 = {0};
float s = 0.0, sws[ws];
  int m = m % ws;

// add by ws-float from [0 n-m]
for (size_t i = 0; i < n - m; i+=ws) {
r128 = _mm_load_ps(&X[i]);
    s128 = _mm_add_ps(r128, s128);
}
// store result of the above code.
  _mm_store_ps(sws, s128);

// sum sws up, and record it into s.
for (size_t i = 0; i < ws; i++) {
    s += sws[i];
  }

// add the rest, [n-m+1 n-1], by a float.
for (size_t i = 1; i <= m; i++) {
    s += X[n - i];
  }

return s;
}

int main(){
  float *X;
  int  n = 10000;
  int i;
  long et;
  int  ws;
  double r1, r2;
  union cpuid_sig {
    unsigned int  i;
    unsigned char  c[4];
  } sig;
  unsigned int highest, eax, ebx, ecx, edx;

  highest = __get_cpuid_max(0x8000000, &sig.i);
  if(0 < __get_cpuid(highest, &eax, &ebx, &ecx, &edx) && 0 != bit_SSE2 & edx){
    printf("SSE2 is supported.\n");
  } else {
    fprintf(stderr, "SSE2 is NOT supported!\n");
    exit(-1);
  }

  ws = sizeof(__m128)/sizeof(float);
  printf("sizeof(__m128)/sizeof(float) = %ld/%ld = %d (%ld)\n", sizeof(__m128), sizeof(float), ws, sizeof(__m128) - sizeof(float)*ws);

// setup X
  X = calloc(n, sizeof(float));
  for(i=0; i<n; i++){
    X[n] = sinf(0.314F*i);
  }

  et = clock();
  r1 = x87_sum(X, n);
  printf("elapsed time of x87 sum: %f\n", (double)(clock() - et)/CLOCKS_PER_SEC);

  et = clock();
  r2 = simd_sum(X, n);
  printf("elapsed time of SSE2 sum: %f\n", (double)(clock() - et)/CLOCKS_PER_SEC);

  printf("the difference of results: %f\n", fabs(r1 - r2));

  return 0;
}