MQL5 OpenCL矩阵乘法深度优化

一、为什么朴素OpenCL优化不够

第一篇文章《OpenCL：连接并行世界的桥梁》介绍了MQL5中OpenCL内核与主机程序交互的基础，并用圆周率计算展示了向量数据类型带来的性能提升。但那种优化是朴素的，因为它完全忽略了执行计算的硬件规格。现代GPU与CPU在存储架构、总线宽度、寄存器数量上差异巨大，只有理解这些规格，才能稳定实现明显超越CPU的提速。本文选择两个大型矩阵相乘这一经典案例，因为它在OpenCL文献中被研究得最充分，且能清晰暴露存储模型的优化空间。

二、OpenCL抽象存储模型与硬件映射

OpenCL定义了一套抽象存储模型以保证代码可移植性：全局内存所有计算单元均可访问，如同主机RAM；常量内存是全局内存一部分，供所有工作单元同时读取；本地内存是每个工作组独占的高速暂存器，组内共享、组间隔离；私有内存则仅供单个工作单元使用，通常位于寄存器。实际硬件如AMD Radeon HD 6970的存储结构与此抽象高度相似，但旧卡如HD 4870的本地内存实际被映射为全局内存，这一点后面会重点说明。

本地内存变量可在内核头声明，但动态数组不能在内核主体中声明，必须指定大小：__local float sharedData[64];。私有变量不含指针且未带__local修饰符时默认私有，通常位于寄存器；私有数组则会溢出到片外高延迟存储。寄存器压力是GPU编程的现实约束，因为大量核心挤在有限芯片面积上，寄存器不可能很多。

