1 files changed, 62 insertions, 3 deletions

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index f6ff411..7bcfe66 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -4,6 +4,10 @@
 #include <omp.h>
 #include "lossdistrib.h"
 
+#ifdef BUILD_TESTS
+#include <gsl/gsl_rng.h>
+#endif
+
 #define MIN(x, y) (((x) < (y)) ? (x) : (y))
 #define USE_BLAS
 
@@ -590,12 +594,67 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
                 free(gshockedbar);
             }
         }
-        dgemv_("t", &T, n, &alpha, Lw, &T, valL, &one, &beta, EL, &one);
-        ELt[t] = ddot_(n, EL, &one, w, &one);
-        ERt[t] = ddot_(n, ER, &one, w, &one);
+        cblas_dgemv(CblasColMajor, CblasTrans, T, *n, alpha, Lw, T, valL, 1, beta, EL, 1);
+        ELt[t] = cblas_ddot(*n, EL, 1, w, 1);
+        ERt[t] = cblas_ddot(*n, ER, 1, w, 1);
     }
     free(Lw);
     free(EL);
     free(ER);
     free(valL);
 }
+
+#ifdef BUILD_TESTS
+int main() {
+    const int nZ = 500;
+    double* Z = malloc(nZ * sizeof(double));
+    double* w = malloc(nZ * sizeof(double));
+
+    GHquad(&nZ, Z, w);
+
+    const int n = 1000;
+    double* p = malloc(n * sizeof(double));
+    double* weights = malloc(n * sizeof(double));
+    double* rho = malloc(n * sizeof(double));
+    gsl_rng* gen = gsl_rng_alloc(gsl_rng_mt19937);
+    for(int i = 0; i<n; i++) {
+        weights[i] = 1./n;
+        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++) {
+        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));
+    double alpha=1, beta=0;
+    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));
+    BCloss_recov_dist(p, &n, &T, weights, S, Z, w, &nZ, rho, &N, &flag, L, R);
+    free(p);
+    free(S);
+    free(weights);
+    free(Z);
+    free(w);
+    free(q);
+    free(rho);
+    free(L);
+    free(R);
+    free(qvec);
+    gsl_rng_free(gen);
+    return 0;
+}
+#endif