1 files changed, 75 insertions, 0 deletions
diff --git a/src/lossdistrib.hpp b/src/lossdistrib.hpp
new file mode 100644
index 0000000..0bde76a
--- /dev/null
+++ b/src/lossdistrib.hpp
@@ -0,0 +1,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;
+};