C++ 反转 4x4 矩阵
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/1148309/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Inverting a 4x4 matrix
提问by clamp
I am looking for a sample code implementation on how to invert a 4x4 matrix. I know there is Gaussian eleminiation, LU decomposition, etc., but instead of looking at them in detail I am really just looking for the code to do this.
我正在寻找有关如何反转 4x4 矩阵的示例代码实现。我知道有高斯消除、LU 分解等，但我没有详细查看它们，我只是在寻找执行此操作的代码。
Language ideally C++, data is available in array of 16 floats in column-major order.
理想的语言是 C++，数据以列优先顺序排列在 16 个浮点数的数组中。
回答by shoosh
here:
这里：
bool gluInvertMatrix(const double m[16], double invOut[16])
{
double inv[16], det;
int i;
inv[0] = m[5] * m[10] * m[15] -
m[5] * m[11] * m[14] -
m[9] * m[6] * m[15] +
m[9] * m[7] * m[14] +
m[13] * m[6] * m[11] -
m[13] * m[7] * m[10];
inv[4] = -m[4] * m[10] * m[15] +
m[4] * m[11] * m[14] +
m[8] * m[6] * m[15] -
m[8] * m[7] * m[14] -
m[12] * m[6] * m[11] +
m[12] * m[7] * m[10];
inv[8] = m[4] * m[9] * m[15] -
m[4] * m[11] * m[13] -
m[8] * m[5] * m[15] +
m[8] * m[7] * m[13] +
m[12] * m[5] * m[11] -
m[12] * m[7] * m[9];
inv[12] = -m[4] * m[9] * m[14] +
m[4] * m[10] * m[13] +
m[8] * m[5] * m[14] -
m[8] * m[6] * m[13] -
m[12] * m[5] * m[10] +
m[12] * m[6] * m[9];
inv[1] = -m[1] * m[10] * m[15] +
m[1] * m[11] * m[14] +
m[9] * m[2] * m[15] -
m[9] * m[3] * m[14] -
m[13] * m[2] * m[11] +
m[13] * m[3] * m[10];
inv[5] = m[0] * m[10] * m[15] -
m[0] * m[11] * m[14] -
m[8] * m[2] * m[15] +
m[8] * m[3] * m[14] +
m[12] * m[2] * m[11] -
m[12] * m[3] * m[10];
inv[9] = -m[0] * m[9] * m[15] +
m[0] * m[11] * m[13] +
m[8] * m[1] * m[15] -
m[8] * m[3] * m[13] -
m[12] * m[1] * m[11] +
m[12] * m[3] * m[9];
inv[13] = m[0] * m[9] * m[14] -
m[0] * m[10] * m[13] -
m[8] * m[1] * m[14] +
m[8] * m[2] * m[13] +
m[12] * m[1] * m[10] -
m[12] * m[2] * m[9];
inv[2] = m[1] * m[6] * m[15] -
m[1] * m[7] * m[14] -
m[5] * m[2] * m[15] +
m[5] * m[3] * m[14] +
m[13] * m[2] * m[7] -
m[13] * m[3] * m[6];
inv[6] = -m[0] * m[6] * m[15] +
m[0] * m[7] * m[14] +
m[4] * m[2] * m[15] -
m[4] * m[3] * m[14] -
m[12] * m[2] * m[7] +
m[12] * m[3] * m[6];
inv[10] = m[0] * m[5] * m[15] -
m[0] * m[7] * m[13] -
m[4] * m[1] * m[15] +
m[4] * m[3] * m[13] +
m[12] * m[1] * m[7] -
m[12] * m[3] * m[5];
inv[14] = -m[0] * m[5] * m[14] +
m[0] * m[6] * m[13] +
m[4] * m[1] * m[14] -
m[4] * m[2] * m[13] -
m[12] * m[1] * m[6] +
m[12] * m[2] * m[5];
inv[3] = -m[1] * m[6] * m[11] +
m[1] * m[7] * m[10] +
m[5] * m[2] * m[11] -
m[5] * m[3] * m[10] -
m[9] * m[2] * m[7] +
m[9] * m[3] * m[6];
inv[7] = m[0] * m[6] * m[11] -
m[0] * m[7] * m[10] -
m[4] * m[2] * m[11] +
m[4] * m[3] * m[10] +
m[8] * m[2] * m[7] -
m[8] * m[3] * m[6];
inv[11] = -m[0] * m[5] * m[11] +
m[0] * m[7] * m[9] +
m[4] * m[1] * m[11] -
m[4] * m[3] * m[9] -
m[8] * m[1] * m[7] +
m[8] * m[3] * m[5];
inv[15] = m[0] * m[5] * m[10] -
m[0] * m[6] * m[9] -
m[4] * m[1] * m[10] +
m[4] * m[2] * m[9] +
m[8] * m[1] * m[6] -
m[8] * m[2] * m[5];
det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
if (det == 0)
return false;
det = 1.0 / det;
for (i = 0; i < 16; i++)
invOut[i] = inv[i] * det;
return true;
}
This was lifted from MESAimplementation of the GLU library.
这是从GLU 库的MESA实现中解除的。
回答by willnode
If anyone looking for more costumized code and "easier to read", then I got this
如果有人在寻找更多服装化的代码和“更容易阅读”，那么我得到了这个
var A2323 = m.m22 * m.m33 - m.m23 * m.m32 ;
var A1323 = m.m21 * m.m33 - m.m23 * m.m31 ;
var A1223 = m.m21 * m.m32 - m.m22 * m.m31 ;
var A0323 = m.m20 * m.m33 - m.m23 * m.m30 ;
var A0223 = m.m20 * m.m32 - m.m22 * m.m30 ;
var A0123 = m.m20 * m.m31 - m.m21 * m.m30 ;
var A2313 = m.m12 * m.m33 - m.m13 * m.m32 ;
var A1313 = m.m11 * m.m33 - m.m13 * m.m31 ;
var A1213 = m.m11 * m.m32 - m.m12 * m.m31 ;
var A2312 = m.m12 * m.m23 - m.m13 * m.m22 ;
var A1312 = m.m11 * m.m23 - m.m13 * m.m21 ;
var A1212 = m.m11 * m.m22 - m.m12 * m.m21 ;
var A0313 = m.m10 * m.m33 - m.m13 * m.m30 ;
var A0213 = m.m10 * m.m32 - m.m12 * m.m30 ;
var A0312 = m.m10 * m.m23 - m.m13 * m.m20 ;
var A0212 = m.m10 * m.m22 - m.m12 * m.m20 ;
var A0113 = m.m10 * m.m31 - m.m11 * m.m30 ;
var A0112 = m.m10 * m.m21 - m.m11 * m.m20 ;
var det = m.m00 * ( m.m11 * A2323 - m.m12 * A1323 + m.m13 * A1223 )
- m.m01 * ( m.m10 * A2323 - m.m12 * A0323 + m.m13 * A0223 )
+ m.m02 * ( m.m10 * A1323 - m.m11 * A0323 + m.m13 * A0123 )
- m.m03 * ( m.m10 * A1223 - m.m11 * A0223 + m.m12 * A0123 ) ;
det = 1 / det;
return new Matrix4x4() {
m00 = det * ( m.m11 * A2323 - m.m12 * A1323 + m.m13 * A1223 ),
m01 = det * - ( m.m01 * A2323 - m.m02 * A1323 + m.m03 * A1223 ),
m02 = det * ( m.m01 * A2313 - m.m02 * A1313 + m.m03 * A1213 ),
m03 = det * - ( m.m01 * A2312 - m.m02 * A1312 + m.m03 * A1212 ),
m10 = det * - ( m.m10 * A2323 - m.m12 * A0323 + m.m13 * A0223 ),
m11 = det * ( m.m00 * A2323 - m.m02 * A0323 + m.m03 * A0223 ),
m12 = det * - ( m.m00 * A2313 - m.m02 * A0313 + m.m03 * A0213 ),
m13 = det * ( m.m00 * A2312 - m.m02 * A0312 + m.m03 * A0212 ),
m20 = det * ( m.m10 * A1323 - m.m11 * A0323 + m.m13 * A0123 ),
m21 = det * - ( m.m00 * A1323 - m.m01 * A0323 + m.m03 * A0123 ),
m22 = det * ( m.m00 * A1313 - m.m01 * A0313 + m.m03 * A0113 ),
m23 = det * - ( m.m00 * A1312 - m.m01 * A0312 + m.m03 * A0112 ),
m30 = det * - ( m.m10 * A1223 - m.m11 * A0223 + m.m12 * A0123 ),
m31 = det * ( m.m00 * A1223 - m.m01 * A0223 + m.m02 * A0123 ),
m32 = det * - ( m.m00 * A1213 - m.m01 * A0213 + m.m02 * A0113 ),
m33 = det * ( m.m00 * A1212 - m.m01 * A0212 + m.m02 * A0112 ),
};
I don't write the code, but my program did. I made a small program to make a programthat calculate the determinant and inverse of any N-matrix.
我不写代码，但我的程序写了。我做了一个小程序来制作一个计算任何N矩阵的行列式和逆的程序。
I do it because once in the past I need a code that inverses 5x5 matrix, but nobody in the earth have done this so I made one.
我这样做是因为过去有一次我需要一个逆 5x5 矩阵的代码，但地球上没有人这样做过，所以我做了一个。
Take a look about the program here.
看看这里的程序。
EDIT: The matrix layout is row-by-row (meaning m01is in the first row and second column). Also the language is C#, but should be easy to convert into C.
编辑：矩阵布局是逐行的（意思m01是在第一行和第二列）。语言也是 C#，但应该很容易转换成 C。
回答by willnode
If you need a C++ matrix library with a lot of functions, have a look at Eigen library - http://eigen.tuxfamily.org
如果您需要具有很多功能的 C++ 矩阵库，请查看 Eigen 库 - http://eigen.tuxfamily.org
回答by willnode
I 'rolled up' the MESA implementation (also wrote a couple of unit tests to ensure it actually works).
我“汇总”了 MESA 实现（还编写了几个单元测试以确保它确实有效）。
Here:
这里：
float invf(int i,int j,const float* m){
int o = 2+(j-i);
i += 4+o;
j += 4-o;
#define e(a,b) m[ ((j+b)%4)*4 + ((i+a)%4) ]
float inv =
+ e(+1,-1)*e(+0,+0)*e(-1,+1)
+ e(+1,+1)*e(+0,-1)*e(-1,+0)
+ e(-1,-1)*e(+1,+0)*e(+0,+1)
- e(-1,-1)*e(+0,+0)*e(+1,+1)
- e(-1,+1)*e(+0,-1)*e(+1,+0)
- e(+1,-1)*e(-1,+0)*e(+0,+1);
return (o%2)?inv : -inv;
#undef e
}
bool inverseMatrix4x4(const float *m, float *out)
{
float inv[16];
for(int i=0;i<4;i++)
for(int j=0;j<4;j++)
inv[j*4+i] = invf(i,j,m);
double D = 0;
for(int k=0;k<4;k++) D += m[k] * inv[k*4];
if (D == 0) return false;
D = 1.0 / D;
for (int i = 0; i < 16; i++)
out[i] = inv[i] * D;
return true;
}
I wrote a little about this and display the pattern of positive/negative factors on my blog.
我写了一些关于这个的文章，并在我的博客上展示了积极/消极因素的模式。
As suggested by @LiraNuna, on many platforms hardware accelerated versions of such routines are available so I'm happy to have a 'backup version' that's readable and concise.
正如@LiraNuna 所建议的那样，在许多平台上都可以使用此类例程的硬件加速版本，因此我很高兴拥有一个可读且简洁的“备份版本”。
Note: this may run 3.5 times slower or worse than the MESA implementation. You can shift the pattern of factors to remove some additions etc... but it would lose in readability and still won't be very fast.
注意：这可能比 MESA 实现慢 3.5 倍或更差。您可以改变因子的模式以删除一些添加等......但它会降低可读性并且仍然不会很快。
回答by Svante
You can use the GNU Scientific Libraryor look the code up in it.
您可以使用GNU Scientific Library或在其中查找代码。
Edit: You seem to want the Linear Algebrasection.
编辑：您似乎想要线性代数部分。
回答by Eugene
Here is a small (just one header) C++ vector mathlibrary (geared towards 3D programming). If you use it, keep in mind that layout of its matrices in memory is inverted comparing to what OpenGL expects, I had fun time figuring it out...
这是一个小的（只有一个头文件）C++矢量数学库（面向 3D 编程）。如果你使用它，请记住它在内存中的矩阵布局与 OpenGL 期望的相比是颠倒的，我很开心地弄清楚了......
回答by Samuel Li
Inspired by @shoosh to check out MESA implementations, I found that matrix inversion looks quite different in more recent mesa releases. I suppose those are good improvements. Here's the matrix inversion code from Mesa-17.3.9:
受@shoosh 的启发，查看 MESA 实现，我发现矩阵求逆在最近的 mesa 版本中看起来大不相同。我想这些都是很好的改进。这是来自Mesa-17.3.9的矩阵求逆代码：
/* Returns true for success, false for failure (singular matrix) */
bool DirectVolumeRenderer::_mesa_invert_matrix_general( GLfloat out[16], const GLfloat in[16] )
{
/**
* References an element of 4x4 matrix.
* Calculate the linear storage index of the element and references it.
*/
#define MAT(m,r,c) (m)[(c)*4+(r)]
/**
* Swaps the values of two floating point variables.
*/
#define SWAP_ROWS(a, b) { GLfloat *_tmp = a; (a)=(b); (b)=_tmp; }
const GLfloat *m = in;
GLfloat wtmp[4][8];
GLfloat m0, m1, m2, m3, s;
GLfloat *r0, *r1, *r2, *r3;
r0 = wtmp[0], r1 = wtmp[1], r2 = wtmp[2], r3 = wtmp[3];
r0[0] = MAT(m,0,0), r0[1] = MAT(m,0,1),
r0[2] = MAT(m,0,2), r0[3] = MAT(m,0,3),
r0[4] = 1.0, r0[5] = r0[6] = r0[7] = 0.0,
r1[0] = MAT(m,1,0), r1[1] = MAT(m,1,1),
r1[2] = MAT(m,1,2), r1[3] = MAT(m,1,3),
r1[5] = 1.0, r1[4] = r1[6] = r1[7] = 0.0,
r2[0] = MAT(m,2,0), r2[1] = MAT(m,2,1),
r2[2] = MAT(m,2,2), r2[3] = MAT(m,2,3),
r2[6] = 1.0, r2[4] = r2[5] = r2[7] = 0.0,
r3[0] = MAT(m,3,0), r3[1] = MAT(m,3,1),
r3[2] = MAT(m,3,2), r3[3] = MAT(m,3,3),
r3[7] = 1.0, r3[4] = r3[5] = r3[6] = 0.0;
/* choose pivot - or die */
if (fabsf(r3[0])>fabsf(r2[0])) SWAP_ROWS(r3, r2);
if (fabsf(r2[0])>fabsf(r1[0])) SWAP_ROWS(r2, r1);
if (fabsf(r1[0])>fabsf(r0[0])) SWAP_ROWS(r1, r0);
if (0.0F == r0[0])
return false;
/* eliminate first variable */
m1 = r1[0]/r0[0]; m2 = r2[0]/r0[0]; m3 = r3[0]/r0[0];
s = r0[1]; r1[1] -= m1 * s; r2[1] -= m2 * s; r3[1] -= m3 * s;
s = r0[2]; r1[2] -= m1 * s; r2[2] -= m2 * s; r3[2] -= m3 * s;
s = r0[3]; r1[3] -= m1 * s; r2[3] -= m2 * s; r3[3] -= m3 * s;
s = r0[4];
if (s != 0.0F) { r1[4] -= m1 * s; r2[4] -= m2 * s; r3[4] -= m3 * s; }
s = r0[5];
if (s != 0.0F) { r1[5] -= m1 * s; r2[5] -= m2 * s; r3[5] -= m3 * s; }
s = r0[6];
if (s != 0.0F) { r1[6] -= m1 * s; r2[6] -= m2 * s; r3[6] -= m3 * s; }
s = r0[7];
if (s != 0.0F) { r1[7] -= m1 * s; r2[7] -= m2 * s; r3[7] -= m3 * s; }
/* choose pivot - or die */
if (fabsf(r3[1])>fabsf(r2[1])) SWAP_ROWS(r3, r2);
if (fabsf(r2[1])>fabsf(r1[1])) SWAP_ROWS(r2, r1);
if (0.0F == r1[1])
return false;
/* eliminate second variable */
m2 = r2[1]/r1[1]; m3 = r3[1]/r1[1];
r2[2] -= m2 * r1[2]; r3[2] -= m3 * r1[2];
r2[3] -= m2 * r1[3]; r3[3] -= m3 * r1[3];
s = r1[4]; if (0.0F != s) { r2[4] -= m2 * s; r3[4] -= m3 * s; }
s = r1[5]; if (0.0F != s) { r2[5] -= m2 * s; r3[5] -= m3 * s; }
s = r1[6]; if (0.0F != s) { r2[6] -= m2 * s; r3[6] -= m3 * s; }
s = r1[7]; if (0.0F != s) { r2[7] -= m2 * s; r3[7] -= m3 * s; }
/* choose pivot - or die */
if (fabsf(r3[2])>fabsf(r2[2])) SWAP_ROWS(r3, r2);
if (0.0F == r2[2])
return false;
/* eliminate third variable */
m3 = r3[2]/r2[2];
r3[3] -= m3 * r2[3], r3[4] -= m3 * r2[4],
r3[5] -= m3 * r2[5], r3[6] -= m3 * r2[6],
r3[7] -= m3 * r2[7];
/* last check */
if (0.0F == r3[3])
return false;
s = 1.0F/r3[3]; /* now back substitute row 3 */
r3[4] *= s; r3[5] *= s; r3[6] *= s; r3[7] *= s;
m2 = r2[3]; /* now back substitute row 2 */
s = 1.0F/r2[2];
r2[4] = s * (r2[4] - r3[4] * m2), r2[5] = s * (r2[5] - r3[5] * m2),
r2[6] = s * (r2[6] - r3[6] * m2), r2[7] = s * (r2[7] - r3[7] * m2);
m1 = r1[3];
r1[4] -= r3[4] * m1, r1[5] -= r3[5] * m1,
r1[6] -= r3[6] * m1, r1[7] -= r3[7] * m1;
m0 = r0[3];
r0[4] -= r3[4] * m0, r0[5] -= r3[5] * m0,
r0[6] -= r3[6] * m0, r0[7] -= r3[7] * m0;
m1 = r1[2]; /* now back substitute row 1 */
s = 1.0F/r1[1];
r1[4] = s * (r1[4] - r2[4] * m1), r1[5] = s * (r1[5] - r2[5] * m1),
r1[6] = s * (r1[6] - r2[6] * m1), r1[7] = s * (r1[7] - r2[7] * m1);
m0 = r0[2];
r0[4] -= r2[4] * m0, r0[5] -= r2[5] * m0,
r0[6] -= r2[6] * m0, r0[7] -= r2[7] * m0;
m0 = r0[1]; /* now back substitute row 0 */
s = 1.0F/r0[0];
r0[4] = s * (r0[4] - r1[4] * m0), r0[5] = s * (r0[5] - r1[5] * m0),
r0[6] = s * (r0[6] - r1[6] * m0), r0[7] = s * (r0[7] - r1[7] * m0);
MAT(out,0,0) = r0[4]; MAT(out,0,1) = r0[5],
MAT(out,0,2) = r0[6]; MAT(out,0,3) = r0[7],
MAT(out,1,0) = r1[4]; MAT(out,1,1) = r1[5],
MAT(out,1,2) = r1[6]; MAT(out,1,3) = r1[7],
MAT(out,2,0) = r2[4]; MAT(out,2,1) = r2[5],
MAT(out,2,2) = r2[6]; MAT(out,2,3) = r2[7],
MAT(out,3,0) = r3[4]; MAT(out,3,1) = r3[5],
MAT(out,3,2) = r3[6]; MAT(out,3,3) = r3[7];
#undef SWAP_ROWS
#undef MAT
return true;
}
Note: you can find this piece of code in the mesa code base: mesa-17.3.9/src/mesa/math/m_matrix.c.
注意：你可以在 mesa 代码库中找到这段代码：mesa-17.3.9/src/mesa/math/m_matrix.c.
回答by kstn
This is the C++ version for @willnode's answer
这是@willnode 答案的 C++ 版本
static inline void InvertMatrix4(const Matrix& m, Matrix& im, double& det)
{
double A2323 = m(2, 2) * m(3, 3) - m(2, 3) * m(3, 2);
double A1323 = m(2, 1) * m(3, 3) - m(2, 3) * m(3, 1);
double A1223 = m(2, 1) * m(3, 2) - m(2, 2) * m(3, 1);
double A0323 = m(2, 0) * m(3, 3) - m(2, 3) * m(3, 0);
double A0223 = m(2, 0) * m(3, 2) - m(2, 2) * m(3, 0);
double A0123 = m(2, 0) * m(3, 1) - m(2, 1) * m(3, 0);
double A2313 = m(1, 2) * m(3, 3) - m(1, 3) * m(3, 2);
double A1313 = m(1, 1) * m(3, 3) - m(1, 3) * m(3, 1);
double A1213 = m(1, 1) * m(3, 2) - m(1, 2) * m(3, 1);
double A2312 = m(1, 2) * m(2, 3) - m(1, 3) * m(2, 2);
double A1312 = m(1, 1) * m(2, 3) - m(1, 3) * m(2, 1);
double A1212 = m(1, 1) * m(2, 2) - m(1, 2) * m(2, 1);
double A0313 = m(1, 0) * m(3, 3) - m(1, 3) * m(3, 0);
double A0213 = m(1, 0) * m(3, 2) - m(1, 2) * m(3, 0);
double A0312 = m(1, 0) * m(2, 3) - m(1, 3) * m(2, 0);
double A0212 = m(1, 0) * m(2, 2) - m(1, 2) * m(2, 0);
double A0113 = m(1, 0) * m(3, 1) - m(1, 1) * m(3, 0);
double A0112 = m(1, 0) * m(2, 1) - m(1, 1) * m(2, 0);
det = m(0, 0) * ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 )
- m(0, 1) * ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 )
+ m(0, 2) * ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 )
- m(0, 3) * ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 );
det = 1 / det;
im(0, 0) = det * ( m(1, 1) * A2323 - m(1, 2) * A1323 + m(1, 3) * A1223 );
im(0, 1) = det * - ( m(0, 1) * A2323 - m(0, 2) * A1323 + m(0, 3) * A1223 );
im(0, 2) = det * ( m(0, 1) * A2313 - m(0, 2) * A1313 + m(0, 3) * A1213 );
im(0, 3) = det * - ( m(0, 1) * A2312 - m(0, 2) * A1312 + m(0, 3) * A1212 );
im(1, 0) = det * - ( m(1, 0) * A2323 - m(1, 2) * A0323 + m(1, 3) * A0223 );
im(1, 1) = det * ( m(0, 0) * A2323 - m(0, 2) * A0323 + m(0, 3) * A0223 );
im(1, 2) = det * - ( m(0, 0) * A2313 - m(0, 2) * A0313 + m(0, 3) * A0213 );
im(1, 3) = det * ( m(0, 0) * A2312 - m(0, 2) * A0312 + m(0, 3) * A0212 );
im(2, 0) = det * ( m(1, 0) * A1323 - m(1, 1) * A0323 + m(1, 3) * A0123 );
im(2, 1) = det * - ( m(0, 0) * A1323 - m(0, 1) * A0323 + m(0, 3) * A0123 );
im(2, 2) = det * ( m(0, 0) * A1313 - m(0, 1) * A0313 + m(0, 3) * A0113 );
im(2, 3) = det * - ( m(0, 0) * A1312 - m(0, 1) * A0312 + m(0, 3) * A0112 );
im(3, 0) = det * - ( m(1, 0) * A1223 - m(1, 1) * A0223 + m(1, 2) * A0123 );
im(3, 1) = det * ( m(0, 0) * A1223 - m(0, 1) * A0223 + m(0, 2) * A0123 );
im(3, 2) = det * - ( m(0, 0) * A1213 - m(0, 1) * A0213 + m(0, 2) * A0113 );
im(3, 3) = det * ( m(0, 0) * A1212 - m(0, 1) * A0212 + m(0, 2) * A0112 );
}
回答by Cauchy Schwarz
You can make it faster according to this blog.
你可以根据这个博客让它更快。
#define SUBP(i,j) input[i][j]
#define SUBQ(i,j) input[i][2+j]
#define SUBR(i,j) input[2+i][j]
#define SUBS(i,j) input[2+i][2+j]
#define OUTP(i,j) output[i][j]
#define OUTQ(i,j) output[i][2+j]
#define OUTR(i,j) output[2+i][j]
#define OUTS(i,j) output[2+i][2+j]
#define INVP(i,j) invP[i][j]
#define INVPQ(i,j) invPQ[i][j]
#define RINVP(i,j) RinvP[i][j]
#define INVPQ(i,j) invPQ[i][j]
#define RINVPQ(i,j) RinvPQ[i][j]
#define INVPQR(i,j) invPQR[i][j]
#define INVS(i,j) invS[i][j]
#define MULTI(MAT1, MAT2, MAT3) \
MAT3(0,0)=MAT1(0,0)*MAT2(0,0) + MAT1(0,1)*MAT2(1,0); \
MAT3(0,1)=MAT1(0,0)*MAT2(0,1) + MAT1(0,1)*MAT2(1,1); \
MAT3(1,0)=MAT1(1,0)*MAT2(0,0) + MAT1(1,1)*MAT2(1,0); \
MAT3(1,1)=MAT1(1,0)*MAT2(0,1) + MAT1(1,1)*MAT2(1,1);
#define INV(MAT1, MAT2) \
_det = 1.0 / (MAT1(0,0) * MAT1(1,1) - MAT1(0,1) * MAT1(1,0)); \
MAT2(0,0) = MAT1(1,1) * _det; \
MAT2(1,1) = MAT1(0,0) * _det; \
MAT2(0,1) = -MAT1(0,1) * _det; \
MAT2(1,0) = -MAT1(1,0) * _det; \
#define SUBTRACT(MAT1, MAT2, MAT3) \
MAT3(0,0)=MAT1(0,0) - MAT2(0,0); \
MAT3(0,1)=MAT1(0,1) - MAT2(0,1); \
MAT3(1,0)=MAT1(1,0) - MAT2(1,0); \
MAT3(1,1)=MAT1(1,1) - MAT2(1,1);
#define NEGATIVE(MAT) \
MAT(0,0)=-MAT(0,0); \
MAT(0,1)=-MAT(0,1); \
MAT(1,0)=-MAT(1,0); \
MAT(1,1)=-MAT(1,1);
void getInvertMatrix(complex<double> input[4][4], complex<double> output[4][4]) {
complex<double> _det;
complex<double> invP[2][2];
complex<double> invPQ[2][2];
complex<double> RinvP[2][2];
complex<double> RinvPQ[2][2];
complex<double> invPQR[2][2];
complex<double> invS[2][2];
INV(SUBP, INVP);
MULTI(SUBR, INVP, RINVP);
MULTI(INVP, SUBQ, INVPQ);
MULTI(RINVP, SUBQ, RINVPQ);
SUBTRACT(SUBS, RINVPQ, INVS);
INV(INVS, OUTS);
NEGATIVE(OUTS);
MULTI(OUTS, RINVP, OUTR);
MULTI(INVPQ, OUTS, OUTQ);
MULTI(INVPQ, OUTR, INVPQR);
SUBTRACT(INVP, INVPQR, OUTP);
}
This is not a complete implementation because Pmay not be invertible, but you can combine this code with MESA implementation to get a better performance.
这不是一个完整的实现，因为P可能不可逆，但是您可以将此代码与 MESA 实现结合起来以获得更好的性能。
