イントリンシック命令（関数）の例外の原因がわかりません

[たろ]

お世話になります。
SIMDプログラミングを勉強しはじめたのですが、エラーの原因がわからず悩んでいます。
助言いただけると助かります。よろしくお願いします。

環境： WinXP Home SP3, VC2008EE, C言語, Win32API

以下サイトを参考に、ビットマップ画像の加工処理（色反転）の最適化・高速化を試すプログラムを作りました。
http://d.hatena.ne.jp/komugi_com/20080323/1206249192

しかし実行したところ、イントリンシック命令（関数）のところでアクセス違反の例外が発生してしまいます。
ポインタの値が不正になっているようには見えず、原因がわかりません。

長いですが、プログラムを貼らせていただきます。
問題の箇所は、DIB32InverseSIMD 関数の中の、_mm_subs_epu8 です。ここで例外が発生します。
SIMD命令を記述しない、同じ処理を行う関数 DIB32Inverse は問題なく動作しています。

コード:

//
//	>cl inverse.c
//
#pragma comment(lib, "user32.lib")
#pragma comment(lib, "gdi32.lib")
#pragma comment(lib, "shell32.lib")
#pragma comment(lib, "winmm.lib")		// timeGetTime

#include <windows.h>
#include <mmsystem.h>					// timeGetTime
#include <emmintrin.h>					// SSE2
#include <stdio.h>						// _snprintf

// 32bitDIB
typedef struct DIB32		DIB32;
struct DIB32
{
	BITMAPINFOHEADER		BIH;
	BYTE					Pixel[1];
};

//
//	色反転(SIMD)
//
void DIB32InverseSIMD( DIB32* dib )
{
	__m128i*	p;
	__m128i*	pEnd;
	__m128i		i255;

	if( !dib ) return;

	p = (__m128i*)dib->Pixel;
	pEnd = (__m128i*)(dib->Pixel + dib->BIH.biSizeImage);	// 最終ピクセル+1

	i255 = _mm_set1_epi8( 255 );

	for( ; p<pEnd; p++ )			// 先頭から128bit(16byte)ずつ
	{
		*p = _mm_subs_epu8( i255, *p );
	}
}

//
//	色反転(C)
//
void DIB32Inverse( DIB32* dib )
{
	LPBYTE	p, pEnd;

	if( !dib ) return;

	p = dib->Pixel;
	pEnd = dib->Pixel + dib->BIH.biSizeImage;	// 最終ピクセル+1

	for( ; p<pEnd; p+=4 )			// 先頭から1ピクセルずつ
	{
		p[0] = 255 - p[0];	// B
		p[1] = 255 - p[1];	// G
		p[2] = 255 - p[2];	// R
	}
}

//
//	24bitBMPファイルピクセルデータを読み込んで32bit領域に格納する
//
void DIB32LoadBMP24Pixel( DIB32* dib32, HANDLE hFile )
{
	DWORD	lineByte24	= (dib32->BIH.biWidth * 3 + 3) & ~3;
	LPBYTE	lineTop24	= (LPBYTE)HeapAlloc( GetProcessHeap(), 0, lineByte24 );
	LPBYTE	lineEnd24	= lineTop24 + dib32->BIH.biWidth * 3;
	LPBYTE	p32			= dib32->Pixel;
	LPBYTE	p32end		= p32 + dib32->BIH.biSizeImage;

	if( !lineTop24 ) return;

	while( p32 < p32end )
	{
		DWORD	dwRead;
		LPBYTE	p24;
		// １行バッファに読み込んで
		if( !ReadFile(hFile, lineTop24, lineByte24, &dwRead, NULL) || dwRead!=lineByte24 ) break;
		// ピクセルコピー
		for( p24=lineTop24; p24<lineEnd24; p24+=3 )
		{
			p32[0] = p24[0];		 // B
			p32[1] = p24[1];		 // G
			p32[2] = p24[2];		 // R
			p32[3] = 127;			 // A(ピクセル透明度固定値)
			p32 += 4;
		}
	}
	HeapFree( GetProcessHeap(), 0, lineTop24 );
}

