1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
|
#include <algorithm>
extern "C" {
int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO);
int dscal_(int* n, double* da, double* dx, int* incx);
int daxpy_(int* n, double* da, double* dx, int* incx, double* dy, int* incy);
void openblas_set_num_threads(int);
int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* lda,
const double* x, int* incx, double* beta, double* y, int* incy);
double ddot_(int* n, const double* dx, int* incx, const double* dy, int* incy);
}
double shockprob(double p, double rho, double Z, bool give_log);
double dshockprob(double p, double rho, double Z);
//simple class which wraps a dynamically allocated array
//but allows it to be externally allocated
struct array {
array(int n) : n(n), extern_alloc(false) {
arr = new double[n];
}
array(int n, double* const ptr) :n(n), extern_alloc(true) {
arr = new(ptr) double[n];
}
~array() {
if(!extern_alloc) {
delete[] arr;
}
}
double* begin() {
return arr;
}
double* begin() const {
return arr;
}
double* end() {
return arr+n;
}
double* end() const {
return arr+n;
}
private:
const int n;
double* arr;
const bool extern_alloc;
};
struct matrix {
matrix(int m, int n) : m(m), n(n), mat(m*n) {};
matrix(int m, int n, double x) : m(m), n(n), mat(m*n) {
std::fill(mat.begin(), mat.end(), x);
}
matrix(int m, int n, double* ptr) : m(m), n(n), mat(m*n, ptr) {};
double& operator()(int i, int j) {
return *(mat.begin()+j*m+i);
}
double operator()(int i, int j) const {
return *(mat.begin()+j*m+i);
}
double* operator()(int i) {
return mat.begin()+i*m;
}
double* data() {
return mat.begin();
}
double& operator[](int i) {
return *(mat.begin()+i);
}
double operator[](int i) const {
return *(mat.begin()+i);
}
private:
int m, n;
array mat;
};
|
