MATLAB+C+Openmp混合编程（4）-openmp基础

本文给出几种常用的并行思路，让大家快速理解并上手。
single、master、nowait等操作需要后续再进行展开，本节不进行分析。
1.parallel for
parallel for将会自动完成线程的运算任务分配，简单但是可操作性低，同时还有一系列的调度命令例如static，dynamic，guided，auto，runtime。对于一般算法来说，调度方式的不同对效率影响不大，多数情况不进行申明。

#include <omp.h> 
void main()
{
#pragma omp parallel for num_threads(6)
			for(int i = 0;i<12;i++)
		{
			printf("OpenMP线程编号为:%d\n", omp_get_thread_num());
		}
		system("pause");
}

如果想改成matlab的mex函数则代码改为

#include <omp.h> 
#include<mex.h>
void mexFunction(int nlhs, mxArray* plhs[],
    int nrhs, const mxArray* prhs[])
{
#pragma omp parallel for num_threads(6)
			for(int i = 0;i<12;i++)
		{
			printf("OpenMP线程编号为:%d\n", omp_get_thread_num());
		}
		system("pause");
}

这个代码编译后运行MATLAB会直接报错，原因就是并行区内不支持printf和mexPrintf函数
本质上说mex项目的编写只是把原cpp项目的main函数改为
void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[])
然后指定输入和输出数组的指针即可
案例代码
将matlab输入的矩阵A分别乘以2和3，输出两个新的矩阵B和C

// mex代码
#include "mex.h"
#include<stdio.h>
#include <omp.h> 
#include<iostream>
#include<ctime>
using namespace std;
// 使用MEX必须包含的头文件
void hh_kk_computing_omp(double* zbd,double*hh, double *kk,int M_zbd, int N_zbd);
void mexFunction(int nlhs, mxArray* plhs[],
    int nrhs, const mxArray* prhs[]) {
    clock_t st1, en1;
    st1 = clock();
    if (nrhs != 1)
        mexErrMsgTxt("Wrong number of input arguments.\n");
    // 检查输入变量数量是否正确，否则报错
#define zbd_IN  prhs[0]
#define hh_OUT plhs[0]
#define kk_OUT plhs[1]
    int M_zbd = mxGetM(zbd_IN);
    int N_zbd = mxGetN(zbd_IN);

    hh_OUT = mxCreateDoubleMatrix(M_zbd, N_zbd, mxREAL);
    kk_OUT = mxCreateDoubleMatrix(M_zbd, N_zbd, mxREAL);
    double* A = mxGetPr(zbd_IN);
    double* B = mxGetPr(hh_OUT);
    double* C = mxGetPr(kk_OUT);
    hh_kk_computing_omp(A, B, C,M_zbd,N_zbd);

}
void hh_kk_computing_omp(double* zbd,double *hh, double *kk,int M_zbd, int N_zbd)
{
#pragma omp parallel for num_threads(6) schedule(guided)
       for (int i = 0; i < M_zbd*N_zbd; i++)
        {
               hh[i] = 2*zbd[i];
               kk[i] = 3*zbd[i];

        }

}

编译命令见MATLAB+C+Openmp混合编程（补充）-mex编译命令
注意修改路径
MATLAB测试代码如下

%MATLAB测试代码
a=rand(2,2)
[b,c]=main(a)

2.parallel
parallel与parallel for不同，parallel的可操作性强，同时是很多算法采用的并行方式，parallel仅申明一个并行区域，我们可以手动分配并行区域内线程的运算任务

#include <omp.h> 
#include<mex.h>
void main()
{
int num_threads=6;
int a[10000]={0}
#pragma omp parallel  num_threads(num_threads)
{
    int tid=omp_get_thread_num();
    int istart=10000*tid/num_threads;
    int iend=10000*(tid+1)/num_threads;
    for (int k=istart;k<iend;++k)
    {
      a[k]+=1;
}
}
}

如上代码实现了a[10000]数组的自增，tid用于获取线程id，假定在并行区申请了6个线程，那么线程的tid分别为0,1,2,3,4,5
线程tid为0的线程负责a[0]-a[1667]的运算任务，剩余线程以此类推
更进一步，我们可以将6个线程的运算任务进一步复杂化，仅前三个线程处理a数组，后三个线程可以处理其他运算任务

#include <omp.h> 
void main()
{
int num_threads=6;
int a[10000]={0}
#pragma omp parallel  num_threads(num_threads)
{
    int tid=omp_get_thread_num();
    int istart=10000*(tid/3);
    int iend=10000*(3+1)/num_threads;
    if (tid<=2)
    {
    for (int k=istart;k<iend;++k)
    {
      a[k]+=1;
}
}
}
}

3.嵌套并行
嵌套并行的实现原理类似于二叉树，一开始仅有一个线程0（理解为根节点），然后根节点分裂出左右子节点(线程0申请线程1,2），如此反复，线程个数不断增加，嵌套并行在实际并行算法中的运用较为复杂，一般涉及到递归算法或者树形结构建立，这里仅给出简单的代码，具体细节请参考《多核异构并行计算OpenMP4.5 C/C++篇》
（vs不支持openmp的嵌套并行!推荐采用code::blocks)

#include <omp.h>
#include<stdio.h>
#include<ctime>
int main()
{
    omp_set_nested(1);
    omp_set_dynamic(0);
#pragma omp parallel num_threads(3)
    {
        if (omp_get_thread_num() == 0)
        {
            omp_set_num_threads(2);       /* line A */

        }
        else
        {
            omp_set_num_threads(5);
            /* line B */
        }
//printf("outer  thread id is %d:threads outer are %d inside threads  is %d\n", omp_get_thread_num(),omp_get_num_threads(),omp_get_max_threads());
//printf("  thread id is %d Inner threads are: %d\n",omp_get_thread_num(),omp_get_max_threads());
//printf("  thread level is %d in this level are %d threads\n",omp_get_active_level(),omp_get_team_size(omp_get_active_level()));

    #pragma omp parallel
    {
        #pragma omp master
        {
            printf("master thread is now in level %d thread id is %d it's team size is %d\n",omp_get_active_level(),omp_get_thread_num(),omp_get_team_size(omp_get_active_level()));
        }
        printf("thread is now in level %d thread id is %d it's team size is %d\n",omp_get_active_level(),omp_get_thread_num(),omp_get_team_size(omp_get_active_level()));

    }


    }
    return 0;
}

练习题
编写mex函数，函数需要实现如下功能:
统计MATLAB矩阵A中小于等于0的数值个数，函数支持openmp并行，当用户输入输入1时为并行统计，输入0时为串行统计
MATLAB测试代码如下

M=500;
N=500;
A=(randi([-50000,58020],M,N));
matlab_result=0;
for i=1:1:M*N
 if(A(i)<=0)
   matlab_result=matlab_result+1;
 end
end
mex_result_parallel=mex_function(A,1);
mex_result_serial=mex_function(A,0);
if(mex_result_parallel==mex_result_serial&&mex_result_parallel==matlab_result)
 fprintf("correct")
end

下一节将给出mex代码