Блог

Этот набор вложенных циклов for корректно работает для значений M = 64 и N = 64, но не работает, когда я делаю M = 128 и N = 64. У меня есть другая программа, которая проверяет правильность значений для умножения матрицы. Интуитивно кажется, что это все еще должно работать, но дает мне неправильный ответ.

 for(int m=64;m<=M;m =64){
for(int n=64;n<=N;n =64){
    for(int i = m-64; i < m; i =16){

        float *A_column_start, *C_column_start;
        __m128 c_1, c_2, c_3, c_4, a_1, a_2, a_3, a_4, mul_1, 
               mul_2, mul_3, mul_4, b_1;
        int j, k;

        for(j = m-64; j < m; j  ){

            //Load 16 contiguous column aligned elements from matrix C in
            //c_1-c_4 registers

            C_column_start = C i j*M;

            c_1 = _mm_loadu_ps(C_column_start);
            c_2 = _mm_loadu_ps(C_column_start 4);
            c_3 = _mm_loadu_ps(C_column_start 8);
            c_4 = _mm_loadu_ps(C_column_start 12);

            for (k=n-64; k < n; k =2){

                //Load 16 contiguous column aligned elements from matrix A to
                //the a_1-a_4 registers

                A_column_start = A k*M;

                a_1 = _mm_loadu_ps(A_column_start i);
                a_2 = _mm_loadu_ps(A_column_start i 4);
                a_3 = _mm_loadu_ps(A_column_start i 8);
                a_4 = _mm_loadu_ps(A_column_start i 12);

                //Load a value to resgister b_1 to act as a "B" or ("A^T") 
                //element to multiply against the A matrix

                b_1 = _mm_load1_ps(A_column_start j);

                mul_1 = _mm_mul_ps(a_1, b_1);
                mul_2 = _mm_mul_ps(a_2, b_1);
                mul_3 = _mm_mul_ps(a_3, b_1);
                mul_4 = _mm_mul_ps(a_4, b_1);

                //Add together all values of the multiplied A and "B"
                //(or "A^T") matrix elements

                c_4 = _mm_add_ps(c_4, mul_4);
                c_3 = _mm_add_ps(c_3, mul_3);
                c_2 = _mm_add_ps(c_2, mul_2);
                c_1 = _mm_add_ps(c_1, mul_1);

                //Move over one column in A, and load the next 16 contiguous 
                //column aligned elements from matrix A to the a_1-a_4 registers

                A_column_start =M;

                a_1 = _mm_loadu_ps(A_column_start i);
                a_2 = _mm_loadu_ps(A_column_start i 4);
                a_3 = _mm_loadu_ps(A_column_start i 8);
                a_4 = _mm_loadu_ps(A_column_start i 12);

                //Load a value to resgister b_1 to act as a "B" or "A^T"
                //element to multiply against the A matrix

                b_1 = _mm_load1_ps(A_column_start j);

                mul_1 = _mm_mul_ps(a_1, b_1);
                mul_2 = _mm_mul_ps(a_2, b_1);
                mul_3 = _mm_mul_ps(a_3, b_1);
                mul_4 = _mm_mul_ps(a_4, b_1);

                //Add together all values of the multiplied A and "B" or
                //("A^T") matrix elements

                c_4 = _mm_add_ps(c_4, mul_4);
                c_3 = _mm_add_ps(c_3, mul_3);
                c_2 = _mm_add_ps(c_2, mul_2);
                c_1 = _mm_add_ps(c_1, mul_1);

            }
            //Store the added up C values back to memory

            _mm_storeu_ps(C_column_start, c_1);
            _mm_storeu_ps(C_column_start 4, c_2);
            _mm_storeu_ps(C_column_start 8, c_3);
            _mm_storeu_ps(C_column_start 12, c_4);

        }

    }
    }
}}
  

Я предполагаю, что вы используете M в коде

 C_column_start = C i j*M;
  

вместо этого необходимо использовать m . Возможно, также в других строках, где вы используете M .
Однако я не совсем понимаю ваш код, поскольку вы не объяснили, для чего предназначен код, а я не программист-математик.

Это работает корректно для M = 64 и N = 64, потому что в этих случаях вы выполняете только одну итерацию в соответствующем цикле (два самых внешних). Когда у вас будет M = 128, вы теперь выполняете два шага во внешнем цикле, и в этом случае строка

 C_column_start = C i j*M;
  

и строка

 A_column_start = A k*M;
  

Для внутреннего цикла будут получены те же результаты, поэтому, по сути, для двух шагов, которые вы выполняете во внешнем цикле (m = 64 128), вы просто удваиваете результат одного шага с m = 128. Исправление так же просто, как изменение M на m, чтобы вы использовали переменную итерации.

Также вам следует рассмотреть возможность выравнивания ваших данных в A и C, чтобы вы могли выполнять загрузки, выровненные по SSE. Это приведет к гораздо более быстрому кодированию.