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

Commit fe6d24d

Browse files
bibikarbibikar
committed
Use Fortran functions for det, inv, cholesky (#12)
Run clang-format
1 parent 683767e commit fe6d24d

9 files changed

Lines changed: 157 additions & 164 deletions

File tree

‎numpy/linalg/cholesky.cc‎

Lines changed: 21 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -5,24 +5,26 @@
55
*/
66

77
#include "cholesky.h"
8-
#include <iostream>
98
#include <cstring>
9+
#include <iostream>
1010

1111
static const double x_mat_test[] = {
12-
5.551063927745538, 0.034194385271978, -0.276508795460738,
13-
0.034194385271978, 4.704686460853461, 0.087572555571367,
14-
-0.276508795460738, 0.087572555571367, 6.07658590927362
15-
};
16-
17-
static const double r_mat_test[] = {
18-
2.356069593145656, 0. , 0. ,
19-
0.014513317166631, 2.16898036516661 , 0. ,
20-
-0.117360198639788, 0.041160281019929, 2.461933858639425
21-
};
12+
5.551063927745538, 0.034194385271978, -0.276508795460738,
13+
0.034194385271978, 4.704686460853461, 0.087572555571367,
14+
-0.276508795460738, 0.087572555571367, 6.07658590927362};
15+
16+
static const double r_mat_test[] = {2.356069593145656,
17+
0.,
18+
0.,
19+
0.014513317166631,
20+
2.16898036516661,
21+
0.,
22+
-0.117360198639788,
23+
0.041160281019929,
24+
2.461933858639425};
2225

2326
static const int test_size = 3;
2427

25-
2628
Cholesky::Cholesky() {
2729
x_mat = r_mat = 0;
2830
}
@@ -43,8 +45,8 @@ void Cholesky::make_args(int size) {
4345
for (int i = 0; i < n; i++) {
4446
r_mat[i * n + i] = 1;
4547
}
46-
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0, x_mat,
47-
n, x_mat, n, n, r_mat, n);
48+
cblas_dgemm(CblasColMajor, CblasNoTrans, CblasTrans, n, n, n, 1.0, x_mat, n,
49+
x_mat, n, n, r_mat, n);
4850