//
//	24ビットBMPファイルを読み込んで32ビットDIBを返す
//	BITMAPINFOHEADER構造体とピクセルデータを連続領域で確保して先頭アドレスを返す
//	+------------------+----------------+
//	| BITMAPINFOHEADER | ピクセルデータ |
//	+------------------+----------------+
//
DIB32* DIB32LoadFile( char* fileName )
{
	DIB32*				dib32 = NULL;
	HANDLE				hFile;
	BITMAPFILEHEADER	bfh;
	BITMAPINFOHEADER	bih;
	DWORD				dwRead, dwSize;

	hFile = CreateFile( fileName, GENERIC_READ, FILE_SHARE_READ, NULL, OPEN_EXISTING, FILE_ATTRIBUTE_NORMAL, NULL );
	if( hFile==INVALID_HANDLE_VALUE ) return NULL;
	// BITMAPFILEHEADER読み込み
	dwSize = sizeof( BITMAPFILEHEADER );
	if( !ReadFile(hFile, (LPBYTE)&bfh, dwSize, &dwRead, NULL) || dwRead!=dwSize ) goto ERR;
	// BMPファイルチェック
	if( bfh.bfType!=0x4D42 ) goto ERR;
	// BITMAPINFOHEADER読み込み
	dwSize = sizeof( BITMAPINFOHEADER );
	if( !ReadFile(hFile, (LPBYTE)&bih, dwSize, &dwRead, NULL) || dwRead!=dwSize ) goto ERR;
	// BMP形式チェック(24ビットWindowsボトムアップビットマップのみ)
	if( bih.biSize!=40 || bih.biHeight<0 || bih.biBitCount!=24 ) goto ERR;
	// ヘッダ情報を32bitに書き換え
	bih.biBitCount	= 32;
	bih.biSizeImage = bih.biHeight * bih.biWidth * 4;
	// メモリ確保
	// ピクセルデータはSIMDで128bit(16byte)単位でループ処理するため16の倍数byteにする
	dwSize = sizeof( BITMAPINFOHEADER ) + (bih.biSizeImage + 15) & ~15;
	dib32 = (DIB32*)HeapAlloc( GetProcessHeap(), 0, dwSize );
	if( !dib32 ) goto ERR;
	// ヘッダコピー
	CopyMemory( dib32, &bih, sizeof(BITMAPINFOHEADER) );
	// ピクセルデータ24bit→32bit変換読み込み
	DIB32LoadBMP24Pixel( dib32, hFile );
	// おわり
	CloseHandle(hFile);
	return dib32;
 ERR:
	if( dib32 ) HeapFree( GetProcessHeap(), 0, dib32 );
	if( hFile!=INVALID_HANDLE_VALUE ) CloseHandle( hFile );
	return NULL;
}

void DIB32Destroy( DIB32* dib )
{
	if( dib ) HeapFree( GetProcessHeap(), 0, dib );
}

// メインウィンドウプロシージャ
LRESULT CALLBACK MainWndProc( HWND hwnd, UINT msg, WPARAM wp, LPARAM lp )
{
	static DIB32* dib = NULL;

	switch( msg )
	{
		case WM_DROPFILES:
		{
			// ファイルのドラッグ＆ドロップ
			HDROP	hDrop = (HDROP)wp;
			UINT	nCount = DragQueryFile( hDrop, (UINT)-1, NULL, 0 );
			if( nCount )
			{
				char	fileName[MAX_PATH] = "";
				DWORD	dwTime;
				char	msg[256];
				DIB32*	tmp;
				// ドロップされた１つ目のファイルを読み込む
				DragQueryFile( hDrop, 0, fileName, MAX_PATH );
				tmp = DIB32LoadFile( fileName );
				if( !tmp ) return 0;
				DIB32Destroy( dib );
				dib = tmp;
				// 色反転、時間計る
				dwTime = timeGetTime();
				DIB32Inverse( dib );
				DIB32InverseSIMD( dib );
				dwTime = timeGetTime() - dwTime;
				// 画像表示
				InvalidateRect( hwnd, NULL, TRUE );
				// 処理時間報告
				_snprintf(msg,sizeof(msg),"画像サイズ ＝ %u x %u\r\n処理時間 ＝ %u.%03u秒",
						dib->BIH.biWidth, dib->BIH.biHeight, dwTime/1000, dwTime%1000);
				MessageBox( hwnd, msg, "情報", MB_OK|MB_ICONINFORMATION );
			}
			return 0;
		}
		case WM_PAINT:
		{
			// クライアント領域描画
			PAINTSTRUCT ps;
			HDC			hDC = BeginPaint( hwnd, &ps );
			if( !hDC ) return 0;
			if( dib )
			{
				// 画像をクライアント領域いっぱいに広げて表示(縦横比無視)
				RECT	rc;
				GetClientRect( hwnd, &rc );
				SetStretchBltMode( hDC, HALFTONE );
				StretchDIBits(	hDC,
								0, 0, rc.right, rc.bottom,		// コピー先X,Y,幅,高さ
								0, 0,							// コピー元左上X,Y
								dib->BIH.biWidth,				// コピー元幅
								dib->BIH.biHeight,				// コピー元高さ
								dib->Pixel,						// コピー元ピクセルデータ
								(LPBITMAPINFO)dib,				// コピー元BITMAPINFO構造体アドレス
								DIB_RGB_COLORS, SRCCOPY );		// RGBデータそのままコピー
			}
			else
			{
				#define USAGE "24bitBMPファイルをドロップしてください"
				TextOut( hDC, 10, 10, USAGE, strlen(USAGE) );
			}
			EndPaint( hwnd, &ps );
			return 0;
		}
		case WM_SIZE:
			// ウィンドウサイズ変更
			InvalidateRect( hwnd, NULL, TRUE );
			break;

		case WM_DESTROY:
			// アプリケーション終了
			DIB32Destroy( dib );
			PostQuitMessage(0);
			break;
	}
	return DefWindowProc( hwnd, msg, wp, lp );
}

