SFITSIO-1.4 / SLLIB-1.4 の紹介

2013年4月 HSCチーム 談話会 資料

山内＠国立天文台 天文データセンター

本談話会の内容

CREDITS

SLLIB，SFITSIO は宇宙研の赤外の皆様，宇宙研のデータセンターの皆様の 貴重な協力と支援のもので開発されてきたソフトウェアです．

特に宇宙研のデータセンターでは，これらのソフトウェアの価値に 多大なるご理解をいただき， 優秀な人材と多額な資金を惜しみなく SLLIB と SFITSIO の開発 (試験・マニュアル作成・翻訳) に投入してくださいました．

感謝を表明したいと思います．

イントロダクション --- SFITSIO，SLLIB とは?

データセンター・解析プロジェクトなどでの業務や， 真面目にデータ解析を行ないたい天文学者のための， 「IDL/Pythonのようなプログラミング言語環境」を提供する， C/C++ のライブラリ・セット．

メモリ管理，配列の選択・コピーなどの 「相対的に科学と距離のある部分」に関する数々の無駄な労力を除去し， プログラマが 科学の本質部分のコーディングに集中でき， 「正しさ」と「透明さ」を保ちつつ， 高いパフォーマンスを狙える開発環境 の整備を目的としている．

SFITSIO，SLLIB の方向性

「C++」の文字を見た瞬間に「そんな難しいもの使えるわけない」 と誤解する人が，かなりの確率で観測される．

イントロダクション --- SFITSIO，SLLIB 開発の動機

開発のキッカケ(問題提起)

ならば，IDL/Python みたいなプログラミング環境をC++上に作ってしまえ:

SFITSIO，SLLIB の開発方針

どんなライブラリを目指しているか?

技術的な方針

最新安定版 SFITSIO-1.4，SLLIB-1.4 の現状

　Ver.1.2→『ただの FITS I/O』
　Ver.1.4→『IDL/Python風のデータ解析向きの開発環境』

SFITSIO-1.4，SLLIB-1.4 に共通する特徴

SFITSIO-1.4 の主な特徴

SLLIB-1.4 の主な特徴

ソフトウェア品質

速度性能

SLLIB，SFITSIOを利用・検討しているプロジェクト

宇宙研ではすでにアーカイブの業務や衛星プロジェクトでの基幹システムにて， 採用実績を有する．

SFITSIO，SLLIB の代表的なメリット

SFITSIO, SLLIB を使うと，開発手順はこうなる

即，本質的なところからコーディングを開始

SLLIB を使ったコードの例

ストリーム

ヘッダ 

#include <sli/stdstreamio.h>
#include <sli/digeststreamio.h>
他

使用例

文字列

ヘッダ 

#include <sli/tstring.h>

使用例

文字列配列

ヘッダ 

#include <sli/tarray_tstring.h>

使用例

多次元配列

内部では単純な1次元バッファを持つ． メンバ関数の引数は 0-indexed である．

基底クラス「mdarray」

型も配列長も可変だが， 型に依存したメンバ関数が使えない．

継承クラス「mdarray_int」「mdarray_float」「mdarray_double」など

通常はこちらを使う． 型が固定である事以外のデメリットはない． 異なる型どおしでも，演算が可能．

ヘッダ 

#include <sli/mdarray.h>             /* mdarrayクラスとその継承クラス */
#include <sli/mdarray_statistics.h>  /* 統計用の関数 */
#include <sli/mdarray_math.h>        /* 配列に対する数学関数 */

使用例

SFITSIO を使ったコードの例

画像，テーブル，ヘッダ，HDUとも， メンバ関数の引数は常に 0-indexed である (ファイル読み取り時の範囲指定時のみ，1-indexed でも記述可能)．

ヘッダ

#include <sli/fitscc.h>
#include <sli/fits_image_statistics.h>

典型的な使い方

