在matlab中, 有n个向量(m维)的矩阵Mat(n, m)
要计算任两个向量间的距离, 即距离矩阵, 可使用下面的并行算法以加速:
#include <iostream>
#include <mex.h>
#include <matrix.h>
#include <thread>
using namespace std;
//预定义线程数
const int nThreads = 4;
//全局变量
int rows, cols, nrow, nw;
double *inVals, *outVals;
//线程执行体定义
void calc(int start, int end) {
double sum, tmp;
int no, i, j;
//计算指定区间
for(no = start; no < end; no++) {
//第i输入向量
i = outVals[no + nrow] - 1; //C索引下标
//第j输入向量
j = outVals[no + 2 * nrow] - 1; //C索引下标
//计算两输入向量间的距离
sum = 0;
for(int k = 0; k < cols; k++)
{
tmp = (inVals[i + k * rows] - inVals[j + k * rows]);
sum += (tmp * tmp);
}
outVals[no + 2 * nrow] = sum;
}
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
if (nrhs != 1) {
mexPrintf("Usage: adjmat(double_features[n_rows * m_cols_features])\n");
}
//指针指向输入数据
inVals = mxGetPr(prhs[0]);
//输入矩阵的行数和列数
rows = mxGetM(prhs[0]);
cols = mxGetN(prhs[0]);
//结果的行数nrow
//结果的列数nw=(i, j, distance)
nrow = (rows * rows - rows)/2, nw = 3;
//分配结果内存
nlhs = 1;
plhs[0] = mxCreateDoubleMatrix(nrow, nw, mxREAL);
outVals = mxGetPr(plhs[0]);
//在结果中分配i和j的组合
int curL = 0;
for(int i = 0; i < rows - 1; i++)
for(int j = i + 1; j < rows; j++) {
outVals[curL] = i + 1; //符合Matlab索引下标规范
outVals[curL + nrow] = j + 1; //符合Matlab索引下标规范
curL++;
}
//按线程数分配计算区间
int seg = nrow / nThreads;
//线程数组
thread threads[nThreads];
//分配每个线程的计算区间,避免冲突
for(int i = 0; i < nThreads; i++) {
if (i == nThreads - 1)
threads[i] = thread(calc, i * seg, nrow);
else
threads[i] = thread(calc, i * seg, (i + 1) * seg);
}
//等待所有线程结束
for (int i = 0; i < nThreads; i++){
threads[i].join();
}
}
编译: (注意:看上一篇博文,如何设置matlab支持C++ 11标准)
mex adjmat.cpp
Matlab中简单测试:
tic; x = rand(5000, 50);
adjmat(x);
toc
笔记本测试时间约:
0.57s