OpenCL compute large matrices

I am trying to compute a big matrix $$C$$, of size $$M \times N$$, using the elements of the input matrices $$A$$ and $$B$$. The matrix $$A$$ is of size $$M\times4$$ and the matrix $$B$$ is of the size $$2\times N$$. Both $$M$$ and $$N$$ are positive even integers values.

The $$C$$ matrix is indexed as:

$$C = \begin{bmatrix} C_{0,0} & C_{0,1} & C_{0,2} & \dots & C_{0,N-1} \\ C_{1,0} & C_{1,1} & C_{1,2} & \dots & C_{1,N-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ C_{M-1,0} & C_{M-1,1} & C_{M-1,2} & \dots & C_{M-1,N-1} \end{bmatrix}$$

The elements of $$C$$ are computed from a function $$f$$, using the elements of $$A$$ and $$B$$, as:

$$C_{i,j} = f(A_{i,0},\ A_{i,1},\ A_{i,2},\ A_{i,3}, \ B_{0,j}, \ B_{1,j} )$$

The OpenCL code, for a simple example function $$f$$, is:

CODE =
" #define TS 2 // Workgroup size TS = 2.

float f( float a,
float b,
float c,
float d,
float e,
float f
)
{
float fvalue =  a + b + c + d + e + f;
return fvalue;
}

__kernel void myFun( const __global float* A,   // Input matrix A of size M x 4.
const __global float* B,   // Input matrix B of size 2 x N.
__global float* C          // Output matrix C.
)
{
const int col = get_local_id(0);                  // Local column ID
const int row = get_local_id(1);                  // Local row ID
const int globalCol = TS*get_group_id(0) + col;   // Global column ID
const int globalRow = TS*get_group_id(1) + row;   // Global row ID
const int ColSizeC = get_global_size(0);          // Column size of the B matrix = Column size of the C matrix.

// Declare sub matrices in the local memory of the workgroup
__local float Asub[TS][4];       // Sub-matrix of A with size 2 x 4
__local float Bsub[TS][TS];      // Sub-matrix of B with size 2 x 2

// Initialise the accumulation register
float acc = 0.0f;

// Load the sub-matrix Asub on the local memory
Asub[row][0] = A[globalRow * 4 + 0];
Asub[row][1] = A[globalRow * 4 + 1];
Asub[row][2] = A[globalRow * 4 + 2];
Asub[row][3] = A[globalRow * 4 + 3];

// Load the submatrix Bsub on the local memory
Bsub[0][col] = B[0 * ColSizeC + globalCol];
Bsub[1][col] = B[1 * ColSizeC + globalCol];

// Synchronise to make sure the submatrices are loaded
barrier(CLK_LOCAL_MEM_FENCE);

// Call the function f
acc = f(Asub[row][0], Asub[row][1], Asub[row][2], Asub[row][3], Bsub[0][col], Bsub[1][col]);

barrier(CLK_LOCAL_MEM_FENCE);

// Store the final result in C
C[globalRow*ColSizeC + globalCol] = acc;

}";


To run this code I use:

(* Load the kernel *)
Cmatrix =
OpenCLFunctionLoad[CODE, "myFun", {{"Float", 2, "Input"}, {"Float", 2, "Input"},
{"Float", 2, "Output"}}, {2, 2}, "ShellOutputFunction" -> Print]

(* Define the input matrices and their dimensions *)
Amatrix = RandomReal[{1, 10}, {8, 4}];
RowSizeA =  First[Dimensions[Amatrix]];
Bmatrix = Transpose[RandomReal[{1, 10}, {4, 2}]];
ColSizeB = Last[Dimensions[Bmatrix]];
ColSizeC = ColSizeB;
RowSizeC = RowSizeA;

(* Compute the result *)
Result = Flatten[Cmatrix[Amatrix, Bmatrix, ConstantArray[0, {RowSizeC, ColSizeC}], {ColSizeC, RowSizeC}], 1];


In the code above, I have only used a simple function $$f$$ to keep the question short, while my original function $$f$$ is much more complicated. My problem arises when I want to compute huge matrices of size ColSizeC = 16384, RowSizeC = 1024 or larger. I seem to run out of memory and my system crashes. I have a 64-bit system with 8.00GB RAM and AMD Graphics processor.

To solve the memory issue I could perhaps call the OpenCl function Cmatrix, write a few rows of the matrix to a file, empty the memory somehow on the fly and continue the process untill I populate/ write all the rows to the file. But, I am not sure how to do it. My ultimate aim is to compute the $$C$$ matrix as fast as possible. How can I do this?.

If memory is your bottleneck then don't use a GPU. The transfer of data between CPU (and its main memory) and the GPU is so slow that it virtually always takes longer than performing the computation directly on the CPU. Writing to file from GPU is even much, much slower.

The last time I tried "OpenCLLink" is quite long ago. But if I am not mistaken, Amatrix, Bmatrix, and the ConstantArray are still on CPU's memory when you call Cmatrix, right? Maybe you should load them with OpenCLMemoryLoad to the GPU first? Even better is to generate them directly on the GPU.

An alternative is to use Compile the code and to run it on the CPU; it needs less than 50 ms.

cf = Compile[{{a, _Real, 1}, {B, _Real, 2}},
Table[
Plus[
CompileGetElement[a, 1],
CompileGetElement[a, 2],
CompileGetElement[a, 3],
CompileGetElement[a, 4],

CompileGetElement[B, i, 1],
CompileGetElement[B, i, 2]
],
{i, 1, Length[B]}
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];

Amatrix = RandomReal[{1, 10}, {1024, 4}];
Bmatrix = RandomReal[{1, 10}, {16384, 2}];
r = cf[Amatrix, Bmatrix]; // AbsoluteTiming // First


0.04768

But of course I don't know if that's within your time constraints.

• Thanks for your ideas. I will try what you suggested with the GPU otherwise a compiled function on the CPU would be the best option. – dykes Nov 28 '18 at 12:29

The problem in my OpenCL code was probably coming due to the local memory constraints of the workgroups. After modifying the code a little bit, I could get it to work for the size 16384 x 1024:

CODE =
"
float f( float a,
float b,
float c,
float d,
float e,
float f
)
{
float fvalue =  a + b + c + d + e + f;
return fvalue;
}

__kernel void myFun( const __global float* A,   // Input matrix A of size M x 4.
const __global float* B,   // Input matrix B of size 2 x N.
__global float* C          // Output matrix C.
)
{
int globalCol = get_global _id(0);
int globalRow = get_global _id(1);
const int ColSizeC = get_global _size(0); // Column size of B = Column size of C.

// Find the offset
int AOffset = globalRow << 2;
int BOffset = globalCol << 1;

// Compute the elements of C
C[globalRow*ColSizeC + globalCol] = f(A[AOffset], A[AOffset + 1], A[AOffset + 2], A[AOffset + 3] , B[BOffset], B[BOffset + 1]);

}";


Run the code:

Cmatrix = OpenCLFunctionLoad[CODE, "myFun", {{"Float", 2, "Input"}, {"Float", 2, "Input"},
{"Float", 2, "Output"}}, {16, 16}, "ShellOutputFunction" -> Print]

(* Define the input matrices and their dimensions *)
Amatrix = RandomReal[{1, 10}, {16384, 4}];
RowSizeA =  First[Dimensions[Amatrix]];
Bmatrix = Transpose[RandomReal[{1, 10}, {1024, 2}]];
ColSizeB = Last[Dimensions[Bmatrix]];
ColSizeC = ColSizeB;
RowSizeC = RowSizeA;

(* Compute the result *)
Result = Flatten[Cmatrix[Amatrix, Bmatrix, ConstantArray[0, {RowSizeC,
ColSizeC}], {ColSizeC, RowSizeC}], 1]; // AbsoluteTiming // First
`

0.56711