第一个改良
和我们的第一个 CUDA 程序一样，我们可以利用 shared memory 来储存每个 row 的数据。不过，因为只有同一个 block 的 threads 可以共享 shared memory，因此现在一个 row 只能由同一个 block 的 threads 来进行计算。另外我们也需要能存放一整个 row 的 shared memory。因此，把先把呼叫 kernel 的部份改成：
matMultCUDA<<>>
(ac, n, bc, n, cc, n, n);
kernel 的部份则改成：
const float* b, size_t ldb, float* c, size_t ldc, int n)
{
extern __shared__ float data[];
const int tid = threadIdx.x;
const int row = blockIdx.x;
int i, j;
for(i = tid; i < n; i += blockDim.x) {
data[i] = a[row * lda + i];
}
__syncthreads();
for(j = tid; j < n; j += blockDim.x) {
float t = 0;
float y = 0;
for(i = 0; i < n; i++) {
float r;
y -= data[i] * b[i * ldb + j];
r = t - y;
y = (r - t) + y;
t = r;
}
c[row * ldc + j] = t;
}
}
第一个部份先把整个 row 读到 shared memory 中，而第二个部份则进行计算，并没有太大的变化。主要的差别是现在一个 row 只由一个 block 进行计算。
在 GeForce 8800GT 上，执行的结果是：
Max error: 1.19209e-007 Average error: 4.22751e-008
Time used: 0.4220 (4.74 GFLOPS)
很明显的，计算的结果并没有改变，不过速度则提高了超过一倍。虽然如此，但是这样的效率仍不尽理想，因为理论上 GeForce 8800GT 有超过 300GFLOPS 的运算性能。即使是把 Kahan's Summation Formula 所需要的额外运算考虑进去，这样的效率仍然连理论最大值的十分之一都不到。
会有这样的结果，原因其实还是同样的：对内存的存取次数太多了。虽然现在 A 矩阵的 row 的数据已经不再需要重复读取，但是 B 矩阵的 column 的数据仍然一直被重复读取。
另一个问题比较不是那么明显：对 B 矩阵的读取，虽然看起来不连续，但实际上它是连续的。这是因为不同的 thread 会读取不同的 column，因此同时间每个 thread 读取的各个 column 加起来，就是一个连续的内存区块。那么，为什么效率还是不佳呢?这是因为，GPU 上的内存控制器，从某个固定的倍数地址开始读取，才会有最高的效率(例如 16 bytes 的倍数)。由于矩阵大小并不是 16 的倍数(这里使用的是 1000x1000 的矩阵)，所以造成效率不佳的情形。
要解决这个问题，我们可以在 cudaMalloc 的时候稍微修改一下，让宽度变成 适当的倍数就可以了。但是，适当的倍数是多少呢?幸运的是，我们并不需要知道这些细节。CUDA 提供了一个 cudaMallocPitch 的函式，可以自动以非常好的的倍数来配置内存。因此，我们可以把 cudaMalloc 的部份改成：
cudaMallocPitch((void**) &ac, &pitch_a, sizeof(float) * n, n);
cudaMallocPitch((void**) &bc, &pitch_b, sizeof(float) * n, n);
cudaMallocPitch((void**) &cc, &pitch_c, sizeof(float) * n, n);
cudaMallocPitch 函式会以适当的倍数配置内存，并把配置的宽度传回。因此，在把矩阵复制到显卡内存上时，要使用它传回的宽度：
sizeof(float) * n, n, cudaMemcpyHostToDevice);
cudaMemcpy2D(bc, pitch_b, b, sizeof(float) * ldb,
sizeof(float) * n, n, cudaMemcpyHostToDevice);
呼叫 kernel 的部份也需要修改：
matMultCUDA<<>>
(ac, pitch_a / sizeof(float), bc, pitch_b / sizeof(float),
cc, pitch_c / sizeof(float), n);
同样的，把计算结果复制回到主内存时，也要使用传回的宽度值：
cudaMemcpy2D(c, sizeof(float) * ldc, cc, pitch_c,
sizeof(float) * n, n, cudaMemcpyDeviceToHost);
这样就完成了。Kernel 部份则不需要修改。
这样的修改有多大的效果呢?在 GeForce 8800GT 上执行，结果如下：
Max error: 1.19209e-007 Average error: 4.22751e-008
Time used: 0.1250 (16.00 GFLOPS)
可以看到，执行速度又再大幅提高了三倍多!而这只是把内存的配置方式稍微修改一下而已。
虽然执行速度提高了很多，但是，和前面提到的理论值相比，其实还是有相当的差距。这是因为，前面也提到过，这样的做法需要 n3+n2 次的内存读取，和 n2 次的内存写入动作。由于 n = 1000，每个数字的大小是 32 bits，所以总共的内存存取数据量约为 4GB。除以实际执行的时间 0.125 秒，得到的带宽数值是约 32GB/s，这已经接近 GeForce 8800GT 显卡内存的带宽了。由于我们计算时间的时候，把配置内存、以及数据的复制动作也计算进去，因此实际上花费在 kernel 的时间是更短的(约 0.09 秒)。因此，可以很明显的看出，这个程序的效率，是受限于内存带宽的。
