From e09bd0a633995ce930fcee6f7aabce72768310ae Mon Sep 17 00:00:00 2001
From: Guillaume Horel <guillaume.horel@gmail.com>
Date: Tue, 12 Dec 2017 18:04:20 -0500
Subject: do not use lapacke which is buggy

---
 src/lossdistrib.c | 31 ++++++++++++++++++-------------
 src/lossdistrib.h |  2 ++
 2 files changed, 20 insertions(+), 13 deletions(-)

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 7028e7f..a3b3866 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -5,7 +5,6 @@
 #include "lossdistrib.h"
 
 #include <cblas.h>
-#include <lapacke.h>
 
 #ifdef BUILD_TESTS
 #include <gsl/gsl_rng.h>
@@ -20,12 +19,15 @@ void GHquad(const int* n, double* Z, double* w) {
     // Initialize array for eigenvectors
     double * V      = malloc(sizeof(double)*(*n)*(*n));
 
-    for(int i = 0; i<(*n)-1; i++){
+    for(int i = 0; i < (*n) - 1; i++){
       w[i] = sqrt((i+1.)/2);
     }
 
     // Run eigen decomposition
-    int INFO = LAPACKE_dstev(LAPACK_COL_MAJOR, JOBZ, *n, Z, w, V, *n);
+    int INFO;
+    double* work = malloc((2 * (*n) - 2) * sizeof(double));
+    dstev_(&JOBZ, n, Z, w, V, n, work, &INFO);
+
 
     for(int i = 0; i<(*n); i++) {
       w[i] = V[i*(*n)] * V[i*(*n)];
@@ -34,6 +36,7 @@ void GHquad(const int* n, double* Z, double* w) {
 
     // Deallocate temporary arrays
     free(V);
+    free(work);
 }
 
 void lossdistrib(const double *p, const int *np, const double *w,
@@ -538,7 +541,7 @@ void BCloss_recov_dist(const double *defaultprob, const int *dim1, const int *di
             }
             lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, N, defaultflag,
                         Lw + i * (*N));
-            lossdistrib(gshocked, dim1, issuerweights, Rshocked , N, N, defaultflag,
+            lossdistrib(gshocked, dim1, issuerweights, Rshocked, N, N, defaultflag,
                         Rw + i * (*N));
         }
         #pragma omp parallel
@@ -648,10 +651,9 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
 
 #ifdef BUILD_TESTS
 int main() {
-    const int nZ = 500;
-    double* Z = malloc(nZ * sizeof(double));
+    const int nZ = 100;
+    double* Z = calloc(nZ, sizeof(double));
     double* w = malloc(nZ * sizeof(double));
-
     GHquad(&nZ, Z, w);
 
     const int n = 1000;
@@ -664,15 +666,12 @@ int main() {
         rho[i] = 0.45;
         p[i] = gsl_rng_uniform(gen);
     }
-
     double* S = malloc(n * nZ * sizeof(double));
-    for(int i = 0; i<nZ*n; i++) {
+    for(int i = 0; i < nZ * n; i++) {
         S[i] = gsl_rng_uniform(gen);
      }
-
     const int N = 301;
     double* q = calloc(N * nZ, sizeof(double));
-
     int flag = 0;
     lossdistrib_Z(p, &n, weights, S, &N, &flag, rho, Z, &nZ, q);
     double *qvec = malloc( N * sizeof(double));
@@ -680,11 +679,17 @@ int main() {
     cblas_dgemv(CblasColMajor, CblasNoTrans, N, nZ, alpha, q, N, w, 1, beta, qvec, 1);
     int T = 25;
     p = realloc(p, n * T * sizeof(double));
+
     for(int i = 0; i<n*T; i++) {
         p[i] = gsl_rng_uniform(gen);
     }
-    double* L = malloc(N * T * sizeof(double));
-    double* R = malloc(N * T * sizeof(double));
+    S = realloc(S, n * sizeof(double));
+    for(int i = 0; i < n; i++) {
+        S[i] = gsl_rng_uniform(gen);
+    }
+    double* L = calloc(N * T, sizeof(double));
+    double* R = calloc(N * T, sizeof(double));
+
     BCloss_recov_dist(p, &n, &T, weights, S, Z, w, &nZ, rho, &N, &flag, L, R);
     free(p);
     free(S);
diff --git a/src/lossdistrib.h b/src/lossdistrib.h
index 557333b..1132178 100644
--- a/src/lossdistrib.h
+++ b/src/lossdistrib.h
@@ -1,5 +1,7 @@
 extern void openblas_set_num_threads(int);
 
+extern void dstev_(char*, const int*, double*, double*, double*, const int*, double*, int*);
+
 void lossdistrib(const double *p, const int *np, const double *w, const double *S,
                  const int *N, const int *T, const int *defaultflag, double *q);
 
-- 
cgit v1.2.3-70-g09d2