C code for 4x4 matrix inversion

C C++ code for row major and column major matrix inversion - 02/2013 - #Jumble, Math

Just leaving some code here to invert either column or row major 4x4 matrices.

The procedure below is designed for column major matrices and comes from the MESA implementation of the glu library. You can find it by downloading the glu-x.x.x.tar.bz2 archive in the file 'src/libutil/project.c'. You'll find also other interesting implems of gluPerspective() gluLookAt() etc.

bool invertColumnMajor(float m[16], float invOut[16])
{
    float 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.f / det;

    for(i = 0; i &lt; 16; i++)
        invOut[i] = inv[i] * det;

    return true;
}

Re-wrote indices to handle row major memory layout (untested):

bool invertRowMajorMESA(float m[16], float invOut[16])
{
    float inv[16], det;
    int i;
 
    inv[ 0] =  m[5] * m[10] * m[15] - m[5] * m[14] * m[11] - m[6] * m[9] * m[15] + m[6] * m[13] * m[11] + m[7] * m[9] * m[14] - m[7] * m[13] * m[10];
    inv[ 1] = -m[1] * m[10] * m[15] + m[1] * m[14] * m[11] + m[2] * m[9] * m[15] - m[2] * m[13] * m[11] - m[3] * m[9] * m[14] + m[3] * m[13] * m[10];
    inv[ 2] =  m[1] * m[ 6] * m[15] - m[1] * m[14] * m[ 7] - m[2] * m[5] * m[15] + m[2] * m[13] * m[ 7] + m[3] * m[5] * m[14] - m[3] * m[13] * m[ 6];
    inv[ 3] = -m[1] * m[ 6] * m[11] + m[1] * m[10] * m[ 7] + m[2] * m[5] * m[11] - m[2] * m[ 9] * m[ 7] - m[3] * m[5] * m[10] + m[3] * m[ 9] * m[ 6];
    inv[ 4] = -m[4] * m[10] * m[15] + m[4] * m[14] * m[11] + m[6] * m[8] * m[15] - m[6] * m[12] * m[11] - m[7] * m[8] * m[14] + m[7] * m[12] * m[10];
    inv[ 5] =  m[0] * m[10] * m[15] - m[0] * m[14] * m[11] - m[2] * m[8] * m[15] + m[2] * m[12] * m[11] + m[3] * m[8] * m[14] - m[3] * m[12] * m[10];
    inv[ 6] = -m[0] * m[ 6] * m[15] + m[0] * m[14] * m[ 7] + m[2] * m[4] * m[15] - m[2] * m[12] * m[ 7] - m[3] * m[4] * m[14] + m[3] * m[12] * m[ 6];
    inv[ 7] =  m[0] * m[ 6] * m[11] - m[0] * m[10] * m[ 7] - m[2] * m[4] * m[11] + m[2] * m[ 8] * m[ 7] + m[3] * m[4] * m[10] - m[3] * m[ 8] * m[ 6];
    inv[ 8] =  m[4] * m[ 9] * m[15] - m[4] * m[13] * m[11] - m[5] * m[8] * m[15] + m[5] * m[12] * m[11] + m[7] * m[8] * m[13] - m[7] * m[12] * m[ 9];
    inv[ 9] = -m[0] * m[ 9] * m[15] + m[0] * m[13] * m[11] + m[1] * m[8] * m[15] - m[1] * m[12] * m[11] - m[3] * m[8] * m[13] + m[3] * m[12] * m[ 9];
    inv[10] =  m[0] * m[ 5] * m[15] - m[0] * m[13] * m[ 7] - m[1] * m[4] * m[15] + m[1] * m[12] * m[ 7] + m[3] * m[4] * m[13] - m[3] * m[12] * m[ 5];
    inv[11] = -m[0] * m[ 5] * m[11] + m[0] * m[ 9] * m[ 7] + m[1] * m[4] * m[11] - m[1] * m[ 8] * m[ 7] - m[3] * m[4] * m[ 9] + m[3] * m[ 8] * m[ 5];
    inv[12] = -m[4] * m[ 9] * m[14] + m[4] * m[13] * m[10] + m[5] * m[8] * m[14] - m[5] * m[12] * m[10] - m[6] * m[8] * m[13] + m[6] * m[12] * m[ 9];
    inv[13] =  m[0] * m[ 9] * m[14] - m[0] * m[13] * m[10] - m[1] * m[8] * m[14] + m[1] * m[12] * m[10] + m[2] * m[8] * m[13] - m[2] * m[12] * m[ 9];
    inv[14] = -m[0] * m[ 5] * m[14] + m[0] * m[13] * m[ 6] + m[1] * m[4] * m[14] - m[1] * m[12] * m[ 6] - m[2] * m[4] * m[13] + m[2] * m[12] * m[ 5];
    inv[15] =  m[0] * m[ 5] * m[10] - m[0] * m[ 9] * m[ 6] - m[1] * m[4] * m[10] + m[1] * m[ 8] * m[ 6] + m[2] * m[4] * m[ 9] - m[2] * m[ 8] * m[ 5];
 
    det = m[0] * inv[0] + m[4] * inv[1] + m[8] * inv[2] + m[12] * inv[3];
 
    if(det == 0)
        return false;
 
    det = 1.f / det;
 
    for(i = 0; i < 16; i++)
        invOut[i] = inv[i] * det;
 
    return true;
}