__kernel void foo( __global float *A ) { /// kernel code }
__kernel void foo( __local float *sharedData ) { }
__kernel void foo( __global float *A )
{
   __local float sharedData[ 64 ];   
}

三、GPU全局内存合并访问机制

GPU全局内存带宽极高但延迟也高，优化首要目标是提高带宽利用率而非降低延迟。假设存储总线宽32字节（256位），访问地址0x00001232的int（4字节）时，硬件会返回所在32字节块的全部8个int，其中7个无用。若16个线程各自独立发请求，共产生16次传输、448字节无用数据；若将它们合并为相干请求，只需3次甚至2次（对齐时）传输。AMD GPU中线程作为波前一部分按SIMD执行，get_global_id(0)排列的16线程刚好是波前四分之一，合并后能高效利用总线。实测相干请求带宽可提升一个数量级。

存储库是数据实际存放单元（通常32位字），串行数据放相邻库，无库冲突。库冲突对本地内存负面影响最大：AMD波前发生本地库冲突时会串行化等待，严重拖慢内核。例外是所有线程访问同一地址时库可广播避免延迟。全局内存冲突影响较小。

四、MQL5串行矩阵乘法基线

我们先在MQL5写朴素CPU三维循环矩阵乘，维度2000时用时约83.7秒（最快循环顺序），而1000维仅10.5秒，印证O(N^3)复杂度。循环顺序对运行时影响达1.73倍，本文取最快的r-cr-c顺序。该代码作为后续所有OpenCL优化的参照基线，CPU最快吞吐约191 MFlops。

//+------------------------------------------------------------------+
//|                                                matr_mul_2dim.mq5 |
//+------------------------------------------------------------------+
#define ROWS1           1000
#define COLSROWS        1000
#define COLS2           1000
float first[ ROWS1  ][ COLSROWS ];
float second[ COLSROWS ][ COLS2 ];
float third[ ROWS1 ][ COLS2 ];
void OnStart()
  {
   MathSrand(GetTickCount());
   genMatrices();
   ArrayInitialize(third,0.0f);
   uint st1=GetTickCount();
   mul();
   double time1=(double)(GetTickCount()-st1)/1000.;
   Print("CPU: time = "+DoubleToString(time1,3)+" s.");
  }
void mul()
  {
   for(int r=0; r<ROWS1; r++)
      for(int cr=0; cr<COLSROWS; cr++)
         for(int c=0; c<COLS2; c++)
            third[r][c]+=first[r][cr]*second[cr][c];
  }

五、OpenCL首次移植与线性缓冲

将算法移植到OpenCL，创建ROWS1*COLS2线程，删除两层外循环，每个线程算一个输出元素。矩阵在GPU中以行优先线性缓冲排列：Matr[row][col]=buff[row*N+col]。首次内核在Radeon HD 4870上反而比CPU慢（12.7秒 vs 9.3秒），因为未做GPU针对性优化。注意REALTYPE需在主机#define和内核#define同步修改；CLExecute是异步的，真正耗时算在CLBufferRead之后。

const string clSrc=
                   "#define COLS2     "+i2s(COLS2)+"                                  \r\n"
                   "#define COLSROWS  "+i2s(COLSROWS)+"                               \r\n"
                   "#define REALTYPE  float                                             \r\n"
                   "__kernel void matricesMul( __global REALTYPE *in1,                  \r\n"
                   "                           __global REALTYPE *in2,                  \r\n"
                   "                           __global REALTYPE *out  )                \r\n"
                   "{                                                                   \r\n"
                   "  int r = get_global_id( 0 );                                       \r\n"
                   "  int c = get_global_id( 1 );                                       \r\n"

六、消除非相干访问与私有累积

首次内核中in2以COLS2为间隔取数，造成非相干访问。修改索引公式让两数组都顺序取数（in2[cr + c*COLSROWS]），并改主机填充为列优先。合并后CPU明显提速，GPU几乎无变化但为后续铺垫。接着引入私有变量sum累积内积，结束后再写全局输出，将全局写延迟隐藏：GPU从12.7秒降至1.3秒，吞吐达12 GFlops。

for( int cr = 0; cr < COLSROWS; cr ++ )
   out[ r * COLS2 + c ] += in1[ r * COLSROWS + cr ] * in2[ cr + c * COLSROWS ];

__kernel void matricesMul( __global REALTYPE *in1, __global REALTYPE *in2, __global REALTYPE *out  )
{
  int r = get_global_id( 0 );
  int c = get_global_id( 1 );
  REALTYPE sum = 0.0;
  for( int cr = 0; cr < COLSROWS; cr ++ )
     sum += in1[ r * COLSROWS + cr ] * in2[ cr + c * COLSROWS ];
  out[ r * COLS2 + c ] = sum;
}

七、整行计算与私有/本地内存迁移

将任务空间从二维（每线程一元素）降为一维（每线程整行），减少间接开销。再把第一矩阵行复制到私有数组rowbuf，避免主循环重复读全局；然后将第二矩阵列由工作组内线程并行搬入__local colbuf，并用barrier(CLK_LOCAL_MEM_FENCE)同步，确保填充完再算。CPU上本地内存因Intel缓存而无益，HD 4870因本地即全局亦无提升，但此模式是标准优化范式。

__kernel void matricesMul( __global REALTYPE *in1, __global REALTYPE *in2, __global REALTYPE *out )
{
  int r = get_global_id( 0 );
  REALTYPE rowbuf[ COLSROWS ];
  for( int col = 0; col < COLSROWS; col ++ )
     rowbuf[ col ] = in1[ r * COLSROWS + col ];
  int idlocal = get_local_id( 0 );
  int nlocal = get_local_size( 0 );
  __local REALTYPE colbuf[ COLSROWS ];
  // ... parallel copy ...
  barrier(CLK_LOCAL_MEM_FENCE);
}

八、内核向量化与最终性能

最后用float4再float8向量化内积。float4初步让GPU再快10%；将类型转换移出主循环后GPU吞吐飙至32 GFlops。但float8在HD 4870上因CL_DEVICE_PREFERRED_VECTOR_WIDTH_FLOAT=4反而略慢，说明须查首选向量宽度。整体看，经优化的GPU内核比朴素CPU代码快约200倍（CPU 115秒 vs GPU 0.5秒量级）。MQL5当前OpenCL API不能显式设工作组大小，限制了一部分极限性能。

inline REALTYPE dot8( REALTYPE8 a, REALTYPE8 b )
{
  REALTYPE8  c = a * b;
  REALTYPE4  _1 = ( REALTYPE4 ) 1.;
  return( dot( c.lo + c.hi, _1 ) );
}
// kernel with float8 rowbuf and dot8 in loop