逆行列さらに再び

http://www.a.mei.titech.ac.jp/~kabe/calcsoft/library/lapack/clapack/usage/index.htmlを参考にfloatの4x4逆行列をAccelerate Frameworkを使って求めてみる。

#import <Accelerate/Accelerate.h>

void InvMatrix4x4( __CLPK_real *matrix )
{
	__CLPK_integer m,n,lda,info,lwork;
	__CLPK_integer piv[4];
	__CLPK_real work[4];
	lwork = lda = m = n = 4;
	sgetrf_(&m, &n, matrix, &lda, piv, &info);
	sgetri_(&n, matrix, &lda, piv, work, &lwork, &info);
}

どの関数使うのか探すときは目眩がしたけど、上記URLのおかげでさくっとできました。

-(void)test
{
	float mat_prime[] = {
		1, 2, 0, 5,
		0, 1, 0, 0,
		3, 0, 1, 0,
		5, 4, 3, 1
	};
	float matrix[16];
	memcpy( matrix, mat_prime, sizeof(float)*16 );
	InvMatrix4x4(matrix);
	for ( int i = 0; i < 4; ++i ) {
		printf("%.2f %.2f %.2f %.2f\n",matrix[i*4],matrix[i*4+1],matrix[i*4+2],matrix[i*4+3]);
	}
}

次は行列ｘベクターをどう扱うのか探そう。

LU分解を使っていて、行列式が0になる、あるいは0に近い行列では失敗するケースがあります。
infoをチェックすれば、成功したかどうかなどわかるので、必要に応じてエラー処理等は書き足す必要があります。