Q:

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.
                       ) 
  {
     // 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?.    

A:

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[
    Compile`GetElement[a, 1],
    Compile`GetElement[a, 2],
    Compile`GetElement[a, 3],
    Compile`GetElement[a, 4],

    Compile`GetElement[B, i, 1],
    Compile`GetElement[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.

A:

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.
                    ) 
{
 // Thread identifiers
 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

