Matlab中特征向量间距离矩阵的并行mex程序

Posted 2020-09-30 clnchanpin

在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