4951
// we now have r_mat = x_mat * x_mat' + n * np.eye(n)
5052
// copy back into x_mat
@@ -57,11 +59,13 @@ void Cholesky::copy_args() {
5759

5860
void Cholesky::compute() {
5961
// compute cholesky decomposition
60-
int info = LAPACKE_dpotrf(LAPACK_COL_MAJOR, 'U', n, r_mat, lda);
62+
int info;
63+
const char uplo = 'U';
64+
dpotrf(&uplo, &n, r_mat, &lda, &info);
6165
assert(info == 0);
6266

6367
// we only want an upper triangular matrix
64-
for (int i = 0; i < n-1; i++) {
68+
for (int i = 0; i < n - 1; i++) {
6569
memset(&r_mat[i * n + i + 1], 0, (n - i - 1) * sizeof(*r_mat));
6670
}
6771
}
@@ -76,10 +80,9 @@ bool Cholesky::test() {
7680
return mat_equal(r_mat, r_mat_test, mat_size);
7781
}
7882

79-
8083
void Cholesky::print_args() {
8184
std::cout << "Cholesky decomposition, A = LL*, of a "
82-
<< "Hermitian positive-definite matrix A." << std::endl;
85+
<< "Hermitian positive-definite matrix A." << std::endl;
8386
std::cout << "A = " << std::endl;
8487
print_mat('c', x_mat, n, n);
8588
}

‎numpy/linalg/det.cc‎

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -5,14 +5,13 @@
55
*/
66

77
#include "det.h"
8-
#include <iostream>
98
#include <cstring>
9+
#include <iostream>
1010

1111
static const double x_mat_test[] = {
12-
0.470442000675409, -0.291482508170914, -0.44183986349643 ,
13-
-0.176333746005435, 0.007410393215614, -0.739195206041762,
14-
0.481736547564898, 0.805743972141035, -0.468344563609981
15-
};
12+
0.470442000675409, -0.291482508170914, -0.44183986349643,
13+
-0.176333746005435, 0.007410393215614, -0.739195206041762,
14+
0.481736547564898, 0.805743972141035, -0.468344563609981};
1615

1716
static const double result_test = 0.4707855751774963;
1817
static const int test_size = 3;
@@ -49,13 +48,14 @@ void Det::copy_args() {
4948

5049
void Det::compute() {
5150
// compute pivoted lu decomposition
52-
int info = LAPACKE_dgetrf(LAPACK_COL_MAJOR, m, n, r_mat, lda, ipiv);
51+
int info;
52+
dgetrf(&n, &n, r_mat, &lda, ipiv, &info);
5353
assert(info == 0);
5454

5555
double t = 1.0;
5656
int i, j;
5757
for (i = 0, j = 0; i < mn_min; i++, j += lda + 1) {
58-
t *= (ipiv[i] == i+1) ? r_mat[j] : -r_mat[j];
58+
t *= (ipiv[i] == i + 1) ? r_mat[j] : -r_mat[j];
5959
}
6060
result = t;
6161
}
@@ -71,7 +71,8 @@ bool Det::test() {
7171
}
7272

7373
void Det::print_args() {
74-
std::cout << "Determinant of " << n << "x" << n << " matrix A." << std::endl;
74+
std::cout << "Determinant of " << n << "x" << n << " matrix A."
75+
<< std::endl;
7576
std::cout << "A =" << std::endl;
7677
print_mat('c', x_mat, n, n);
7778
}

‎numpy/linalg/dot.cc‎

Lines changed: 10 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -5,26 +5,23 @@
55
*/
66

77
#include "dot.h"
8-
#include <iostream>
98
#include <cstring>
9+
#include <iostream>
1010

1111
static const double a_mat_test[] = {
12-
0.470442000675409, -0.176333746005435, 0.481736547564898,
13-
-0.291482508170914, 0.007410393215614, 0.805743972141035,
14-
-0.44183986349643 , -0.739195206041762, -0.468344563609981
15-
};
12+
0.470442000675409, -0.176333746005435, 0.481736547564898,
13+
-0.291482508170914, 0.007410393215614, 0.805743972141035,
14+
-0.44183986349643, -0.739195206041762, -0.468344563609981};
1615

1716
static const double b_mat_test[] = {
18-
-0.279551412836935, -1.866235595807669, 0.949267732307811,
19-
0.393910888693485, 0.357832357041521, 1.430099195743549,
20-
-0.202579028296422, -1.225349132812327, 0.535350173863021
21-
};
17+
-0.279551412836935, -1.866235595807669, 0.949267732307811,
18+
0.393910888693485, 0.357832357041521, 1.430099195743549,
19+
-0.202579028296422, -1.225349132812327, 0.535350173863021};
2220

2321
static const double r_mat_test[] = {
24-
-0.298562230242867, -1.531348988185161, 0.452298407313714,
25-
-0.078823449378468, -0.44069096675697 , 0.165257833413781,
26-
-0.072783295837747, 1.133954922888981, -1.727275138478663
27-
};
22+
-0.298562230242867, -1.531348988185161, 0.452298407313714,
23+
-0.078823449378468, -0.44069096675697, 0.165257833413781,
24+
-0.072783295837747, 1.133954922888981, -1.727275138478663};
2825

2926
static const int test_size = 3;
3027

@@ -82,7 +79,6 @@ bool Dot::test() {
8279
return mat_equal(r_mat, r_mat_test, m * n);
8380
}
8481

85-
8682
void Dot::print_args() {
8783
std::cout << "Matrix multiplication A * B." << std::endl;
8884
std::cout << "A =" << std::endl;
@@ -95,4 +91,3 @@ void Dot::print_result() {
9591
std::cout << "A * B =" << std::endl;
9692
print_mat('r', r_mat, m, n);
9793
}
98-

‎numpy/linalg/eig.cc‎

Lines changed: 24 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -9,39 +9,32 @@
99
#include <iostream>
1010

1111
static const double a_mat_test[] = {
12-
0.470442000675409, -0.291482508170914, -0.44183986349643 ,
13-
-0.176333746005435, 0.007410393215614, -0.739195206041762,
14-
0.481736547564898, 0.805743972141035, -0.468344563609981
15-
};
12+
0.470442000675409, -0.291482508170914, -0.44183986349643,
13+
-0.176333746005435, 0.007410393215614, -0.739195206041762,
14+
0.481736547564898, 0.805743972141035, -0.468344563609981};
1615

1716
static const std::complex<double> w_vec_complex_test[] = {
18-
{ 0.558344640162537, 0. },
19-
{-0.274418404940748, 0.876285061400947},
20-
{-0.274418404940748, -0.876285061400947}
21-
};
17+
{0.558344640162537, 0.},
18+
{-0.274418404940748, 0.876285061400947},
19+
{-0.274418404940748, -0.876285061400947}};
2220