観測フレームの1次処理のコードの例

  fitscc bias_fits, flat_fits, obj1_fits;             /* FITSオブジェクト */
  const char *bsec;
  double bsec_mean;
  bias_fits.read_stream("bias.fits.gz");              /* 処理済みバイアスを読む */
  flat_fits.read_stream("flat.fits.gz");              /* 処理済みフラットを読む */
  obj1_fits.read_stream("obj1.fits");                 /* 観測フレームを読む */
  fits_image &obj1pri = obj1_fits.image(0L);          /* 長いのでaliasを定義 */
  obj1pri.convert_type(FITS::FLOAT_T);                /* データ型を変換する */
  mdarray_float &obj1_array = obj1pri.float_array();  /* 内部配列objectの参照を取得 */
  bsec = obj1pri.header("BIASSEC").svalue();          /* overscan情報を取得 */
  bsec_mean = md_mean(obj1_array.sectionf(bsec));     /* overscan領域の平均値 */
  obj1_array -= bsec_mean;                            /* overscan(平均値)引き */
  obj1_array -= bias_fits.image(0L).data_array();     /* バイアス引き */
  obj1_array /= flat_fits.image(0L).data_array();     /* フラット補正 */
  obj1_fits.write_stream("obj1_done.fits");           /* 保存 */

新規 FITS を作るのに便利な「FITSテンプレート」

「FITSテンプレート」とは，新規 FITS の元となる， FITSヘッダのようなテキストファイルである． SFITSIO は，曖昧な FITS テンプレートから， 正しい FITS を生成してくれる．

ヘッダだけの高速読み取り

データセンターなどでは，大量の FITS ファイルを管理するために， FITS ヘッダを高速に読み取る必要がある．

FITS ヘッダをマネージするための fits_header クラスが持つ read_stream() を使えば， Data Unit を読まずにファイルへのアクセスを終了する事が可能． もちろん，圧縮ファイルでも高速アクセスが可能．

今後の予定

おねがい

まとめ

SIMD命令: イントロ

SIMD命令とは?

天文業界では?

性能は?

SIMD命令: 実際のコード

2バイトのバイトスワップのための関数 s_byteswap2()

処理の流れ: 

inline static void s_byteswap2( void *buf, size_t len_elements )
{
    unsigned char *d_ptr = (unsigned char *)buf;
    unsigned char tmp;

#ifdef USE_OPTIMIZED_BYTE_SWAP
#ifdef _SSE2_IS_OK
    if ( _SSE2_MIN_NBYTES <= 2 * len_elements ) {
	bool is_aligned = ( ((size_t)d_ptr & 0x01) == 0 );
	if ( is_aligned == true ) {			/* align if possible */
	    while ( ((size_t)d_ptr & 0x0f) != 0 ) {
		len_elements --;
		tmp = d_ptr[0];  d_ptr[0] = d_ptr[1];  d_ptr[1] = tmp;
		d_ptr += 2;
	    }
	}
	/* NOTE: Do not use _mm_stream...() that causes slow down.      */
	/*       Using load+store improves slightly, but u+u is enough. */
	while ( 8 <= len_elements ) {
	    len_elements -= 8;
	    __m128i r0 = _mm_loadu_si128((__m128i *)d_ptr);
	    __m128i r1 = r0;
	    /* 1-byte shift */
	    r0 = _mm_srli_epi16(r0, 8);
	    r1 = _mm_slli_epi16(r1, 8);
	    /* or */
	    r0 = _mm_or_si128(r0, r1);
	    _mm_storeu_si128((__m128i *)d_ptr, r0);
	    d_ptr += 2 * 8;
	}
    }
#else
    while ( 8 <= len_elements ) {
	len_elements -= 8;
	tmp = d_ptr[0];  d_ptr[0] = d_ptr[1];  d_ptr[1] = tmp;
	tmp = d_ptr[2];  d_ptr[2] = d_ptr[3];  d_ptr[3] = tmp;
	tmp = d_ptr[4];  d_ptr[4] = d_ptr[5];  d_ptr[5] = tmp;
	tmp = d_ptr[6];  d_ptr[6] = d_ptr[7];  d_ptr[7] = tmp;
	tmp = d_ptr[8];  d_ptr[8] = d_ptr[9];  d_ptr[9] = tmp;
	tmp = d_ptr[10];  d_ptr[10] = d_ptr[11];  d_ptr[11] = tmp;
	tmp = d_ptr[12];  d_ptr[12] = d_ptr[13];  d_ptr[13] = tmp;
	tmp = d_ptr[14];  d_ptr[14] = d_ptr[15];  d_ptr[15] = tmp;
	d_ptr += 2 * 8;
    }
#endif	/* _SSE2_IS_OK */
#endif	/* USE_OPTIMIZED_BYTE_SWAP */

    while ( 0 < len_elements ) {
	len_elements --;
	tmp = d_ptr[0];  d_ptr[0] = d_ptr[1];  d_ptr[1] = tmp;
	d_ptr += 2;
    }

    return;
}

