Skip to content
This repository was archived by the owner on Mar 24, 2026. It is now read-only.

Commit f6ddc50

Browse files
bibikaranton-malakhov
bibikar
authored and
anton-malakhov
committed
Add native benchmarks for linear algebra (#7)
Add native benchmarks for linear algebra
1 parent d84b089 commit f6ddc50

17 files changed

Lines changed: 812 additions & 0 deletions

File tree

‎numpy/linalg/Makefile‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
# Copyright (C) 2019 Intel Corporation
2+
#
3+
# SPDX-License-Identifier: MIT
4+
5+
CXX = icc
6+
CXXFLAGS = -O3 -g -xCORE-AVX2 -axCOMMON-AVX512 -qopenmp
7+
override CXXFLAGS += -std=c++11
8+
LDFLAGS = -lmkl_rt -qopenmp
9+
10+
TARGET = linalg
11+
BENCHES = cholesky det dot inv lu qr svd
12+
SOURCES = $(addsuffix .cc,$(BENCHES)) linalg.cc
13+
14+
$(TARGET): $(SOURCES:.cc=.o)
15+
$(CXX) $^ -o $@ $(LDFLAGS)
16+
17+
clean:
18+
rm -f $(TARGET) *.o
19+

‎numpy/linalg/bench.h‎

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
/*
2+
* Copyright (C) 2016, 2019 Intel Corporation
3+
*
4+
* SPDX-License-Identifier: MIT
5+
*/
6+
7+
#pragma once
8+
#include <algorithm>
9+
10+
using namespace std;
11+
12+
#include "assert.h"
13+
#include "stdlib.h"
14+
15+
#if defined(__INTEL_COMPILER)
16+
17+
#include "mkl.h"
18+
19+
class Random {
20+
private:
21+
enum { SEED = 77777 };
22+
static double const d_zero = 0.0, d_one = 1.0;
23+
VSLStreamStatePtr stream;
24+
25+
public:
26+
Random() {
27+
int err = vslNewStream(&stream, VSL_BRNG_MT19937, SEED);
28+
assert(err == VSL_STATUS_OK);
29+
}
30+
~Random() {
31+
int err = vslDeleteStream(&stream);
32+
assert(err == VSL_STATUS_OK);
33+
}
34+
void init_mat(double *mat, int size) {
35+
int err = vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF, stream, size, mat,
36+
d_zero, d_one);
37+
assert(err == VSL_STATUS_OK);
38+
}
39+
};
40+
41+
#else
42+
43+
#include "cblas.h"
44+
#include "lapacke.h"
45+
46+
class Random {
47+
public:
48+
void init_mat(double *mat, int size) {
49+
}
50+
};
51+
52+
static void *mkl_malloc(int size, int align) {
53+
return malloc(size);
54+
}
55+
56+
static void mkl_free(void *p) {
57+
free(p);
58+
}
59+
#endif
60+
61+
class Bench {
62+
private:
63+
Random random;
64+
65+
public:
66+
double *make_random_mat(int size) {
67+
double *mat = make_mat(size);
68+
random.init_mat(mat, size);
69+
return mat;
70+
}
71+
72+
double *make_mat(int mat_size) {
73+
double *mat = (double *) mkl_malloc(mat_size * sizeof(double), 64);
74+
assert(mat);
75+
return mat;
76+
}
77+
78+
virtual void make_args(int size) = 0;
79+
virtual void copy_args() = 0;
80+
virtual void compute() = 0;
81+
};

‎numpy/linalg/cholesky.cc‎

