- SSE - Matrix inverse with cramer 4x4, How do extends NxN?

Question

With the follow code, I calculate the inverse matrix 4x4 with cramer rules, but how do extend this code for NxN matrix?

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);
}

I searched on google, but I have not found anything useful... I have searched also the gauss-jordan method for inverse matrix, but nothing...

Intel has an example for 6x6 here: http://www.intel.com/design/pentiumiii/sml/245044.htm — Paul R, Jun 18 '14 at 09:12
Elements of a matrix are not necessarily numbers, they could be matrices as well. E.g. you may think of 16x16 matrix of numbers as a 4x4 matrix with 4x4 matrix elements. Then you could use Cramer's rule to invert it. It is trivial to extend this principle to larger matrices. — Marat Dukhan, Jun 18 '14 at 09:44
@MaratDukhan: how can I do what you say? do you think it is possible to extend this concept with cramer? — user3661321, Jun 18 '14 at 10:17
How does this code behave when no inverse exists? i.e., if the determinant is zero, or evaluates as zero through ill-conditioning or precision limits? — Brett Hale, Jun 19 '14 at 06:30

score 3 · Answer 1 · edited May 23 '17 at 11:57

3

For matrix inversion it's going to be more efficient to use Choleksy decomposition than Gauss–Jordan elimination See this paper Matrix Inversion Using Cholesky Decomposition.

I implemented Choleksy decomposition with SSE, AVX, and FMA as well as with multiple threads. I get about 50% of the performance of the Intel MKL. I'm currently rewritting my kernel code so I hope to get closer to MKL performance soon.

Note that the inversion is often unnecessary to solve a system of equations. In my case I only need the Cholesky decomposition and then I use backward and forward substitution to solve.

edited May 23 '17 at 11:57

Community

1
1

answered Jun 18 '14 at 09:41

Z boson

32,619
11
123
226

thanks, but I need use only SIMD technology, without OPEN MP,AVX, etc. – user3661321 Jun 18 '14 at 10:19
3

@user3661321, AVX is SIMD. You mean you need SSE only. It should be easy to convert the AVX code to SSE. What size N are you thinking of? For smaller values of N (whatever smaller means) it make make sense to use a different method than Cholesky Decomposition. – Z boson Jun 18 '14 at 10:38

- SSE - Matrix inverse with cramer 4x4, How do extends NxN?

1 Answers1