1 files changed, 48 insertions, 22 deletions

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 27b4a28..78a3bc2 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -6,6 +6,19 @@
 
 #ifdef USE_BLAS
 #include <cblas.h>
+#else
+#include <R_ext/Lapack.h>
+double cblas_ddot(const int N, const double* X, const int incX, const double* Y, const int incY) {
+    double r = 0.;
+    const double* Xptr = X;
+    const double* Yptr = Y;
+    for(int i = 0; i < N; i++) {
+        r += *Xptr * *Yptr;
+        Xptr += incX;
+        Yptr += incY;
+    }
+    return r;
+}
 #endif
 
 #ifdef BUILD_TESTS
@@ -28,8 +41,11 @@ void GHquad(const int* n, double* Z, double* w) {
     // Run eigen decomposition
     int INFO;
     double* work = malloc((2 * (*n) - 2) * sizeof(double));
+    #if USE_BLAS
     dstev_(&JOBZ, n, Z, w, V, n, work, &INFO);
-
+    #else
+    F77_NAME(dstev)(&JOBZ, n, Z, w, V, n, work, &INFO);
+    #endif
 
     for(int i = 0; i<(*n); i++) {
       w[i] = V[i*(*n)] * V[i*(*n)];
@@ -54,7 +70,9 @@ void lossdistrib(const double *p, const int *np, const double *w,
        T where we truncate the distribution.
        defaultflag if true compute the default distribution
        q the loss distribution */
+    #ifdef USE_BLAS
     openblas_set_num_threads(1);
+    #endif
     int d1, d2, bound;
     double d, p1, p2, pbar;
     double *qtemp = calloc(*T, sizeof(double));
@@ -133,14 +151,7 @@ void exp_trunc(const double *p, const int *np, const double *w, const double *S,
     lossdistrib(p, np, w, S, N, &T, &flag, qtemp);
     double* val = malloc(T * sizeof(double));
     posK(T, *K, lu, val);
-    #if USE_BLAS
     *r = cblas_ddot(T, val, 1, qtemp, 1);
-    #else
-    r = 0;
-    for(int i=0; i < T; i++) {
-        *r += val[i] * qtemp[i];
-    }
-    #endif
     free(qtemp);
 }
 
@@ -434,13 +445,13 @@ void fitprob(const double* Z, const double* w, const int* nZ, const double* rho,
     double *q = malloc( 2 * (*nZ) * sizeof(double));
     double dp, p, phi;
 
-    if(*p0 == 0){
+    if(*p0 == 0) {
         *result = 0;
-    }else if(*rho == 1){
+    } else if(*rho == 1) {
         *result = pnorm(quantile(Z, w, *nZ, *p0), 0, 1, 1, 0);
-    }else{
+    } else {
         shockprobvec2(*p0, *rho, Z, *nZ, q);
-        dp = (cblas_ddot(*nZ, q, 1, w, 1) - *p0)/cblas_ddot(*nZ, q + *nZ, 1, w, 1);
+        dp = (cblas_ddot(*nZ, q, 1, w, 1) - *p0) / cblas_ddot(*nZ, q + *nZ, 1, w, 1);
         p = *p0;
         while(fabs(dp) > eps){
             phi = 1;
@@ -449,7 +460,8 @@ void fitprob(const double* Z, const double* w, const int* nZ, const double* rho,
             }
             p -= phi * dp;
             shockprobvec2(p, *rho, Z, *nZ, q);
-            dp = (cblas_ddot(*nZ, q, 1, w, 1) - *p0)/cblas_ddot(*nZ, q + *nZ, 1, w, 1);
+
+            dp = (cblas_ddot(*nZ, q, 1, w, 1) - *p0) / cblas_ddot(*nZ, q + *nZ, 1, w, 1);
         }
         *result = p;
     }
@@ -485,7 +497,7 @@ void lossdistrib_joint_Z(const double *dp, const double *pp, const int *ndp,
     double* qmat = calloc(N2 * (*nZ), sizeof(double));
 
     const double alpha = 1;
-    const double beta = 0;
+    const double beta_ = 0;
     #pragma omp parallel for
     for(int i = 0; i < *nZ; i++){
         double* dpshocked = malloc(sizeof(double) * (*ndp));
@@ -507,9 +519,12 @@ void lossdistrib_joint_Z(const double *dp, const double *pp, const int *ndp,
             free(ppshocked);
         }
     }
-
-    cblas_dgemv(CblasColMajor, CblasNoTrans, N2, *nZ, alpha, qmat, N2, wZ, 1, beta, q, 1);
-
+    #ifdef USE_BLAS
+    cblas_dgemv(CblasColMajor, CblasNoTrans, N2, *nZ, alpha, qmat, N2, wZ, 1, beta_, q, 1);
+    #else
+    int one = 1;
+    F77_NAME(dgemv)("N", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta_, q, &one);
+    #endif
     free(qmat);
 }
 
@@ -543,7 +558,7 @@ void BCloss_recov_dist(const double *defaultprob, const int *dim1, const int *di
     double *Lw = malloc((*N) * (*n) * sizeof(double));
     double *Rw = malloc((*N) * (*n) * sizeof(double));
     const double alpha = 1;
-    const double beta = 0;
+    const double beta_ = 0;
 
     for(int t = 0; t < (*dim2); t++) {
         memset(Lw, 0, (*N) * (*n) * sizeof(double));
@@ -570,8 +585,14 @@ void BCloss_recov_dist(const double *defaultprob, const int *dim1, const int *di
             }
             free(gshocked);
         }
-        cblas_dgemv(CblasColMajor, CblasNoTrans, *N, *n, alpha, Lw, *N, w, 1, beta, L + t * (*N), 1);
-        cblas_dgemv(CblasColMajor, CblasNoTrans, *N, *n, alpha, Rw, *N, w, 1, beta, R + t * (*N), 1);
+        #if USE_BLAS
+        cblas_dgemv(CblasColMajor, CblasNoTrans, *N, *n, alpha, Lw, *N, w, 1, beta_, L + t * (*N), 1);
+        cblas_dgemv(CblasColMajor, CblasNoTrans, *N, *n, alpha, Rw, *N, w, 1, beta_, R + t * (*N), 1);
+        #else
+        int one = 1;
+        F77_NAME(dgemv)("N", N, n, &alpha, Lw, N, w, &one, &beta_, L + t * (*N), &one);
+        F77_NAME(dgemv)("N", N, n, &alpha, Rw, N, w, &one, &beta_, R + t * (*N), &one);
+        #endif
     }
     free(Lw);
     free(Rw);
@@ -611,7 +632,7 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
     double* ER = malloc( (*n) * sizeof(double));
     double* Lw = malloc(T * (*n) * sizeof(double));
     const double alpha = 1;
-    const double beta = 0;
+    const double beta_ = 0;
     for(int t = 0; t < (*dim2); t++) {
         memset(Lw, 0, T * (*n) * sizeof(double));
 
@@ -655,7 +676,12 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
             }
             free(gshocked);
         }
-        cblas_dgemv(CblasColMajor, CblasTrans, T, *n, alpha, Lw, T, valL, 1, beta, EL, 1);
+        #if USE_BLAS
+        cblas_dgemv(CblasColMajor, CblasTrans, T, *n, alpha, Lw, T, valL, 1, beta_, EL, 1);
+        #else
+        int one = 1;
+        F77_NAME(dgemv)("T", &T, n, &alpha, Lw, &T, valL, &one, &beta_, EL, &one);
+        #endif
 
         ELt[t] = cblas_ddot(*n, EL, 1, w, 1);
         ERt[t] = cblas_ddot(*n, ER, 1, w, 1);