Lines changed: 53 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,53 @@
1+
/*
2+
* Copyright (C) 2019 Intel Corporation
3+
*
4+
* SPDX-License-Identifier: MIT
5+
*/
6+
7+
#include "cholesky.h"
8+
#include <cstring>
9+
10+
Cholesky::Cholesky() {
11+
x_mat = r_mat = 0;
12+
}
13+
14+
void Cholesky::make_args(int size) {
15+
n = lda = size;
16+
17+
mat_size = n * n;
18+
int r_size = mat_size;
19+
20+
// input matrix
21+
x_mat = make_random_mat(mat_size);
22+
23+
// matrix for result
24+
r_mat = make_mat(r_size);
25+
memset(r_mat, 0, r_size * sizeof(*r_mat));
26+
// Set r_mat to identity matrix as in python bench
27+
for (int i = 0; i < n; i++) {
28+
r_mat[i * n + i] = 1;
29+
}
30+
cblas_dsyrk(CblasColMajor, CblasUpper, CblasNoTrans, n, n, 1.0, x_mat, lda,
31+
n, r_mat, lda);
32+
33+
// we now have r_mat = x_mat * x_mat' + n * np.eye(n)
34+
// copy back into x_mat
35+
memcpy(x_mat, r_mat, mat_size * sizeof(*x_mat));
36+
}
37+
38+
void Cholesky::copy_args() {
39+
memcpy(r_mat, x_mat, mat_size * sizeof(*x_mat));
40+
}
41+
42+
void Cholesky::compute() {
43+
// compute cholesky decomposition
44+
int info = LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'U', n, r_mat, lda);
45+
assert(info == 0);
46+
}
47+
48+
Cholesky::~Cholesky() {
49+
if (r_mat)
50+
mkl_free(r_mat);
51+
if (x_mat)
52+
mkl_free(x_mat);
53+
}

‎numpy/linalg/cholesky.h‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
/*
2+
* Copyright (C) 2019 Intel Corporation
3+
*
4+
* SPDX-License-Identifier: MIT
5+
*/
6+
7+
#include "bench.h"
8+
9+
class Cholesky : public Bench {
10+
public:
11+
Cholesky();
12+
~Cholesky();
13+
void make_args(int size);
14+
void copy_args();
15+
void compute();
16+
17+
private:
18+
double *x_mat, *r_mat;
19+
int n, lda, mat_size;
20+
};

‎numpy/linalg/det.cc‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
/*
2+
* Copyright (C) 2016, 2019 Intel Corporation
3+
*
4+
* SPDX-License-Identifier: MIT
5+
*/
6+
7+
#include "det.h"
8+
#include <cstring>
9+
10+
Det::Det() {
11+
r_mat = x_mat = 0;
12+
ipiv = 0;
13+
}
14+
15+
void Det::make_args(int size) {
16+
n = size;
17+
m = size;
18+
mn_min = min(m, n);
19+
lda = size;
20+
int mat_size = m * n;
21+
assert(m == n);
22+
23+
// input matrix
24+
x_mat = make_random_mat(mat_size);
25+
26+
// list of pivots
27+
ipiv = (int *) mkl_malloc(mn_min * sizeof(int), 64);
28+
assert(ipiv);
29+
30+
// matrix for result
31+
r_mat = make_random_mat(mat_size);
32+
33+
copy_args();
34+
}
35+
36+
void Det::copy_args() {
37+
memcpy(r_mat, x_mat, mat_size * sizeof(*r_mat));
38+
}
39+
40+
void Det::compute() {
41+
// compute pivoted lu decomposition
42+
int info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, m, n, r_mat, lda, ipiv);
43+
assert(info == 0);
44+
45+
double t = 1.0;
46+
int i, j;
47+
for (i = 0, j = 0; i < mn_min; i++, j += lda + 1) {
48+
t *= (ipiv[i] == i) ? r_mat[j] : -r_mat[j];
49+
}
50+
result = t;
51+
}
52+
53+
Det::~Det() {
54+
if (r_mat)
55+
mkl_free(r_mat);
56+
if (ipiv)
57+
mkl_free(ipiv);
58+
if (x_mat)
59+
mkl_free(x_mat);
60+
}

‎numpy/linalg/det.h‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,22 @@
1+
/*
2+
* Copyright (C) 2016, 2019 Intel Corporation
3+
*
4+
* SPDX-License-Identifier: MIT
5+
*/
6+
7+
#include "bench.h"
8+
9+
class Det : public Bench {
10+
public:
11+
Det();
12+
~Det();
13+
void make_args(int size);
14+
void copy_args();
15+
void compute();
16+
17+
private:
18+
double *x_mat, *r_mat;
19+
int *ipiv;
20+
int m, n, lda, mn_min, mat_size;
21+
double result;
22+
};