2321
static const std::complex<double> vr_mat_complex_test[] = {
24-
{-0.870718618937641, 0. },
25-
{ 0.491334396303628, 0. },
26-
{ 0.020966583991578, 0. },
27-
{ 0.12292498568887 , 0.296216896839797},
28-
{ 0.107429842088307, 0.640393071981176},
29-
{-0.689565472096294, 0. },
30-
{ 0.12292498568887 , -0.296216896839797},
31-
{ 0.107429842088307, -0.640393071981176},
32-
{-0.689565472096294, -0. }
33-
};
22+
{-0.870718618937641, 0.},
23+
{0.491334396303628, 0.},
24+
{0.020966583991578, 0.},
25+
{0.12292498568887, 0.296216896839797},
26+
{0.107429842088307, 0.640393071981176},
27+
{-0.689565472096294, 0.},
28+
{0.12292498568887, -0.296216896839797},
29+
{0.107429842088307, -0.640393071981176},
30+
{-0.689565472096294, -0.}};
3431

3532
static const int test_size = 3;
3633

3734
// We set these to zero, because the expectation is complex
3835
// eigenvalues and eigenvectors for this test input
3936
static const double wr_vec_test[] = {0., 0., 0.};
40-
static const double vr_mat_test[] = {
41-
0., 0., 0.,
42-
0., 0., 0.,
43-
0., 0., 0.
44-
};
37+
static const double vr_mat_test[] = {0., 0., 0., 0., 0., 0., 0., 0., 0.};
4538

4639
Eig::Eig() {
4740
a_mat = r_mat = vl_mat = vr_mat = wr_vec = wi_vec = 0;
@@ -68,8 +61,8 @@ void Eig::make_args(int size) {
6861
// complex eigenvalues and eigenvectors
6962
w_vec_complex =
7063
(std::complex<double> *) mkl_malloc(n * sizeof(*w_vec_complex), 64);
71-
vr_mat_complex =
72-
(std::complex<double> *) mkl_malloc(mat_size * sizeof(*vr_mat_complex), 64);
64+
vr_mat_complex = (std::complex<double> *) mkl_malloc(
65+
mat_size * sizeof(*vr_mat_complex), 64);
7366
}
7467

7568
void Eig::copy_args() {
@@ -129,17 +122,16 @@ bool Eig::test() {
129122
compute();
130123

131124
if (only_real)
132-
return mat_equal(wr_vec, wr_vec_test, n)
133-
&& mat_equal(vr_mat, vr_mat_test, mat_size);
125+
return mat_equal(wr_vec, wr_vec_test, n) &&
126+
mat_equal(vr_mat, vr_mat_test, mat_size);
134127
else
135-
return mat_equal(w_vec_complex, w_vec_complex_test, n)
136-
&& mat_equal(vr_mat_complex, vr_mat_complex_test, mat_size);
128+
return mat_equal(w_vec_complex, w_vec_complex_test, n) &&
129+
mat_equal(vr_mat_complex, vr_mat_complex_test, mat_size);
137130
}
138131

139-
140132
void Eig::print_args() {
141-
std::cout << "Eigenvalues and eigenvectors of " <<
142-
n << "*" << n << " matrix A." << std::endl;
133+
std::cout << "Eigenvalues and eigenvectors of " << n << "*" << n
134+
<< " matrix A." << std::endl;
143135
std::cout << "A =" << std::endl;
144136
print_mat('c', a_mat, n, n);
145137
}

‎numpy/linalg/inv.cc‎

Lines changed: 33 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -5,25 +5,25 @@
55
*/
66

77
#include "inv.h"
8-
#include <iostream>
98
#include <cstring>
9+
#include <iostream>
1010

1111
static const double x_mat_test[] = {
12-
0.470442000675409, -0.291482508170914, -0.44183986349643 ,
13-
-0.176333746005435, 0.007410393215614, -0.739195206041762,
14-
0.481736547564898, 0.805743972141035, -0.468344563609981
15-
};
12+
0.470442000675409, -0.176333746005435, 0.481736547564898,
13+
-0.291482508170914, 0.007410393215614, 0.805743972141035,
14+
-0.44183986349643, -0.739195206041762, -0.468344563609981};
1615

1716
static const double r_mat_test[] = {
18-
1.257751926455496, -1.046174905778333, 0.46462060722515 ,
19-
-0.931809131348847, -0.01588524263938 , 0.904147816602771,
20-
-0.309375898184227, -1.103418649248418, -0.101770412852239
21-
};
17+
1.257751926455496, -0.931809131348847, -0.309375898184227,
18+
-1.046174905778332, -0.01588524263938, -1.103418649248418,
19+
0.46462060722515, 0.904147816602771, -0.101770412852238};
2220
static const int test_size = 3;
2321

