summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/lossdistrib.c94
-rw-r--r--src/lossdistrib.h13
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);
generated by cgit v1.2.3-70-g09d2 (git 2.52.0) at 2025-12-21 18:24:21 +0000