‎numpy/linalg/dot.cc‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,47 @@
1+
/*
2+
* Copyright (C) 2016, 2019 Intel Corporation
3+
*
4+
* SPDX-License-Identifier: MIT
5+
*/
6+
7+
#include "dot.h"
8+
#include <cstring>
9+
10+
Dot::Dot() {
11+
a_mat = b_mat = c_mat = r_mat = 0;
12+
}
13+
14+
Dot::~Dot() {
15+
if (a_mat)
16+
mkl_free(a_mat);
17+
if (b_mat)
18+
mkl_free(b_mat);
19+
if (c_mat)
20+
mkl_free(c_mat);
21+
if (r_mat)
22+
mkl_free(r_mat);
23+
}
24+
25+
void Dot::make_args(int size) {
26+
m = n = k = size;
27+
28+
a_mat = make_random_mat(m * k);
29+
b_mat = make_random_mat(k * n);
30+
c_mat = make_random_mat(m * n);
31+
32+
r_mat = make_mat(m * n);
33+
34+
copy_args();
35+
}
36+
37+
void Dot::copy_args() {
38+
memcpy(r_mat, c_mat, m * n * sizeof(*r_mat));
39+
}
40+
41+
void Dot::compute() {
42+
double alpha = 1.0;
43+
double beta = 0.0;
44+
45+
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, n, n, k, alpha,
46+
a_mat, k, b_mat, n, beta, r_mat, n);
47+
}

‎numpy/linalg/dot.h‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
/*
2+
* Copyright (C) 2016, 2019 Intel Corporation
3+
*
4+
* SPDX-License-Identifier: MIT
5+
*/
6+
7+
#include "bench.h"
8+
9+
class Dot : public Bench {
10+
public:
11+
Dot();
12+
~Dot();
13+
void make_args(int size);
14+
void copy_args();
15+
void compute();
16+
17+
private:
18+
double *a_mat, *b_mat, *c_mat, *r_mat;
19+
int m, n, k;
20+
};

‎numpy/linalg/inv.cc‎

Lines changed: 57 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,57 @@
1+
/*
2+
* Copyright (C) 2016, 2019 Intel Corporation
3+
*
4+
* SPDX-License-Identifier: MIT
5+
*/
6+
7+
#include "inv.h"
8+
#include <cstring>
9+
10+
Inv::Inv() {
11+
x_mat = 0;
12+
r_mat = 0;
13+
ipiv = 0;
14+
}
15+
16+
Inv::~Inv() {
17+
if (r_mat)
18+
mkl_free(r_mat);
19+
if (ipiv)
20+
mkl_free(ipiv);
21+
if (x_mat)
22+
mkl_free(x_mat);
23+
}
24+
25+
void Inv::make_args(int size) {
26+
n = size;
27+
m = size;
28+
lda = size;
29+
mat_size = m * n;
30+
int mn_min = min(m, n);
31+
32+
assert(m == n);
33+
34+
// input matrix
35+
x_mat = make_random_mat(mat_size);
36+
37+
// list of pivots
38+
ipiv = (int *) mkl_malloc(mn_min * sizeof(int), 64);
39+
assert(ipiv);
40+
41+
// matrix for result
42+
r_mat = make_mat(mat_size);
43+
copy_args();
44+
}
45+
46+
void Inv::copy_args() {
47+
memcpy(r_mat, x_mat, mat_size * sizeof(*r_mat));
48+
}
49+
50+
void Inv::compute() {
51+
// compute pivoted lu decomposition
52+
int info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, m, n, r_mat, lda, ipiv);
53+
assert(info == 0);
54+
55+
info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, r_mat, lda, ipiv);
56+
assert(info == 0);
57+
}

0 commit comments

Comments
 (0)