diff options
| -rw-r--r-- | src/lossdistrib.c | 22 |
1 files changed, 22 insertions, 0 deletions
diff --git a/src/lossdistrib.c b/src/lossdistrib.c index 494fed7..27b4a28 100644 --- a/src/lossdistrib.c +++ b/src/lossdistrib.c @@ -4,7 +4,9 @@ #include <omp.h> #include "lossdistrib.h" +#ifdef USE_BLAS #include <cblas.h> +#endif #ifdef BUILD_TESTS #include <gsl/gsl_rng.h> @@ -131,7 +133,14 @@ 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); } @@ -292,14 +301,27 @@ void joint_default_averagerecov_distrib(const double *p,const int *np, index[n] = (unsigned int) temp; } if(j < i+1) { + #if USE_BLAS cblas_dscal(*N, 1-p[i], q+(*N)*j, 1); + #else + double* temp = q + (*N) * j; + for( int k = 0; k < *N; k++) { + temp[k] *= 1 - p[i]; + } + #endif } for(int k = 0; k < *N; k++) { q[j*(*N)+index[k]+1] += weights[k] * p[i] * q[(j-1)*(*N)+k]; q[j*(*N)+index[k]] += (1-weights[k]) * p[i] * q[(j-1)*(*N)+k]; } } + #if USE_BLAS cblas_dscal(*N, 1-p[i], q, 1); + #else + for(int k = 0; k < *N; k++) { + q[k] *= 1 - p[i]; + } + #endif } free(newrecov); free(index); |