2422
Inv::Inv() {
2523
x_mat = 0;
24+
x_mat_init = 0;
2625
r_mat = 0;
26+
identity = 0;
2727
ipiv = 0;
2828
}
2929

@@ -34,6 +34,10 @@ void Inv::clean_args() {
3434
mkl_free(ipiv);
3535
if (x_mat)
3636
mkl_free(x_mat);
37+
if (x_mat_init)
38+
mkl_free(x_mat_init);
39+
if (identity)
40+
mkl_free(identity);
3741
}
3842

3943
Inv::~Inv() {
@@ -42,55 +46,58 @@ Inv::~Inv() {
4246

4347
void Inv::make_args(int size) {
4448
n = size;
45-
m = size;
4649
lda = size;
47-
mat_size = m * n;
48-
int mn_min = min(m, n);
49-
50-
assert(m == n);
50+
mat_size = n * n;
5151

5252
// input matrix
53-
x_mat = make_random_mat(mat_size);
53+
x_mat_init = make_random_mat(mat_size);
54+
x_mat = make_mat(mat_size);
5455

5556
// list of pivots
56-
ipiv = (int *) mkl_malloc(mn_min * sizeof(int), 64);
57+
ipiv = (int *) mkl_malloc(n * sizeof(int), 64);
5758
assert(ipiv);
5859

5960
// matrix for result
6061
r_mat = make_mat(mat_size);
62+
63+
// identity matrix
64+
identity = make_mat(mat_size);
65+
memset(identity, 0, mat_size * sizeof(*identity));
66+
for (int i = 0; i < mat_size; i += n + 1)
67+
identity[i] = 1;
68+
6169
copy_args();
6270
}
6371

6472
void Inv::copy_args() {
65-
memcpy(r_mat, x_mat, mat_size * sizeof(*r_mat));
73+
memcpy(x_mat, x_mat_init, mat_size * sizeof(*x_mat));
74+
memcpy(r_mat, identity, mat_size * sizeof(*r_mat));
6675
}
6776

6877
void Inv::compute() {
69-
// compute pivoted lu decomposition
70-
int info = LAPACKE_dgetrf(LAPACK_ROW_MAJOR, m, n, r_mat, lda, ipiv);
71-
assert(info == 0);
72-
73-
info = LAPACKE_dgetri(LAPACK_ROW_MAJOR, n, r_mat, lda, ipiv);
78+
// Solve the equation X * X**-1 = I for X**-1.
79+
int info;
80+
dgesv(&n, &n, x_mat, &n, ipiv, r_mat, &n, &info);
7481
assert(info == 0);
7582
}
7683

7784
bool Inv::test() {
7885
clean_args();
7986
make_args(test_size);
80-
memcpy(x_mat, x_mat_test, mat_size * sizeof(*x_mat));
8187
copy_args();
88+
memcpy(x_mat, x_mat_test, mat_size * sizeof(*x_mat));
8289
compute();
8390

8491
return mat_equal(r_mat, r_mat_test, mat_size);
8592
}
8693

8794
void Inv::print_args() {
88-
std::cout << "Inverse of " << m << "*" << n << " matrix A." << std::endl;
95+
std::cout << "Inverse of " << n << "*" << n << " matrix A." << std::endl;
8996
std::cout << "A =" << std::endl;
90-
print_mat('r', x_mat, m, n);
97+
print_mat('c', x_mat, n, n);
9198
}
9299

93100
void Inv::print_result() {
94101
std::cout << "A**-1 =" << std::endl;
95-
print_mat('r', r_mat, m, n);
102+
print_mat('c', r_mat, n, n);
96103
}

‎numpy/linalg/inv.h‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,7 +19,7 @@ class Inv : public Bench {
1919
void compute();
2020

2121
private:
22-
double *x_mat, *r_mat;
22+
double *x_mat, *r_mat, *x_mat_init, *identity;
2323
int *ipiv;
24-
int m, n, lda, mat_size;
24+
int n, lda, mat_size;
2525
};

0 commit comments

Comments
 (0)