Как исправить проблему в умножении матриц на OpenCL?

Question

[tj57]

Делаю умножение двух матриц на OpenCL, используя векторный тип данных. Опирался на этот пример: https://cnugteren.github.io/tutorial/pages/page6.html. В моём случае размер локальной группы DX = 8, размер векторного типа WIDTH = 4. Умножаю 2 квадратные матрицы 16 х 16.

#define DX 8
#define WIDTH 4

kernel void mul(global const int4 *a, global const int4 *b, global int4 *c, int l, int m, int n) {

    size_t i = get_global_id(0);
    size_t j = get_global_id(1);

	const int row = get_local_id(0); 
	const int col = get_local_id(1);

	const int globalRow = (DX/WIDTH)*get_group_id(0) + row; 
	const int globalCol = DX*get_group_id(1) + col; 

	__local int4 A_part[DX][DX/WIDTH];
    __local int4 B_part[DX][DX/WIDTH];

	const int numTiles = l/DX;

	int4 value = {
		0, 0, 0, 0 
	};

	int localIdx = DX/4;

	for (int i = 0; i < numTiles; i++) {

		const int localRow = DX/WIDTH*i + row;
		const int localCol = DX*i + col;
		A_part[col][row] = a[localCol*(l/WIDTH) + globalRow];
		B_part[col][row] = b[globalCol*(m/WIDTH) + localRow];

		barrier(CLK_LOCAL_MEM_FENCE);

		int4 A_vector, B_vector;
        int valB;

		for (int k=0; k < DX / WIDTH; k++) {

           B_vector = B_part[col][k];

		   for (int w=0; w<WIDTH; w++) {

				A_vector = A_part[WIDTH*k + w][row];
                 switch (w) {
		case 0: valB = B_vector.x; break;
                case 1: valB = B_vector.y; break;
                case 2: valB = B_vector.z; break;
                case 3: valB = B_vector.w; break;
                }
				value.x += A_vector.x * valB;
				value.y += A_vector.y * valB;
				value.z += A_vector.z * valB;
				value.w += A_vector.w * valB;
		   }
		   
        }
		
		barrier(CLK_LOCAL_MEM_FENCE);
	}

    c[globalCol*(m/WIDTH) + globalRow] = value;
}

Такой вариант успешно и верно работает. Теперь мне нужно улучшить алгоритм, чтобы не было внутреннего цикла по WIDTH, так как это лишние вычисления. Я сделал union из матрицы А, чтобы при её записи в локальную память обращаться к ней как к векторному типу, а при вычислениях - как к скалярному. Не могу до конца понять, как правильно обращаться к матрице А при непосредственном умножении:

#define DX 8
#define WIDTH 4

kernel void mul(global const int4 *a, global const int4 *b, global int4 *c, int l, int m, int n) {

    size_t i = get_global_id(0);
    size_t j = get_global_id(1);

	const int row = get_local_id(0); 
	const int col = get_local_id(1);

	const int globalRow = (DX/WIDTH)*get_group_id(0) + row; 
	const int globalCol = DX*get_group_id(1) + col; 

	union {
		int A_part[DX][DX]; 
		int4 A_part_vector[DX][DX/WIDTH];
	} A_union;

    __local int4 B_part[DX][DX/WIDTH];

	const int numTiles = l/DX;

	int4 value = {
		0, 0, 0, 0 
	};

	int localIdx = DX/4;

	for (int i = 0; i < numTiles; i++) {

		const int localRow = DX/WIDTH*i + row;
		const int localCol = DX*i + col;
		A_union.A_part_vector[col][row] = a[localCol*(l/WIDTH) + globalRow]; // обращение к А как к векторному типу
		B_part[col][row] = b[globalCol*(m/WIDTH) + localRow];

		barrier(CLK_LOCAL_MEM_FENCE);

		int4 B_vector;
        int valB;

		for (int k=0; k < DX; k++) { // теперь цикл не по DX/WIDTH, а только по DX

           B_vector = B_part[col][k];

		   // внутренний цикл по ширине убран

			value.x +=  A_union.A_part[WIDTH][row] * B_vector.x; // обращение к А как к скалярному типу
			value.y +=  A_union.A_part[WIDTH][row] * B_vector.y;
			value.z +=  A_union.A_part[WIDTH][row] * B_vector.z;
			value.w +=  A_union.A_part[WIDTH][row] * B_vector.w;
		   
        }
		
		barrier(CLK_LOCAL_MEM_FENCE);
	}

    c[globalCol*(m/WIDTH) + globalRow] = value;
}

Код компилируется, но выводит неверные значения умноженной матрицы:



Как исправить ошибку в перемножении матрицы ?