【发布时间】：2014-08-08 12:06:29
【问题描述】：
使用下面的代码，我用 cramer 规则计算逆矩阵 4x4，但是如何将此代码扩展为 NxN 矩阵？
void PIII_Inverse_4x4(float* src) {
__m128 minor0,minor1,minor2,minor3;
__m128 row0,row1,row2,row3;
__m128 det,tmp1;
tmp1= _mm_loadh_pi(_mm_loadl_pi(tmp1, (__m64*)(src)), (__m64*)(src+ 4));
row1= _mm_loadh_pi(_mm_loadl_pi(row1, (__m64*)(src+8)), (__m64*)(src+12));
row0= _mm_shuffle_ps(tmp1, row1, 0x88);
row1= _mm_shuffle_ps(row1, tmp1, 0xDD);
tmp1= _mm_loadh_pi(_mm_loadl_pi(tmp1, (__m64*)(src+ 2)), (__m64*)(src+ 6));
row3= _mm_loadh_pi(_mm_loadl_pi(row3, (__m64*)(src+10)), (__m64*)(src+14));
row2= _mm_shuffle_ps(tmp1, row3, 0x88);
row3= _mm_shuffle_ps(row3, tmp1, 0xDD);
tmp1= _mm_mul_ps(row2, row3);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor0= _mm_mul_ps(row1, tmp1);
minor1= _mm_mul_ps(row0, tmp1);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor0= _mm_sub_ps(_mm_mul_ps(row1, tmp1), minor0);
minor1= _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor1);
minor1= _mm_shuffle_ps(minor1, minor1, 0x4E);
tmp1= _mm_mul_ps(row1, row2);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor0= _mm_add_ps(_mm_mul_ps(row3, tmp1), minor0);
minor3= _mm_mul_ps(row0, tmp1);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor0= _mm_sub_ps(minor0, _mm_mul_ps(row3, tmp1));
minor3= _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor3);
minor3= _mm_shuffle_ps(minor3, minor3, 0x4E);
tmp1= _mm_mul_ps(_mm_shuffle_ps(row1, row1, 0x4E), row3);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0xB1);
row2= _mm_shuffle_ps(row2, row2, 0x4E);
minor0= _mm_add_ps(_mm_mul_ps(row2, tmp1), minor0);
minor2= _mm_mul_ps(row0, tmp1);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor0= _mm_sub_ps(minor0, _mm_mul_ps(row2, tmp1));
minor2= _mm_sub_ps(_mm_mul_ps(row0, tmp1), minor2);
minor2= _mm_shuffle_ps(minor2, minor2, 0x4E);
tmp1= _mm_mul_ps(row0, row1);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor2= _mm_add_ps(_mm_mul_ps(row3, tmp1), minor2);
minor3= _mm_sub_ps(_mm_mul_ps(row2, tmp1), minor3);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor2= _mm_sub_ps(_mm_mul_ps(row3, tmp1), minor2);
minor3= _mm_sub_ps(minor3, _mm_mul_ps(row2, tmp1));
tmp1= _mm_mul_ps(row0, row3);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor1= _mm_sub_ps(minor1, _mm_mul_ps(row2, tmp1));
minor2= _mm_add_ps(_mm_mul_ps(row1, tmp1), minor2);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor1= _mm_add_ps(_mm_mul_ps(row2, tmp1), minor1);
minor2= _mm_sub_ps(minor2, _mm_mul_ps(row1, tmp1));
tmp1= _mm_mul_ps(row0, row2);
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0xB1);
minor1= _mm_add_ps(_mm_mul_ps(row3, tmp1), minor1);
minor3= _mm_sub_ps(minor3, _mm_mul_ps(row1, tmp1));
tmp1= _mm_shuffle_ps(tmp1, tmp1, 0x4E);
minor1= _mm_sub_ps(minor1, _mm_mul_ps(row3, tmp1));
minor3= _mm_add_ps(_mm_mul_ps(row1, tmp1), minor3);
// -----------------------------------------------
// -----------------------------------------------
// -----------------------------------------------
det= _mm_mul_ps(row0, minor0);
det= _mm_add_ps(_mm_shuffle_ps(det, det, 0x4E), det);
det= _mm_add_ss(_mm_shuffle_ps(det, det, 0xB1), det);
tmp1= _mm_rcp_ss(det);
det= _mm_sub_ss(_mm_add_ss(tmp1, tmp1), _mm_mul_ss(det, _mm_mul_ss(tmp1, tmp1)));
det= _mm_shuffle_ps(det, det, 0x00);
minor0 = _mm_mul_ps(det, minor0);
_mm_storel_pi((__m64*)(src), minor0);
_mm_storeh_pi((__m64*)(src+2), minor0);
minor1 = _mm_mul_ps(det, minor1);
_mm_storel_pi((__m64*)(src+4), minor1);
_mm_storeh_pi((__m64*)(src+6), minor1);
minor2 = _mm_mul_ps(det, minor2);
_mm_storel_pi((__m64*)(src+ 8), minor2);
_mm_storeh_pi((__m64*)(src+10), minor2);
minor3 = _mm_mul_ps(det, minor3);
_mm_storel_pi((__m64*)(src+12), minor3);
_mm_storeh_pi((__m64*)(src+14), minor3);
}
我在谷歌上搜索过，但我没有找到任何有用的东西......我也搜索了逆矩阵的高斯乔丹方法，但没有......
【问题讨论】：
-
英特尔在这里有一个 6x6 的例子：intel.com/design/pentiumiii/sml/245044.htm
-
矩阵的元素不一定是数字，它们也可以是矩阵。例如。您可能会将 16x16 数字矩阵视为具有 4x4 矩阵元素的 4x4 矩阵。然后你可以使用克莱默规则来反转它。将此原理扩展到更大的矩阵是微不足道的。
-
@MaratDukhan：我该怎么做？你认为有可能用 cramer 来扩展这个概念吗？
-
当不存在逆时，这段代码的行为如何？即，如果行列式为零，或者通过病态或精度限制评估为零？
标签： matrix x86 sse simd matrix-inverse