Alternate implementation

Here we implement the Laplace expansion technique to recursively compute the inverse.Use this routine to invert a row major matrix:

float MINOR(float m[16], int r0, int r1, int r2, int c0, int c1, int c2)
{
    return m[4*r0+c0] * (m[4*r1+c1] * m[4*r2+c2] - m[4*r2+c1] * m[4*r1+c2]) -
           m[4*r0+c1] * (m[4*r1+c0] * m[4*r2+c2] - m[4*r2+c0] * m[4*r1+c2]) +
           m[4*r0+c2] * (m[4*r1+c0] * m[4*r2+c1] - m[4*r2+c0] * m[4*r1+c1]);
}


void adjoint(float m[16], float adjOut[16])
{
    adjOut[ 0] =  MINOR(m,1,2,3,1,2,3); adjOut[ 1] = -MINOR(m,0,2,3,1,2,3); adjOut[ 2] =  MINOR(m,0,1,3,1,2,3); adjOut[ 3] = -MINOR(m,0,1,2,1,2,3);
    adjOut[ 4] = -MINOR(m,1,2,3,0,2,3); adjOut[ 5] =  MINOR(m,0,2,3,0,2,3); adjOut[ 6] = -MINOR(m,0,1,3,0,2,3); adjOut[ 7] =  MINOR(m,0,1,2,0,2,3);
    adjOut[ 8] =  MINOR(m,1,2,3,0,1,3); adjOut[ 9] = -MINOR(m,0,2,3,0,1,3); adjOut[10] =  MINOR(m,0,1,3,0,1,3); adjOut[11] = -MINOR(m,0,1,2,0,1,3);
    adjOut[12] = -MINOR(m,1,2,3,0,1,2); adjOut[13] =  MINOR(m,0,2,3,0,1,2); adjOut[14] = -MINOR(m,0,1,3,0,1,2); adjOut[15] =  MINOR(m,0,1,2,0,1,2);
}

float det(float m[16])
{
    return m[0] * MINOR(m, 1, 2, 3, 1, 2, 3) -
           m[1] * MINOR(m, 1, 2, 3, 0, 2, 3) +
           m[2] * MINOR(m, 1, 2, 3, 0, 1, 3) -
           m[3] * MINOR(m, 1, 2, 3, 0, 1, 2);
}


