What is the fastest way to perform vector-by-vector dot products for two MxN matrices with small M in CUDA?

I have two matrices, each of which are MxN, where M = 16 and N is much larger (say n = 262144, for example). My goal is to produce a vector of length N, where each element corresponds to the dot product of the nth vector in each of the matrices.

I have attempted the following approach, where cIdx corresponds to the column index of a column vector in each matrix. Not surprisingly, NVIDIA Visual Profiler is telling me that this approach is mostly memory-bandwidth bound.

    public static void MatrixDotProduct(
        float* matrix1,
        float* matrix2,
        float* dotProduct,
        int2 matrixDimensions)
    {
        int i = blockIdx.x * blockDim.x + threadIdx.x;
        int stride = gridDim.x * blockDim.x;
        float sum;

        for (int cIdx = i; cIdx < matrixDimensions.y; cIdx += stride)
        {
            int ci = cIdx * matrixDimensions.x;
            sum = 0f;

            for (int j = 0; j < matrixDimensions.x; j++)
            {
                sum += matrix1[ci + j] * matrix2[ci + j];
            }

            dotProduct[cIdx] = sum;
        }
    }

I have also found a version of a vector-by-vector dot product, which I tried to parallelize as well. Unfortunately, this ran 20% slower than the above implementation (perhaps it's because my M = 16?). Is there a better way to approach this problem that I'm missing here?