利用VS2019+MKL加速矩阵运算。

软件版本

Intel OneAP 2022.0.2
Visual Studio 2019 Community

工程配置

• 新建空工程

  
  
  image.png

• 打开OneAPI选项

  
  
  image.png

• 加入代码（官方示例）

#include <stdio.h>
#include <time.h>
#include <stdlib.h>
#include "mkl.h"

void print_arr(int N, char* name, double* array);
void init_arr(int N, double* a);
void Dgemm_multiply(double* a, double* b, double* c, int N);

int main(int argc, char* argv[])
{
    clock_t start, stop;
    int i, j;
    int N;
    double* a;
    double* b;
    double* c;
    if (argc < 2)
    {
        printf("Enter matrix size N=");
        //please enter small number first to ensure that the 
        //multiplication is correct! and then you may enter 
        //a "reasonably" large number say like 500 or even 1000
        scanf("%d", &N);
    }
    else
    {
        N = atoi(argv[1]);
    }

    a = (double*)malloc(sizeof(double) * N * N);
    b = (double*)malloc(sizeof(double) * N * N);
    c = (double*)malloc(sizeof(double) * N * N);

    init_arr(N, a);
    init_arr(N, b);

    //DGEMM Multiply
    //reallocate to force cash to be flushed
    a = (double*)malloc(sizeof(double) * N * N);
    b = (double*)malloc(sizeof(double) * N * N);
    c = (double*)malloc(sizeof(double) * N * N);
    init_arr(N, a);
    init_arr(N, b);

    start = clock();
    //for(i=0;i<1000;i++)
    Dgemm_multiply(a, b, c, N);
    stop = clock();

    printf("Dgemm_multiply(). Elapsed time = %g seconds\n",
        ((double)(stop - start)) / CLOCKS_PER_SEC);
    //print simple test case of data to be sure multiplication is correct
    if (N < 7) {
        print_arr(N, "a", a);
        print_arr(N, "b", b);
        print_arr(N, "c", c);
    }

    free(a);
    free(b);
    free(c);

    return 0;
}


//DGEMM way. The PREFERED way, especially for large matrices
void Dgemm_multiply(double* a, double* b, double* c, int N)
{

    double alpha = 1.0, beta = 0.;
    int incx = 1;
    int incy = N;
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, N, N, N, alpha, b, N, a, N, beta, c, N);
}

//initialize array with random data
void init_arr(int N, double* a)
{
    int i, j;
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            a[i * N + j] = (i + j + 1) % 10; //keep all entries less than 10. pleasing to the eye!
        }
    }
}

//print array to std out
void print_arr(int N, char* name, double* array)
{
    int i, j;
    printf("\n%s\n", name);
    for (i = 0; i < N; i++) {
        for (j = 0; j < N; j++) {
            printf("%g\t", array[N * i + j]);
        }
        printf("\n");
    }
}

运行和调试

• 偶然出现的奇怪问题，在设置OneAPI方式为Sequential的时候，程序一直可以运行，但是OneAPI设置为Parallel后，输入的N只能是小于15（多次实测），不然就会报错！

  
  
  Parallel N=14
  
  
  Parallel N=15
  
  
  错误

• 解决方法，在debug出来的exe文件夹下放入libiomp5md.dll

点亮“在看”
欢迎关注、收藏、转载哈!
最有意思的GNSS NEWS…不定期更新…