From 0d5b7f9966ffa8cd7b2531a19f4e22e175ac9f4d Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Fri, 16 Feb 2018 12:58:38 -0500 Subject: use internal Blas and Lapack on Windows --- src/lossdistrib.c | 70 ++++++++++++++++++++++++++++++++++++++----------------- 1 file changed, 48 insertions(+), 22 deletions(-) (limited to 'src/lossdistrib.c') 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 +#else +#include +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); -- cgit v1.2.3-70-g09d2