- SSE - Обратная матрица с крамером 4x4, Как расширить NxN?

С помощью следующего кода я вычисляю обратную матрицу 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);
}

Я искал в Google, но не нашел ничего полезного... Я также искал метод Гаусса-Джордана для обратной матрицы, но ничего...


x86 matrix simd matrix-inverse sse
person user3661321    schedule 18.06.2014    source источник
comment
У Intel есть пример для 6x6 здесь: intel.com/design/pentiumiii/sml/ 245044.htm   -  person Paul R    schedule 18.06.2014
comment
Элементы матрицы не обязательно являются числами, они также могут быть матрицами. Например. вы можете думать о матрице чисел 16x16 как о матрице 4x4 с элементами матрицы 4x4. Тогда вы можете использовать правило Крамера, чтобы инвертировать его. Тривиально распространить этот принцип на большие матрицы.   -  person Marat Dukhan    schedule 18.06.2014
comment
@MaratDukhan: как я могу делать то, что ты говоришь? как вы думаете, можно ли расширить эту концепцию с помощью Крамера?   -  person user3661321    schedule 18.06.2014
comment
Как ведет себя этот код, когда не существует обратного? т. е. если определитель равен нулю или оценивается как нуль из-за плохой обусловленности или пределов точности?   -  person Brett Hale    schedule 19.06.2014

Ответы (1)


arrow_upward
3
arrow_downward

Для обращения матриц более эффективно будет использовать разложение Холекси, чем метод исключения Гаусса-Жордана. статья Обращение матрицы с использованием разложения Холецкого.

Я реализовал декомпозицию Холекси с помощью SSE, AVX и FMA, а также с несколькими потоками. Я получаю около 50% производительности Intel MKL. В настоящее время я переписываю код своего ядра, поэтому надеюсь вскоре приблизиться к производительности MKL.

Обратите внимание, что инверсия часто не нужна для решения системы уравнений. В моем случае мне нужно только разложение Холецкого, а затем я использую обратную и прямую замену для решения.

person Z boson    schedule 18.06.2014
comment
спасибо, но мне нужно использовать только технологию SIMD, без OPEN MP, AVX и т. д. - person user3661321; 18.06.2014
comment
@ user3661321, AVX — это SIMD. Вы имеете в виду, что вам нужен только SSE. Должно быть легко преобразовать код AVX в SSE. Какой размер N вы имеете в виду? Для меньших значений N (независимо от того, что меньше означает) имеет смысл использовать метод, отличный от разложения Холецкого. - person Z boson; 18.06.2014