Inverting a 4x4 matrix

Question

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.

Language ideally C++, data is available in array of 16 floats in column-major order.

Answer 1

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 MESA implementation of the GLU library.

Answer 2

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 program that calculate the determinant and inverse of any N-matrix.

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.

Take a look about the program here.

EDIT: The matrix layout is row-by-row (meaning m01 is in the first row and second column). Also the language is C#, but should be easy to convert into C.