int WINAPI WinMain( HINSTANCE hInstance, HINSTANCE hPrevInstance, LPSTR pCmdLine, int nCmdShow )
{
	WNDCLASSEX	wc;
	MSG			msg;
	HWND		hwnd;
	// ウィンドウクラス登録
	ZeroMemory( &wc, sizeof(WNDCLASSEX) );
	wc.cbSize			= sizeof(WNDCLASSEX);
	wc.lpfnWndProc		= MainWndProc;
	wc.hInstance		= hInstance;
	wc.hCursor			= LoadCursor( NULL, IDC_ARROW );
	wc.hbrBackground	= (HBRUSH)( COLOR_APPWORKSPACE+1 );
	wc.lpszClassName	= "MainWndClass";
	RegisterClassEx( &wc );
	// ウィンドウ作成
	hwnd = CreateWindowEx(	WS_EX_ACCEPTFILES, "MainWndClass", "メインウィンドウ", WS_OVERLAPPEDWINDOW,
							CW_USEDEFAULT, CW_USEDEFAULT, /*CW_USEDEFAULT*/400, /*CW_USEDEFAULT*/400,
							NULL, (HMENU)NULL, hInstance, NULL );
	// 表示
	ShowWindow( hwnd, nCmdShow );
	UpdateWindow( hwnd );
	// メッセージループ
	for( ;; )
	{
		int ret = GetMessage( &msg, NULL, 0, 0 );
		// WM_QUIT(0)orエラー(-1)なら抜ける
		if( ret==0 || ret==-1 ) break;
		// メッセージ処理
		TranslateMessage( &msg );
		DispatchMessage( &msg );
	}
	return msg.wParam;
}

なお、アセンブラの知識は、書籍「はじめて読む8086」をざっと一読した程度です。
自分でニーモニック記述のコードを書いたり、コンパイラの作る中間ファイルを眺める
といったことはまだできないレベルです・・。

よろしくお願いします。

[たろ]

書き忘れたことがありましたので追記します。

CPUは Core Duo T2300 1.66GHz です。
CPU-Zというソフトで確認すると、SSE, SSE2, SSE3 までサポートしていると出ます。
今回問題の _mm_subs_epu8 は、SSE2 の命令だから使えるはず・・と思っています。

以下の自作SSE2テストプログラムは動作することを確認しています。

コード:

#include <emmintrin.h>
#include <stdio.h>

int main( void )
{
	{
		__m128i a, b, c;

		a = _mm_set_epi32( 255, 255, 255, 255 );
		b = _mm_set_epi32( 1, 5, 50, 100 );
		c = _mm_sub_epi32( a, b );

		printf("%d %d %d %d\n", ((int*)&c)[0], ((int*)&c)[1], ((int*)&c)[2], ((int*)&c)[3] );
	}
	{
		int x[4] = { 255, 255, 255, 255 };
		int y[4] = { 1, 5, 50, 100 };

		*(__m128i*)x = _mm_sub_epi32( *(__m128i*)x, *(__m128i*)y );

		printf("%d %d %d %d\n", x[0], x[1], x[2], x[3]);
	}
	return 0;
}

[たろ]

自己レスです。
いきなりですみませんが、原因というか例外発生のきっかけ？が判明したかもしれません。

HeapAlloc で確保した領域を、_mm_sub_epi32 などのイントリンシック関数に渡してはいけないのでしょうか？
同じような例外が発生するパターンを見つけました。

以下プログラムは、_mm_sub_epi32 で例外が発生してしまいます。

コード:

#include <windows.h>
#include <emmintrin.h>
#include <stdio.h>