void invertRowMajor(float m[16], float invOut[16])
{
    adjoint(m, invOut);

    float inv_det = 1.0f / det(m);
    for(int i = 0; i < 16; ++i)
        invOut[i] = invOut[i] * inv_det;
}

4x4 matrix inversion with SSE / SIMD

More complex but faster implementation is available here:

// Copy paste of the code in the pdf.
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);
}

If you need to invert larger matrices I recommend using Eigen.

HLSL

Inversion routines for HLSL shader code

float3x3 MatrixTo3x3(float4x3 InMat)
{
    return float3x3(
        InMat._11_12_13,
        InMat._21_22_23,
        InMat._31_32_33);
}

float3x3 MatrixTo3x3(float4x4 InMat)
{
    return float3x3(
        InMat._11_12_13,
        InMat._21_22_23,
        InMat._31_32_33);
}

float4x3 MatrixTo4x3(float4x4 InMat)
{
    return float4x3(
        InMat._11_12_13,
        InMat._21_22_23,
        InMat._31_32_33,
        InMat._41_42_43);
}

float4x4 MatrixTo4x4(float3x3 InMat)
{
    return float4x4(
        float4(InMat._11_12_13, 0),
        float4(InMat._21_22_23, 0),
        float4(InMat._31_32_33, 0),
        float4(0, 0, 0, 1));
}

float4x4 MatrixTo4x4(float4x3 InMat)
{
    return float4x4(
        float4(InMat._11_12_13, 0),
        float4(InMat._21_22_23, 0),
        float4(InMat._31_32_33, 0),
        float4(InMat._41_42_43, 1));
}

float3x3 MatrixFromScale(float3 InScaleXYZ)
{
    return float3x3(
		InScaleXYZ.x, 0, 0,
		0, InScaleXYZ.y, 0,
		0, 0, InScaleXYZ.z);
}

void MatrixDecompose(float3x3 InMat, inout float3x3 OutRotation, inout float3 OutScaleXYZ)
{
    OutScaleXYZ = float3(
        length(InMat._11_12_13),
        length(InMat._21_22_23),
        length(InMat._31_32_33));

    OutScaleXYZ.x *= determinant(InMat) < 0 ? -1 : 1;
    
    OutRotation = float3x3(
        InMat._11_12_13 / OutScaleXYZ.x,
        InMat._21_22_23 / OutScaleXYZ.y,
        InMat._31_32_33 / OutScaleXYZ.z);
}

float2x2 MatrixAdjoint(float2x2 InMat)
{
    return float2x2(
         InMat._22, 
        -InMat._12, 
        -InMat._21, 
         InMat._11);
}

float3x3 MatrixAdjoint(float3x3 InMat)
{
    return float3x3(
         determinant(float2x2(InMat._22_23, InMat._32_33)),
        -determinant(float2x2(InMat._12_13, InMat._32_33)),
         determinant(float2x2(InMat._12_13, InMat._22_23)),
        -determinant(float2x2(InMat._21_23, InMat._31_33)),
         determinant(float2x2(InMat._11_13, InMat._31_33)),
        -determinant(float2x2(InMat._11_13, InMat._21_23)),
         determinant(float2x2(InMat._21_22, InMat._31_32)),
        -determinant(float2x2(InMat._11_12, InMat._31_32)),
         determinant(float2x2(InMat._11_12, InMat._21_22)));
}

