From 0997866e260e31244729466eae548a860f9a96ae Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Mon, 7 Mar 2016 15:39:21 -0500 Subject: switch to cblas and lapacke interface --- src/lossdistrib.c | 94 +++++++++++++++++++++++++------------------------------ 1 file changed, 42 insertions(+), 52 deletions(-) (limited to 'src/lossdistrib.c') 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 #include "lossdistrib.h" +#include +#include + #ifdef BUILD_TESTS #include #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 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); -- cgit v1.2.3-70-g09d2