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.
)
{
// Thread identifiers
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]);
// Synchronise before loading the next tile
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?.
