C++ - 高效计算矢量矩阵乘积

我需要尽可能高效地计算产品向量矩阵。具体来说，给定一个向量s和一个矩阵A，我需要计算s * A。我有一个类Vector，它包装了一个std::vector和一个类Matrix，它也包装了一个std::vector（为了提高效率）。

天真的方法（我目前正在使用的方法）是拥有类似的东西

Vector<T> timesMatrix(Matrix<T>& matrix)
{
    Vector<unsigned int> result(matrix.columns());
    // constructor that does a resize on the underlying std::vector
    for(unsigned int i = 0 ; i < vector.size() ; ++i)
    {
        for(unsigned int j = 0 ; j < matrix.columns() ; ++j)
        {
            result[j] += (vector[i] * matrix.getElementAt(i, j));
            // getElementAt accesses the appropriate entry
            // of the underlying std::vector
        }
    }
    return result;
}

它工作正常，需要近 12000 微秒。请注意，向量s有 499 个元素，而A是499 x 15500 。

下一步是尝试并行化计算：如果我有N线程，那么我可以为每个线程提供向量s的一部分和矩阵A的"对应"行。每个线程将计算一个 499 大小的Vector，最终结果将是它们的条目总和。
首先，在类Matrix中，我添加了一个方法，用于从Matrix中提取一些行并构建一个较小的行：

Matrix<T> extractSomeRows(unsigned int start, unsigned int end)
{
    unsigned int rowsToExtract = end - start + 1;
    std::vector<T> tmp;
    tmp.reserve(rowsToExtract * numColumns);
    for(unsigned int i = start * numColumns ; i < (end+1) * numColumns ; ++i)
    {
        tmp.push_back(matrix[i]);
    }
    return Matrix<T>(rowsToExtract, numColumns, tmp);
}

然后我定义了一个线程例程

void timesMatrixThreadRoutine
    (Matrix<T>& matrix, unsigned int start, unsigned int end, Vector<T>& newRow)
{
    // newRow is supposed to contain the partial result
    // computed by a thread
    newRow.resize(matrix.columns());
    for(unsigned int i = start ; i < end + 1 ; ++i)
    {
        for(unsigned int j = 0 ; j < matrix.columns() ; ++j)
        {
            newRow[j] += vector[i] * matrix.getElementAt(i - start, j);
        }
    }
}

最后，我修改了上面显示的timesMatrix方法的代码：

Vector<T> timesMatrix(Matrix<T>& matrix)
{
    static const unsigned int NUM_THREADS = 4;
    unsigned int matRows = matrix.rows();
    unsigned int matColumns = matrix.columns();
    unsigned int rowsEachThread = vector.size()/NUM_THREADS;
    std::thread threads[NUM_THREADS];
    Vector<T> tmp[NUM_THREADS];
    unsigned int start, end;
    // all but the last thread
    for(unsigned int i = 0 ; i < NUM_THREADS - 1 ; ++i)
    {
        start = i*rowsEachThread;
        end = (i+1)*rowsEachThread - 1;
        threads[i] = std::thread(&Vector<T>::timesMatrixThreadRoutine, this,
            matrix.extractSomeRows(start, end), start, end, std::ref(tmp[i]));
    }
    // last thread
    start = (NUM_THREADS-1)*rowsEachThread;
    end = matRows - 1;
    threads[NUM_THREADS - 1] = std::thread(&Vector<T>::timesMatrixThreadRoutine, this,
        matrix.extractSomeRows(start, end), start, end, std::ref(tmp[NUM_THREADS-1]));
    for(unsigned int i = 0 ; i < NUM_THREADS ; ++i)
    {
        threads[i].join();
    }
    Vector<unsigned int> result(matColumns);
    for(unsigned int i = 0 ; i < NUM_THREADS ; ++i)
    {
        result = result + tmp[i];    // the operator+ is overloaded
    }
    return result;
}

它仍然有效，但现在需要近 30000 微秒，几乎是以前的三倍。

我做错了什么吗？你认为有更好的方法吗？

编辑 - 使用"轻量级"VirtualMatrix

按照 Ilya Ovodov 的建议，我定义了一个包装T* matrixData的类VirtualMatrix，它在构造函数中初始化

VirtualMatrix(Matrix<T>& m)
{
    numRows = m.rows();
    numColumns = m.columns();
    matrixData = m.pointerToData();
    // pointerToData() returns underlyingVector.data();
}

然后有一个方法可以检索矩阵的特定条目：

inline T getElementAt(unsigned int row, unsigned int column)
{
    return *(matrixData + row*numColumns + column);
}

现在执行时间更好（大约 8000 微秒），但也许需要做出一些改进。特别是线程例程现在是

void timesMatrixThreadRoutine
    (VirtualMatrix<T>& matrix, unsigned int startRow, unsigned int endRow, Vector<T>& newRow)
{
    unsigned int matColumns = matrix.columns();
    newRow.resize(matColumns);
    for(unsigned int i = startRow ; i < endRow + 1 ; ++i)
    {
        for(unsigned int j = 0 ; j < matColumns ; ++j)
        {
            newRow[j] += (vector[i] * matrix.getElementAt(i, j));
        }
    }
}

真正慢的部分是嵌套for循环的部分。如果我删除它，结果显然是错误的，但在不到 500 微秒的时间内"计算"。也就是说，现在传递参数几乎不需要时间，而沉重的部分实际上是计算。

据你说，有没有办法让它更快？