SIMD命令: SIMD利用のための準備とコンパイル

例えば，次のようなコードで SSE2 の利用が可能かを判定． 

#define USE_SIMD 1

#if defined(USE_SIMD) && defined(__SSE2__)
#if (defined(__GNUC__) && __GNUC__ >= 4) || defined(__INTEL_COMPILER)
#define _SSE2_IS_OK 1
#endif
#endif

命令セットに応じたヘッダを include する． 

#include <xmmintrin.h>    /* SSE */
#include <emmintrin.h>    /* SSE2 */
#include <immintrin.h>    /* AVX */

4.0 以降の gcc でコンパイルする． 64-bitの場合はオプションは必要なし． 32-bitの場合は「-msse2」をつける． 

$ gcc -Wall -O2 hoge.c -o hoge
$ g++ -Wall -O2 hoge.cc -o hoge

SIMD命令: ロード命令・ストア命令

ロード命令

メモリから 128-bit レジスタ (XMMレジスタ) に16バイト分転送する．

もし「_mm_load_...」を アラインしてないアドレス(例えば0xbffff4e1)に対して使うと， SEGVで落ちる．

最近のCPUでは，「_mm_load_...」，「_mm_loadu_...」との 性能差はかなり小さくなっている． 可能ならばアラインさせておいて，「_mm_loadu_...」で十分かと．

ストア命令

128-bit レジスタ (XMMレジスタ) からメモリへ16バイト分転送する．

ストア命令も最近のCPUでは，「_mm_store_...」，「_mm_storeu_...」との 性能差はかなり小さくなっている． 可能ならばアラインさせておいて，「_mm_storeu_...」で十分．

ストリーム命令

128-bit レジスタ (XMMレジスタ) からメモリへ16バイト分， キャッシュを経由せずに書き込む． DMA のようなもの．

ロード命令・ストア命令のまとめ

基本的には，今のCPUなら，可能ならアラインさせておいて， 「_mm_loadu_...」 「_mm_storeu_...」 でOK．ストリーム命令については，この後で解説．

SIMD命令: 使い方が難しいストリーム命令

キャッシュに入りきらないようなデータのコピーや初期化などに利用する．

メモリ初期化のための関数 s_memset()

inline static void *s_memset( void *s, int c, size_t n, size_t total_n )
{
    unsigned char *d_ptr = (unsigned char *)s;

#ifdef _SSE2_IS_OK
    if ( _SSE2_MIN_NBYTES <= n ) {	/* use SSE2 if n is large enough */
	const unsigned char src[16] __attribute__((aligned(16))) = 
					  {c,c,c,c, c,c,c,c, c,c,c,c, c,c,c,c};
	size_t aa = ((size_t)d_ptr & 0x0f);
	bool use_mm_stream;
	if ( _SSE2_CPU_CACHE_SIZE <= n || _SSE2_CPU_CACHE_SIZE <= total_n ) {
	    use_mm_stream = true;
	}
	else use_mm_stream = false;
	/* set start */
	if ( 0 < aa ) {			/* align with 16-byte */
	    aa = 16 - aa;
	    n -= aa;
	    memset(d_ptr, c, aa);
	    d_ptr += aa;
	}
	__m128i r0 = _mm_load_si128((__m128i *)src);
	if ( use_mm_stream == true ) {
	    while ( 16 <= n ) {
		n -= 16;
		/* without polluting the cache */
		_mm_stream_si128((__m128i *)d_ptr, r0);
		d_ptr += 16;
	    }
	}
	else {
	    while ( 16 <= n ) {
		n -= 16;
		_mm_store_si128((__m128i *)d_ptr, r0);
		d_ptr += 16;
	    }
	}
    }
#endif	/* _SSE2_IS_OK */

    memset(d_ptr, c, n);

    return s;
}

SIMD命令: どんな演算命令がある?

インテルコンパイラのリファレンス がわかりやすい(AVXは若干説明が手抜きっぽい)．

自分のマシンがどの命令セットに対応しているかを知るには， 

$ cat /proc/cpuinfo

し，flags のところをみる．

SIMD命令: 注意点

SIMD命令: コンパイラによるSIMD命令の埋め込み