From 544eb106db66e859a0933244bf46624eab1ab555 Mon Sep 17 00:00:00 2001 From: Guillaume Horel Date: Tue, 10 Mar 2015 14:17:51 -0400 Subject: fix lossdistrib_truncated --- src/lossdistrib.c | 38 ++++++++++++++++++-------------------- 1 file changed, 18 insertions(+), 20 deletions(-) (limited to 'src') diff --git a/src/lossdistrib.c b/src/lossdistrib.c index 61f3c84..d28b1b2 100644 --- a/src/lossdistrib.c +++ b/src/lossdistrib.c @@ -145,7 +145,7 @@ void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaul int i, j; double* pshocked = malloc(sizeof(double) * (*np) * (*nZ)); -#pragma omp parallel for private(j) + #pragma omp parallel for private(j) for(i = 0; i < *nZ; i++){ for(j = 0; j < *np; j++){ pshocked[j + (*np) * i] = shockprob(p[j], rho[j], Z[i], 0); @@ -172,34 +172,32 @@ void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N, int i, j, d1, d2, M; double lu, d, p1, p2; - double *q1, *q2; + double *qtemp; + int one = 1; + double pbar; + int bound; lu = 1./(*N-1); - q1 = calloc( *T, sizeof(double)); - q2 = calloc( *T, sizeof(double)); + qtemp = malloc( *T * sizeof(double)); q[0] = 1; M = 1; for(i=0; i<(*np); i++){ d = (*defaultflag)? w[i] / lu : S[i] * w[i] / lu; - d1 = floor(d); - d2 = ceil(d); - p1 = p[i] * (d2-d); + d1 = (int)floor(d); + d2 = d1 + 1; + p1 = p[i] * ((double)d2-d); p2 = p[i] - p1; - for(j=0; j < MIN(M, *T); j++){ - q1[j] = p1 * q[j]; - q2[j] = p2 * q[j]; - q[j] = (1-p[i]) * q[j]; - } - for(j=0; j < MIN(M, *T-d1); j++){ - q[d1+j] += q1[j]; - }; - for(j=0; j < MIN(M, *T-d2); j++){ - q[d2+j] += q2[j]; - }; + memcpy(qtemp, q, MIN(M, *T) * sizeof(double)); + pbar = 1-p[i]; + bound = MIN(M, *T); + dscal_(&bound, &pbar, q, &one); + bound = MIN(M, *T-d1); + daxpy_(&bound, &p1, qtemp, &one, q+d1, &one); + bound = MIN(M, *T-d2); + daxpy_(&bound, &p2, qtemp, &one, q+d2, &one); M += d2; } - free(q1); - free(q2); + free(qtemp); } void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) { -- cgit v1.2.3-70-g09d2