【发布时间】：2015-03-19 16:35:40
【问题描述】：
我正在对串行和 Cilk 数组表示法版本的矩阵乘法性能进行基准测试。 Cilk 实现所需的时间几乎是串行的两倍，我不明白为什么。
这是单核执行
这是 Cilk 乘法的精髓。由于排名限制，我必须在设置目标矩阵元素值之前将每个乘法存储在一个数组中，然后 __sec_reduce_add 这个数组。
int* sum = new int[VEC_SIZE];
for (int r = 0; (r < rowsThis); r++) {
for (int c = 0; (c < colsOther); c++) {
sum[0:VEC_SIZE] = myData[r][0:VEC_SIZE] * otherData[0:VEC_SIZE][c];
product.data[r][c] = __sec_reduce_add(sum[0:VEC_SIZE]);
}
}
我了解缓存问题，并且认为 Cilk 版本的缓存命中次数少于串行版本的任何原因，因为它们都访问希望在缓存中的列数组以及一系列未命中的行元素。
是否存在我忽略的明显依赖项或我应该使用的 Cilk 句法元素？我应该以不同的方式使用 Cilk 来实现使用 SIMD 操作的非块矩阵乘法的最大性能吗？
我是 Cilk 的新手，因此欢迎任何帮助/建议。
编辑：
这是串行实现：
for (int row = 0; (row < rowsThis); row++) {
for (int col = 0; (col < colsOther); col++) {
int sum = 0;
for (int i = 0; (i < VEC_SIZE); i++) {
sum += (matrix1[row][i] * matrix2[i][col]);
}
matrix3[row][col] = sum;
}
}
内存被适当地（取消）分配，并且乘法结果对于两种实现都是正确的。矩阵大小在编译时是未知的，并且在整个基准测试中使用了多种矩阵大小。
编译器：icpc (ICC) 15.0.0 20140723
编译标志：icpc -g Wall -O2 -std=c++11
忽略分配内存的使用、从向量到普通数组的转换等。我凭直觉破解了另一个程序来运行它，假设它比事实证明的要简单。我无法让编译器在二维向量的两个维度上都接受 cilk 表示法，所以我决定使用传统数组来进行基准测试。
以下是所有适当的文件：
矩阵测试.cpp
#include <string>
#include <fstream>
#include <stdlib.h>
#include "Matrix.h"
#define MATRIX_SIZE 2000
#define REPETITIONS 1
int main(int argc, char** argv) {
auto init = [](int row, int col) { return row + col; };
const Matrix matrix1(MATRIX_SIZE, MATRIX_SIZE, init);
const Matrix matrix2(MATRIX_SIZE, MATRIX_SIZE, init);
for (size_t i = 0; (i < REPETITIONS); i++) {
const Matrix matrix3 = matrix1 * matrix2;
std::cout << "Diag sum: " << matrix3.sumDiag() << std::endl;
}
return 0;
}
矩阵.h
#ifndef MATRIX_H
#define MATRIX_H
#include <iostream>
#include <functional>
#include <vector>
using TwoDVec = std::vector<std::vector<int>>;
class Matrix {
public:
Matrix();
~Matrix();
Matrix(int rows, int cols);
Matrix(int rows, int cols, std::function<Val(int, int)> init);
Matrix(const Matrix& src);
Matrix(Matrix&& src);
Matrix operator*(const Matrix& src) const throw(std::exception);
Matrix& operator=(const Matrix& src);
int sumDiag() const;
protected:
int** getRawData() const;
private:
TwoDVec data;
};
#endif
矩阵.cpp
#include <iostream>
#include <algorithm>
#include <stdexcept>
#include "Matrix.h"
#if defined(CILK)
Matrix
Matrix::operator*(const Matrix& other) const throw(std::exception) {
int rowsOther = other.data.size();
int colsOther = rowsOther > 0 ? other.data[0].size() : 0;
int rowsThis = data.size();
int colsThis = rowsThis > 0 ? data[0].size() : 0;
if (colsThis != rowsOther) {
throw std::runtime_error("Invalid matrices for multiplication.");
}
int** thisRaw = this->getRawData(); // held until d'tor
int** otherRaw = other.getRawData();
Matrix product(rowsThis, colsOther);
const int VEC_SIZE = colsThis;
for (int r = 0; (r < rowsThis); r++) {
for (int c = 0; (c < colsOther); c++) {
product.data[r][c] = __sec_reduce_add(thisRaw[r][0:VEC_SIZE]
* otherRaw[0:VEC_SIZE][c]);
}
}
delete[] thisRaw;
delete[] otherRaw;
return product;
}
#elif defined(SERIAL)
Matrix
Matrix::operator*(const Matrix& other) const throw(std::exception) {
int rowsOther = other.data.size();
int colsOther = rowsOther > 0 ? other.data[0].size() : 0;
int rowsThis = data.size();
int colsThis = rowsThis > 0 ? data[0].size() : 0;
if (colsThis != rowsOther) {
throw std::runtime_error("Invalid matrices for multiplication.");
}
int** thisRaw = this->getRawData(); // held until d'tor
int** otherRaw = other.getRawData();
Matrix product(rowsThis, colsOther);
const int VEC_SIZE = colsThis;
for (int r = 0; (r < rowsThis); r++) {
for (int c = 0; (c < colsOther); c++) {
int sum = 0;
for (int i = 0; (i < VEC_SIZE); i++) {
sum += (thisRaw[r][i] * otherRaw[i][c]);
}
product.data[r][c] = sum;
}
}
delete[] thisRaw;
delete[] otherRaw;
return product;
}
#endif
// Default c'tor
Matrix::Matrix()
: Matrix(1,1) { }
Matrix::~Matrix() { }
// Rows/Cols c'tor
Matrix::Matrix(int rows, int cols)
: data(TwoDVec(rows, std::vector<int>(cols, 0))) { }
// Init func c'tor
Matrix::Matrix(int rows, int cols, std::function<int(int, int)> init)
: Matrix(rows, cols) {
for (int r = 0; (r < rows); r++) {
for (int c = 0; (c < cols); c++) {
data[r][c] = init(r,c);
}
}
}
// Copy c'tor
Matrix::Matrix(const Matrix& other) {
int rows = other.data.size();
int cols = rows > 0 ? other.data[0].size() : 0;
this->data.resize(rows, std::vector<int>(cols, 0));
for(int r = 0; (r < rows); r++) {
this->data[r] = other.data[r];
}
}
// Move c'tor
Matrix::Matrix(Matrix&& other) {
if (this == &other) return;
this->data.clear();
int rows = other.data.size();
for (int r = 0; (r < rows); r++) {
this->data[r] = std::move(other.data[r]);
}
}
Matrix&
Matrix::operator=(const Matrix& other) {
int rows = other.data.size();
this->data.resize(rows, std::vector<int>(0));
for (int r = 0; (r < rows); r++) {
this->data[r] = other.data[r];
}
return *this;
}
int
Matrix::sumDiag() const {
int rows = data.size();
int cols = rows > 0 ? data[0].size() : 0;
int dim = (rows < cols ? rows : cols);
int sum = 0;
for (int i = 0; (i < dim); i++) {
sum += data[i][i];
}
return sum;
}
int**
Matrix::getRawData() const {
int** rawData = new int*[data.size()];
for (int i = 0; (i < data.size()); i++) {
rawData[i] = const_cast<int*>(data[i].data());
}
return rawData;
}
【问题讨论】：
-
VEC_SIZE 有多大，您使用的是什么编译器？序列号是什么样的？
-
rowsThis 和 colsOther 有多大？ VEC_SIZE 是常数吗？你还记得释放使用“new int[VEC_SIZE]”创建的临时数组吗？
-
我已经进行了适当的修改。矩阵大小的编译时感知是否为这些 Cilk 操作提供了改进的优化机会？
-
大小可能有缓存效果。例如，如果 VEC_SIZE 很大，则临时数组总和可能会溢出 L1 缓存。
-
另一个问题：是否涉及C99变长数组？ IE。每个矩阵的最后一维是编译时常数吗？
标签： arrays caching vectorization matrix-multiplication cilk-plus