float4x4 MatrixAdjoint(float4x4 InMat)
{
    return float4x4(
         determinant(float3x3(InMat._22_23_24, InMat._32_33_34, InMat._42_43_44)), 
        -determinant(float3x3(InMat._12_13_14, InMat._32_33_34, InMat._42_43_44)),
         determinant(float3x3(InMat._12_13_14, InMat._22_23_24, InMat._42_43_44)),
        -determinant(float3x3(InMat._12_13_14, InMat._22_23_24, InMat._32_33_34)),
        -determinant(float3x3(InMat._21_23_24, InMat._31_33_34, InMat._41_43_44)),
         determinant(float3x3(InMat._11_13_14, InMat._31_33_34, InMat._41_43_44)),
        -determinant(float3x3(InMat._11_13_14, InMat._21_23_24, InMat._41_43_44)),
         determinant(float3x3(InMat._11_13_14, InMat._21_23_24, InMat._31_33_34)),
         determinant(float3x3(InMat._21_22_24, InMat._31_32_34, InMat._41_42_44)),
        -determinant(float3x3(InMat._11_12_14, InMat._31_32_34, InMat._41_42_44)),
         determinant(float3x3(InMat._11_12_14, InMat._21_22_24, InMat._41_42_44)),
        -determinant(float3x3(InMat._11_12_14, InMat._21_22_24, InMat._31_32_34)),
        -determinant(float3x3(InMat._21_22_23, InMat._31_32_33, InMat._41_42_43)),
         determinant(float3x3(InMat._11_12_13, InMat._31_32_33, InMat._41_42_43)),
        -determinant(float3x3(InMat._11_12_13, InMat._21_22_23, InMat._41_42_43)),
         determinant(float3x3(InMat._11_12_13, InMat._21_22_23, InMat._31_32_33)));
}

float2x2 MatrixInverse(float2x2 InMat)
{
    return MatrixAdjoint(InMat) / determinant(InMat);
}

float3x3 MatrixInverse(float3x3 InMat)
{
    return MatrixAdjoint(InMat) / determinant(InMat);
}

// Full matrix inversion, no assumption done
// (slowest)
float4x4 MatrixInverse(float4x4 InMat)
{
    return MatrixAdjoint(InMat) / determinant(InMat);
}

float3x3 MatrixInverse_Affine(float3x3 InMat)
{
    float3 Scale; // out
    float3x3 Rotation; // out
    MatrixDecompose(InMat, Rotation, Scale);
    
    return mul(transpose(Rotation), MatrixFromScale(rcp(Scale)));
}

// warning: assumes left-hand multiplication xA = b
// (last row is the translation)
float4x3 MatrixInverse_Affine(float4x3 InMat)
{
    float3x3 Inverse3x3 = MatrixInverse_Affine(MatrixTo3x3(InMat));
    float3 PosInverse = mul(-InMat._41_42_43, Inverse3x3);
        
    return float4x3(Inverse3x3, PosInverse);
}

// Optimized inversion for when the matrix is an affine transformation.
// which for instance our last column is always 0 0 0 1
// warning: assumes left-hand multiplication xA = b
// (last row is the translation)
float4x4 MatrixInverse_Affine(float4x4 InMat)
{
    return MatrixTo4x4(MatrixInverse_Affine(MatrixTo4x3(InMat)));
}

float3x3 MatrixInverse_Orthonormal(float3x3 InMat)
{
    return transpose(InMat);
}

// Optimized inversion for when there is no scale and shearpresent - 
// only pure translation and rotation.
// (i.e. rigid deformation)
// warning: assumes left-hand multiplication xA = b
// (last row is the translation)
float4x3 MatrixInverse_Orthonormal(float4x3 InMat)
{
    float3x3 Inverse3x3 = MatrixInverse_Orthonormal(MatrixTo3x3(InMat));
    float3 PosInverse = mul(-InMat._41_42_43, Inverse3x3);
        
    return float4x3(Inverse3x3, PosInverse);
}

// warning: assumes left-hand multiplication xA = b
// (last row is the translation)
float4x4 MatrixInverse_Orthonormal(float4x4 InMat)
{
    return MatrixTo4x4(MatrixInverse_Orthonormal(MatrixTo4x3(InMat)));
}

No comments

Rodolphe Vaillant's homepage

Research, teaching and more...

C code for 4x4 matrix inversion

Alternate implementation

4x4 matrix inversion with SSE / SIMD

HLSL