switch to cblas and lapacke interface

2 files changed, 45 insertions, 62 deletions

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 7bcfe66..7c0e997 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -4,19 +4,18 @@
 #include <omp.h>
 #include "lossdistrib.h"
 
+#include <cblas.h>
+#include <lapacke.h>
+
 #ifdef BUILD_TESTS
 #include <gsl/gsl_rng.h>
 #endif
 
 #define MIN(x, y) (((x) < (y)) ? (x) : (y))
-#define USE_BLAS
 
 void GHquad(const int* n, double* Z, double* w) {
     // Setup for eigenvalue computations
-    char JOBZ   = 'V'; // Compute eigenvalues & vectors
-    int INFO;
-    // Initialize array for workspace
-    double * WORK   = malloc(sizeof(double)*(2*(*n)-2));
+    char JOBZ   = 'v'; // Compute eigenvalues & vectors
 
     // Initialize array for eigenvectors
     double * V      = malloc(sizeof(double)*(*n)*(*n));
@@ -26,7 +25,7 @@ void GHquad(const int* n, double* Z, double* w) {
     }
 
     // Run eigen decomposition
-    dstev_(&JOBZ, n, Z, w, V, n, WORK, &INFO);
+    int INFO = LAPACKE_dstev(LAPACK_COL_MAJOR, JOBZ, *n, Z, w, V, *n);
 
     for(int i = 0; i<(*n); i++) {
       w[i] = V[i*(*n)] * V[i*(*n)];
@@ -34,7 +33,6 @@ void GHquad(const int* n, double* Z, double* w) {
     }
 
     // Deallocate temporary arrays
-    free(WORK);
     free(V);
 }
 
@@ -57,7 +55,6 @@ void lossdistrib(const double *p, const int *np, const double *w, const double *
     double lu = 1./(*N-1);
     q[0] = 1;
     int M = 1;
-    const int one = 1;
     for(int i=0; i<(*np); i++){
         d = (*defaultflag)? w[i]/lu : S[i] * w[i]/ lu;
         d1 = floor(d);
@@ -68,7 +65,7 @@ void lossdistrib(const double *p, const int *np, const double *w, const double *
         pbar = 1-p[i];
         bound = MIN(M, *T);
 #ifdef USE_BLAS
-        dscal_(&bound, &pbar, q, &one);
+        cblas_dscal(bound, pbar, q, 1);
 #else
         for(int j=0; j < bound; j++){
             q[j] *= pbar;
@@ -76,8 +73,8 @@ void lossdistrib(const double *p, const int *np, const double *w, const double *
 #endif
         bound = MIN(M, *T-d2);
 #ifdef USE_BLAS
-        daxpy_(&bound, &p1, qtemp, &one, q+d1, &one);
-        daxpy_(&bound, &p2, qtemp, &one, q+d2, &one);
+        cblas_daxpy(bound, p1, qtemp, 1, q+d1, 1);
+        cblas_daxpy(bound, p2, qtemp, 1, q+d2, 1);
 #else
         for(int j=0; j < MIN(M, *T-d2); j++){
             q[d1+j] += p1 * qtemp[j];
@@ -99,7 +96,7 @@ void lossdistrib(const double *p, const int *np, const double *w, const double *
 }
 
 void lossdistrib_Z(const double *p, const int *np, const double *w, const double *S, const int *N,
-                   const int *defaultflag, const double *rho, const double *Z, const int *nZ, double *q){
+                   const int *defaultflag, const double *rho, const double *Z, const int *nZ, double *q) {
     double* pshocked = malloc(sizeof(double) * (*np) * (*nZ));
 
     #pragma omp parallel for
@@ -113,9 +110,8 @@ void lossdistrib_Z(const double *p, const int *np, const double *w, const double
     free(pshocked);
 }
 
-static inline void posK(int T, double K, double lu, double* val){
-    int i = 0;
-    for(i = 0; i < T; i++){
+static inline void posK(int T, double K, double lu, double* val) {
+    for(int i = 0; i < T; i++) {
         val[i] = K-lu*i;
     }
 }
@@ -126,12 +122,11 @@ void exp_trunc(const double *p, const int *np, const double *w, const double *S,
     double lu = 1./(*N+1);
     int T = (int) floor((*K) * (*N))+1;
     const int flag = 0;
-    const int one = 1;
     double* qtemp = calloc( T, sizeof(double));
     lossdistrib(p, np, w, S, N, &T, &flag, qtemp);
     double* val = malloc(T * sizeof(double));
     posK(T, *K, lu, val);
-    *r = ddot_(&T, val, &one, qtemp, &one);
+    *r = cblas_ddot(T, val, 1, qtemp, 1);
     free(qtemp);
 }
 
@@ -163,9 +158,8 @@ void lossdistrib_joint(const double *p, const double* pp, const int *np,
     double* qtemp = calloc( (*N) * (*N), sizeof(double));
     q[0] = 1;
 
-    int Mx=1, My=1;
-    const int one = 1;
-    for(int k = 0; k<(*np); k++) {
+    int Mx = 1, My = 1;
+    for(int k = 0; k < (*np); k++) {
         y1 = (1-S[k]) * w[k] / lu;
         y2 = w[k]/lu;
         x = (*defaultflag)? y2 : y2 - y1;
@@ -198,7 +192,7 @@ void lossdistrib_joint(const double *p, const double* pp, const int *np,
         for(int n = 0; n < MIN(My, *N); n++) {
             memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
 #ifdef USE_BLAS
-            dscal_(&bound, &pbar, q+(*N)*n, &one);
+            cblas_dscal(bound, pbar, q+(*N)*n, 1);
 #else
             for(int m=0; m < bound; m++) {
                 q[m+(*N)*n] *= pbar;
@@ -209,20 +203,20 @@ void lossdistrib_joint(const double *p, const double* pp, const int *np,
         for(int n=0; n < MIN(My, *N-j2-1); n++) {
             begin = qtemp+(*N)*n;
 #ifdef USE_BLAS
-            daxpy_(&bound, &w1, begin, &one, q+i+(*N)*(j1+n), &one);
-            daxpy_(&bound, &w2, begin, &one, q+i+(*N)*(j1+1+n), &one);
-            daxpy_(&bound, &w3, begin, &one, q+i+1+(*N)*(j1+1+n), &one);
-            daxpy_(&bound, &w4, begin, &one, q+i+1+(*N)*(j1+n), &one);
+            cblas_daxpy(bound, w1, begin, 1, q+i+(*N)*(j1+n), 1);
+            cblas_daxpy(bound, w2, begin, 1, q+i+(*N)*(j1+1+n), 1);
+            cblas_daxpy(bound, w3, begin, 1, q+i+1+(*N)*(j1+1+n), 1);
+            cblas_daxpy(bound, w4, begin, 1, q+i+1+(*N)*(j1+n), 1);
 
             if(pp) {
                 if(defaultflag) {
-                    daxpy_(&bound, &ppw1, begin, &one, q+i+(*N)*(j2+n), &one);
-                    daxpy_(&bound, &ppw2, begin, &one, q+i+(*N)*(j2+1+n), &one);
-                    daxpy_(&bound, &ppw3, begin, &one, q+i+1+(*N)*(j2+1+n), &one);
-                    daxpy_(&bound, &ppw2, begin, &one, q+i+1+(*N)*(j2+n), &one);
+                    cblas_daxpy(bound, ppw1, begin, 1, q+i+(*N)*(j2+n), 1);
+                    cblas_daxpy(bound, ppw2, begin, 1, q+i+(*N)*(j2+1+n), 1);
+                    cblas_daxpy(bound, ppw3, begin, 1, q+i+1+(*N)*(j2+1+n), 1);
+                    cblas_daxpy(bound, ppw2, begin, 1, q+i+1+(*N)*(j2+n), 1);
                 } else {
-                    daxpy_(&bound, &ppw1, begin, &one, q+(*N)*(j2+n), &one);
-                    daxpy_(&bound, &ppw2, begin, &one, q+(*N)*(j2+1+n), &one);
+                    cblas_daxpy(bound, ppw1, begin, 1, q+(*N)*(j2+n), 1);
+                    cblas_daxpy(bound, ppw2, begin, 1, q+(*N)*(j2+1+n), 1);
                 }
             }
 #else
@@ -309,7 +303,7 @@ void recovdist(const double *dp, const double *pp, const int *n, const double *w
         M += d2u;
     }
     /* correction for weight loss */
-    if(M>*N){
+    if(M > *N){
         double sum = 0;
         for(int j=0; j<MIN(M, *N); j++){
             sum += q[j];
@@ -335,12 +329,11 @@ double dshockprob(double p, double rho, double Z){
     return( dnorm((qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1-rho), 0, 1, 0) * dqnorm(p)/sqrt(1-rho) );
 }
 
-void shockprobvec2(double p, double rho, double* Z, int nZ, double *q){
+void shockprobvec2(const double p, const double rho, const double* Z, const int nZ, double *q){
     /* return a two column vectors with shockprob in the first column
        and dshockprob in the second column*/
-    int i;
     #pragma omp parallel for
-    for(i = 0; i < nZ; i++){
+    for(int i = 0; i < nZ; i++){
         q[i] = shockprob(p, rho, Z[i], 0);
         q[i + nZ] = dshockprob(p, rho, Z[i]);
     }
@@ -354,7 +347,7 @@ double shockseverity(double S, double Z, double rho, double p){
     }
 }
 
-double quantile(double* Z, double* w, int nZ, double p0){
+double quantile(const double* Z, const double* w, const int nZ, const double p0){
     double cumw;
     int i;
     cumw = 0;
@@ -367,9 +360,9 @@ double quantile(double* Z, double* w, int nZ, double p0){
     return( Z[i] );
 }
 
-void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result){
+void fitprob(const double* Z, const double* w, const int* nZ, const double* rho,
+             const double* p0, double* result){
     double eps = 1e-12;
-    int one = 1;
     double *q = malloc( 2 * (*nZ) * sizeof(double));
     double dp, p, phi;
 
@@ -379,7 +372,7 @@ void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* res
         *result = pnorm(quantile(Z, w, *nZ, *p0), 0, 1, 1, 0);
     }else{
         shockprobvec2(*p0, *rho, Z, *nZ, q);
-        dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
+        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;
@@ -388,7 +381,7 @@ void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* res
             }
             p -= phi * dp;
             shockprobvec2(p, *rho, Z, *nZ, q);
-            dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
+            dp = (cblas_ddot(*nZ, q, 1, w, 1) - *p0)/cblas_ddot(*nZ, q + *nZ, 1, w, 1);
         }
         *result = p;
     }
@@ -398,16 +391,16 @@ void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* res
 void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ,
                      double* rho, double* porig, double* pmod, double* q){
     double ptemp, ptilde;
-    int i;
+
     if(*porig==0){
-        for(i = 0; i < *nZ; i++){
+        for(int i = 0; i < *nZ; i++){
             q[i] = *R;
         }
     }else{
         ptemp = (1 - *R) / (1 - *Rtilde) * *porig;
         fitprob(Z, w, nZ, rho, &ptemp, &ptilde);
         #pragma omp parallel for
-        for(i = 0; i < *nZ; i++){
+        for(int i = 0; i < *nZ; i++){
             q[i] = fabs(1 - (1 - *Rtilde) * exp( shockprob(ptilde, *rho, Z[i], 1) -
                                                  shockprob(*pmod, *rho, Z[i], 1)));
         }
@@ -425,7 +418,6 @@ void lossdistrib_joint_Z(const double *dp, const double *pp, const int *ndp,
 
     const double alpha = 1;
     const double beta = 0;
-    const int one = 1;
 
     #pragma omp parallel for
     for(int i = 0; i < *nZ; i++){
@@ -449,7 +441,7 @@ void lossdistrib_joint_Z(const double *dp, const double *pp, const int *ndp,
         }
     }
 
-    dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);
+    cblas_dgemv(CblasColMajor, CblasNoTrans, N2, *nZ, alpha, qmat, N2, wZ, 1, beta, q, 1);
 
     free(qmat);
 }
@@ -478,7 +470,6 @@ 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 int one = 1;
     const double alpha = 1;
     const double beta = 0;
 
@@ -494,7 +485,7 @@ void BCloss_recov_dist(const double *defaultprob, const int *dim1, const int *di
                 double g = defaultprob[j + (*dim1) * t];
                 gshocked[j] = shockprob(g, rho[j], Z[i], 0);
                 Sshocked[j] = shockseverity(1-recov[j], Z[i], rho[j], g);
-                Rshocked[j] = 1 - Sshocked[j+(*dim1)*i];
+                Rshocked[j] = 1 - Sshocked[j];
             }
             lossdistrib(gshocked, dim1, issuerweights, Sshocked, N, N, defaultflag,
                         Lw + i * (*N));
@@ -504,8 +495,8 @@ void BCloss_recov_dist(const double *defaultprob, const int *dim1, const int *di
             free(Sshocked);
             free(Rshocked);
         }
-        dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
-        dgemv_("n", N, n, &alpha, Rw, N, w, &one, &beta, R + t * (*N), &one);
+        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);
     }
     free(Lw);
     free(Rw);
@@ -536,7 +527,6 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
       ELt: vector of length dim2
       ERt: vector of length dim2
     */
-    const int one = 1;
     int T = (int) floor((*K) * (*N))+1;
     double lu = 1./(*N+1);
     double* valL = malloc( T * sizeof(double));
@@ -572,7 +562,7 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
                 lossdistrib(gshocked, dim1, issuerweights, Sshocked,
                             N, &T, defaultflag, Lw + i * T);
                 ER[i] = 0;
-                Ktilde = *K - ddot_(dim1, issuerweights, &one, Sshocked, &one);
+                Ktilde = *K - cblas_ddot(*dim1, issuerweights, 1, Sshocked, 1);
                 if(Ktilde > 0){
                     Ttilde = (int) floor(Ktilde * (*N))+1;
                     Rw = calloc(Ttilde, sizeof(double));
@@ -580,7 +570,7 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
                                 N, &Ttilde, defaultflag, Rw);
                     valR = malloc(Ttilde * sizeof(double));
                     posK(Ttilde, Ktilde, lu, valR);
-                    ER[i] = ddot_(&Ttilde, Rw, &one, valR, &one);
+                    ER[i] = cblas_ddot(Ttilde, Rw, 1, valR, 1);
                 }
                 if(Rw) {
                     free(Rw);
diff --git a/src/lossdistrib.h b/src/lossdistrib.h
index 3a65800..557333b 100644
--- a/src/lossdistrib.h
+++ b/src/lossdistrib.h
@@ -1,10 +1,3 @@
-extern int dgemv_(char* trans, const int *m, const int *n, const double* alpha, double* A, const int* lda,
-                  const double* x, const int* incx, const double* beta, double* y, const int* incy);
-extern double ddot_(const int* n, const double* dx, const int* incx, const double* dy, const int* incy);
-extern int dscal_(int* n, double* da, double* dx, const int* incx);
-extern int daxpy_(int* n, double* da, double* dx, const int* incx, double* dy, const int* incy);
-extern int dstev_(char* JOBZ, const int* n, double* D, double* E, double* Z, const int* ldz,
-                  double* WORK, int* INFO);
 extern void openblas_set_num_threads(int);
 
 void lossdistrib(const double *p, const int *np, const double *w, const double *S,
@@ -31,11 +24,11 @@ double dqnorm(double x);
 
 double dshockprob(double p, double rho, double Z);
 
-void shockprobvec2(double p, double rho, double* Z, int nZ, double *q);
+void shockprobvec2(const double p, const double rho, const double* Z, const int nZ, double *q);
 
 double shockseverity(double S, double Z, double rho, double p);
 
-void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result);
+void fitprob(const double* Z, const double* w, const int* nZ, const double* rho, const double* p0, double* result);
 
 void stochasticrecov(double* R, double* Rtilde, double* Z, double* w, int* nZ,
                      double* rho, double* porig, double* pmod, double* q);
@@ -56,6 +49,6 @@ void BCloss_recov_trunc(const double *defaultprob, const int *dim1, const int *d
                         const double * K, const int *defaultflag,
                         double *ELt, double *ERt);
 
-double quantile(double* Z, double* w, int nZ, double p0);
+double quantile(const double* Z, const double* w, const int nZ, const double p0);
 
 void GHquad(const int *n, double* Z, double* w);