int main( void )
{
	int* x = (int*)HeapAlloc( GetProcessHeap(), 0, sizeof(int)*4 );
	int* y = (int*)HeapAlloc( GetProcessHeap(), 0, sizeof(int)*4 );

	x[0]=255, x[1]=255, x[2]=255, x[3]=255;
	y[0]=1,   y[1]=5,   y[2]=50,  y[3]=100;

	*(__m128i*)x = _mm_sub_epi32( *(__m128i*)x, *(__m128i*)y );

	printf("%d %d %d %d\n", x[0], x[1], x[2], x[3]);

	HeapFree( GetProcessHeap(), 0, x );
	HeapFree( GetProcessHeap(), 0, y );
	return 0;
}

HeapAlloc でなく、_mm_malloc を使ってメモリ確保すれば動きます。

コード:

int main( void )
{
	int* x = (int*)_mm_malloc( sizeof(int)*4, sizeof(__m128i) );
	int* y = (int*)_mm_malloc( sizeof(int)*4, sizeof(__m128i) );

	x[0]=255, x[1]=255, x[2]=255, x[3]=255;
	y[0]=1,   y[1]=5,   y[2]=50,  y[3]=100;

	*(__m128i*)x = _mm_sub_epi32( *(__m128i*)x, *(__m128i*)y );

	printf("%d %d %d %d\n", x[0], x[1], x[2], x[3]);

	_mm_free( x );
	_mm_free( y );
	return 0;
}

不思議なのですが、静的な配列なら問題なく動くように見えます。

コード:

int main( void )
{
	int x[4] = { 255, 255, 255, 255 };
	int y[4] = { 1, 5, 50, 100 };

	*(__m128i*)x = _mm_sub_epi32( *(__m128i*)x, *(__m128i*)y );

	printf("%d %d %d %d\n", x[0], x[1], x[2], x[3]);
	return 0;
}

これはいったいどういう事なのでしょうか・・・。

HeapAlloc と _mm_malloc はなにが違うのか？「メモリを16バイトでアライメントする」とは、
私は「16の倍数バイト確保しておく」という意味だと思っていましたが、違うのでしょうか・・？

[御津凪]

少なくとも、SIMD系演算でメモリ上にあるデータを操作する場合は 16 バイトアラインメントになっている必要があります。
16ビットアラインメントとは、操作するデータのあるアドレス値が16の倍数（16進数で1桁目が常に0）である必要があるということです。
それに合わないデータを渡すと例外が発生します。

なので、確保するタイミングによっては HeapAlloc や malloc でも（偶然16バイトアラインメントになって）通ることがあるでしょうし、逆にローカル変数で例外が発生する時があるでしょう。

解決策は、_mm_malloc のように、16バイトアラインメントされた領域を確保してくれる関数を使用して領域を用意するか、
自前で16バイトアラインメントに合わせる様にオフセットして使用するかのどちらかにする必要があります。

[たろ]

御津凪さん、ありがとうございます。

メモリのアライメントとは、アドレス値のことだったんですね。
そんな制約があるとは想像できず、入門サイトにも書いてありましたがスルーしてました・・勉強になりました。

HeapAllocでも、15バイト余分に確保した上で、16バイト境界のアドレスを使うようにしたら、ちゃんと動作しました。

ただ、_mm_malloc や _aligned_malloc に相当する、HeapAlloc 系のWin32APIは存在しないんでしょうか・・。
自前でオフセットするのも面倒なので、_mm_malloc を使っておきたいと思います。

また、ローカル変数が問題なかったのは、たまたま16バイト境界のアドレス値になっていたからなんですね。
危ないので、_declspec(align(16)) というのをつけて定義すべき、というのも知りました。

解決とさせて頂きます。ありがとうございました。

[ISLe]

ただ、_mm_malloc や _aligned_malloc に相当する、HeapAlloc 系のWin32APIは存在しないんでしょうか・・。
自前でオフセットするのも面倒なので、_mm_malloc を使っておきたいと思います。

VirtualAllocはページ単位でメモリを確保するAPIでページ境界(64KB)にアラインメントされます。
HeapAllocはもちろん、_mm_malloc や _aligned_malloc も内部でVirtualAllocを呼び出していると思われます。

[たろ]

ISLeさんありがとうございます。遅い返信ですみません。

VirtualAllocは、名前は聞いたことありましたが、アライメントの事は気づきませんでした。
64KBブロック（ページ）単位で確保されるということは、アドレス値も64KBアライメントになる、ということなんですね。

このAPIは、メモリの断片化を防ぐとか自分でメモリの配置まで管理したい人向けの難易度の高いもの・・というイメージで、
使ったことがなかったですが、ピクセルデータをVirtualAllocで確保すると都合がよい、という単純な使い道でもいいのかな。

ありがとうございました。