From c85640622f4686db60f3c6b469e007798f43ad53 Mon Sep 17 00:00:00 2001
From: Guillaume Horel <guillaume.horel@serenitascapital.com>
Date: Tue, 10 Mar 2015 12:55:29 -0400
Subject: reindent everything

---
 src/lossdistrib.c | 1682 ++++++++++++++++++++++++++---------------------------
 1 file changed, 841 insertions(+), 841 deletions(-)

(limited to 'src')

diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index e5f2c53..61f3c84 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -1,841 +1,841 @@
-#include <R.h>
-#include <Rmath.h>
-#include <string.h>
-#include <omp.h>
-#include "lossdistrib.h"
-
-#define MIN(x, y) (((x) < (y)) ? (x) : (y))
-
-extern int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* lda,
-                  double* x, int* incx, double* beta, double* y, int* incy);
-extern double ddot_(int* n, double* dx, int* incx, double* dy, int* incy);
-extern int dscal_(int* n, double* da, double* dx, int* incx);
-extern int daxpy_(int* n, double* da, double* dx, int* incx, double* dy, int* incy);
-extern int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO);
-
-extern void openblas_set_num_threads(int);
-
-void GHquad(int* n, double* Z, double* w) {
-    // Setup for eigenvalue computations
-    char JOBZ   = 'V'; // Compute eigenvalues & vectors
-    int INFO;
-    int i;
-    // Initialize array for workspace
-    double * WORK   = malloc(sizeof(double)*(2*(*n)-2));
-
-    // Initialize array for eigenvectors
-    double * V      = malloc(sizeof(double)*(*n)*(*n));
-
-    for(i = 0; i<(*n)-1; i++){
-      w[i] = sqrt((i+1.)/2);
-    }
-
-    // Run eigen decomposition
-    dstev_(&JOBZ, n, Z, w, V, n, WORK, &INFO);
-
-    for (i=0; i<(*n); i++) {
-      w[i] = V[i*(*n)] * V[i*(*n)];
-      Z[i] *= sqrt(2);
-    }
-
-    // Deallocate temporary arrays
-    free(WORK);
-    free(V);
-}
-
-void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
-		 double *q) {
-  /* recursive algorithm with first order correction for computing
-     the loss distribution.
-     p vector of default probabilities
-     np length of p
-     w issuer weights
-     S vector of severities (should be same length as p)
-     N number of ticks in the grid
-     defaultflag if true compute the default distribution
-     q the loss distribution */
-
-  int i, j, d1, d2, M;
-  double lu, d, p1, p2, sum;
-  double *qtemp;
-
-  lu = 1./(*N-1);
-  qtemp = calloc(*N, 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);
-    p2 = p[i] - p1;
-    memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
-    for(j=0; j < MIN(M, *N); j++){
-      q[j] = (1-p[i]) * q[j];
-    }
-    for(j=0; j < MIN(M, *N-d2); j++){
-      q[d1+j] += p1 * qtemp[j];
-      q[d2+j] += p2 * qtemp[j];
-    };
-    M+=d2;
-  }
-
-  /* correction for weight loss */
-  if(M > *N){
-    sum = 0;
-    for(j=0; j<MIN(M, *N); j++){
-      sum += q[j];
-    }
-    q[*N-1] += 1-sum;
-  }
-  free(qtemp);
-}
-
-void lossdistrib_blas(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
-		      double *q) {
-  /* recursive algorithm with first order correction for computing
-     the loss distribution.
-     p: vector of default probabilities
-     np: length of p
-     w: issuer weights
-     S: vector of severities (should be same length as p)
-     N: number of ticks in the grid
-     defaultflag: if true compute the default distribution
-     q: the loss distribution */
-
-  int i, j, d1, d2, M;
-  double lu, d, p1, p2, sum;
-  double *qtemp;
-  int bound;
-  double pbar;
-  int one = 1;
-  openblas_set_num_threads(1);
-  lu = 1./(*N-1);
-  qtemp = calloc(*N, 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);
-    p2 = p[i] - p1;
-    memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
-    pbar = 1-p[i];
-    bound = MIN(M, *N);
-    dscal_(&bound, &pbar, q, &one);
-    bound = MIN(M, *N-d2);
-    daxpy_(&bound, &p1, qtemp, &one, q+d1, &one);
-    daxpy_(&bound, &p2, qtemp, &one, q+d2, &one);
-    M += d2;
-  }
-  /* correction for weight loss */
-  if(M > *N){
-    sum = 0;
-    for(j=0; j<MIN(M, *N); j++){
-      sum += q[j];
-    }
-    q[*N-1] += 1-sum;
-  }
-  free(qtemp);
-}
-
-void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
-                   double *rho, double *Z, int *nZ, double *q){
-  int i, j;
-  double* pshocked = malloc(sizeof(double) * (*np) * (*nZ));
-
-  #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);
-    }
-    lossdistrib_blas(pshocked + (*np) * i, np, w, S + (*np) * i, N,
-		     defaultflag, q + (*N) * i);
-  }
-  free(pshocked);
-}
-
-void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
-                           int *T, int *defaultflag, double *q) {
-  /* recursive algorithm with first order correction for computing
-     the loss distribution.
-     input:
-     p vector of default probabilities
-     np length of p
-     S vector of severities (should be same length as p)
-     N number of ticks in the grid
-     T where to truncate the distribution
-     defaultflag if true computes the default distribution
-     output:
-     q the loss distribution */
-
-  int i, j, d1, d2, M;
-  double lu, d, p1, p2;
-  double *q1, *q2;
-
-  lu = 1./(*N-1);
-  q1 = calloc( *T, sizeof(double));
-  q2 = calloc( *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);
-    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];
-    };
-    M += d2;
-  }
-  free(q1);
-  free(q2);
-}
-
-void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
-  /* recursive algorithm with first order correction
-     computes jointly the loss and recovery distribution
-     p vector of default probabilities
-     np length of p
-     w vector of issuer weights (length np)
-     S vector of severities (should be same length as p)
-     N number of ticks in the grid
-     defaultflag if true computes the default distribution
-     q the joint probability distribution */
-
-  int i, j, k, m, n;
-  int Mx, My;
-  double lu, x, y, sum;
-  double alpha1, alpha2, beta1, beta2;
-  double w1, w2, w3, w4;
-  double *qtemp;
-
-  lu = 1./(*N-1);
-  qtemp = calloc( (*N) * (*N), sizeof(double));
-  q[0] = 1;
-  Mx=1;
-  My=1;
-  for(k=0; k<(*np); k++){
-    x = (*defaultflag)? w[k] /lu : S[k] * w[k] / lu;
-    y = (1-S[k]) * w[k] / lu;
-    i = floor(x);
-    j = floor(y);
-    alpha1 = i + 1 - x;
-    alpha2 = 1 - alpha1;
-    beta1 = j + 1 - y;
-    beta2 = 1 - beta1;
-    w1 = alpha1 * beta1;
-    w2 = alpha1 * beta2;
-    w3 = alpha2 * beta2;
-    w4 = alpha2 * beta1;
-
-    for(n=0; n<MIN(My, *N); n++){
-      memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
-    }
-    for(n=0; n<MIN(My, *N); n++){
-      for(m=0; m<MIN(Mx, *N); m++){
-	q[m+(*N)*n] = (1-p[k])* q[m+(*N)*n];
-      }
-    }
-    for(n=0; n < MIN(My, *N-j-1); n++){
-      for(m=0; m < MIN(Mx, *N-i-1); m++){
-	q[i+m+(*N)*(j+n)] += w1 * p[k] * qtemp[m+(*N)*n];
-	q[i+m+(*N)*(j+1+n)] += w2 * p[k] * qtemp[m+(*N)*n];
-	q[i+1+m+(*N)*(j+1+n)] += w3 * p[k] * qtemp[m+(*N)*n];
-	q[i+1+m+(*N)*(j+n)] += w4 * p[k] *qtemp[m+(*N)*n];
-      }
-    }
-    Mx += i+1;
-    My += j+1;
-  }
-  /* correction for weight loss */
-  if(Mx>*N || My>*N){
-    sum = 0;
-    for(m=0; m < MIN(Mx, *N); m++){
-      for(n=0; n < MIN(My, *N); n++){
-  	sum += q[m+n*(*N)];
-      }
-    }
-    q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
-  }
-  free(qtemp);
-}
-
-void lossdistrib_joint_blas(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
-  /* recursive algorithm with first order correction
-     computes jointly the loss and recovery distribution
-     p vector of default probabilities
-     np length of p
-     w vector of issuer weights (length np)
-     S vector of severities (should be same length as p)
-     N number of ticks in the grid
-     defaultflag if true computes the default distribution
-     q the joint probability distribution */
-
-  int i, j, k, m, n;
-  int Mx, My;
-  double lu, x, y, sum, pbar;
-  double alpha1, alpha2, beta1, beta2;
-  double w1, w2, w3, w4;
-  double *qtemp;
-  int bound;
-  int one = 1;
-
-  /* only use one thread, performance is horrible if use multiple threads */
-  openblas_set_num_threads(1);
-
-  lu = 1./(*N-1);
-  qtemp = calloc( (*N) * (*N), sizeof(double));
-  q[0] = 1;
-  Mx=1;
-  My=1;
-  for(k=0; k<(*np); k++){
-    x = (*defaultflag)? w[k] /lu : S[k] * w[k] / lu;
-    y = (1-S[k]) * w[k] / lu;
-    i = floor(x);
-    j = floor(y);
-    alpha1 = i + 1 - x;
-    alpha2 = 1 - alpha1;
-    beta1 = j + 1 - y;
-    beta2 = 1 - beta1;
-    w1 = alpha1 * beta1 * p[k];
-    w2 = alpha1 * beta2 * p[k];
-    w3 = alpha2 * beta2 * p[k];
-    w4 = alpha2 * beta1 * p[k];
-
-    for(n=0; n<MIN(My, *N); n++){
-      memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
-    }
-
-    bound = MIN(Mx, *N);
-    pbar = 1-p[k];
-    for(n=0; n<MIN(My, *N); n++){
-      dscal_(&bound, &pbar, q+(*N)*n, &one);
-    }
-    bound = MIN(Mx, *N-i-1);
-    for(n=0; n < MIN(My, *N-j-1); n++){
-      daxpy_(&bound, &w1, qtemp+(*N)*n, &one, q+i+(*N)*(j+n), &one);
-      daxpy_(&bound, &w2, qtemp+(*N)*n, &one, q+i+(*N)*(j+1+n), &one);
-      daxpy_(&bound, &w3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j+1+n), &one);
-      daxpy_(&bound, &w4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j+n), &one);
-    }
-    Mx += i+1;
-    My += j+1;
-  }
-  /* correction for weight loss */
-  if(Mx>*N || My>*N){
-    sum = 0;
-    for(m=0; m < MIN(Mx, *N); m++){
-      for(n=0; n < MIN(My, *N); n++){
-  	sum += q[m+n*(*N)];
-      }
-    }
-    q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
-  }
-  free(qtemp);
-}
-
-void recovdist(double *dp, double *pp, int *n, double *w, double *S, int *N, double *q) {
-  /* recursive algorithm with first order correction for computing
-     the recovery distribution in case of prepayment.
-     dp vector of default probabilities
-     pp vector of prepay probabilities
-     n length of p
-     S vector of severities (should be same length as p)
-     w vector of weights
-     N number of ticks in the grid
-     q the loss distribution */
-
-  int i, j, d1l, d1u, d2l, d2u;
-  int M;
-  double lu, d1, d2, dp1, dp2, pp1, pp2, sum;
-  double *qtemp;
-
-  lu = 1./(*N - 1);
-  qtemp = calloc( (*N), sizeof(double));
-  q[0] = 1;
-  M=1;
-  for(i=0; i<(*n); i++){
-    d1 = w[i] * (1-S[i]) /lu;
-    d2 = w[i]/lu;
-    d1l = floor(d1);
-    d1u = d1l + 1;
-    d2l = floor(d2);
-    d2u = d2l + 1;
-    dp1 = dp[i] * (d1u - d1);
-    dp2 = dp[i] - dp1;
-    pp1 = pp[i] * (d2u - d2);
-    pp2 = pp[i] - pp1;
-    memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
-    for(j = 0; j< MIN(M, *N); j++){
-      q[j] = (1-dp[i]-pp[i]) * q[j];
-    }
-    for(j=0; j < MIN(M, *N-d2u); j++){
-      q[d1l+j] += dp1 * qtemp[j];
-      q[d1u+j] += dp2 * qtemp[j];
-      q[d2l+j] += pp1 * qtemp[j];
-      q[d2u+j] += pp2 * qtemp[j];
-    };
-    M += d2u;
-  }
-  /* correction for weight loss */
-  if(M>*N){
-    sum = 0;
-    for(j=0; j<MIN(M, *N); j++){
-      sum += q[j];
-    }
-    q[*N-1] += 1-sum;
-  }
-  free(qtemp);
-}
-
-void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
-			      double *S, int *N, int *defaultflag, double *q) {
-  int i, j1, j2, k, m, n;
-  double lu, x, y1, y2, sum;
-  double alpha1, alpha2, beta1, beta2;
-  double dpw1, dpw2, dpw3, dpw4;
-  double ppw1, ppw2, ppw3;
-  double *qtemp;
-  int Mx, My;
-
-  lu = 1./(*N-1);
-  qtemp = calloc((*N) * (*N), sizeof(double));
-  q[0] = 1;
-  Mx=1;
-  My=1;
-
-  for(k=0; k<(*ndp); k++){
-    y1 = (1-S[k]) * w[k]/lu;
-    y2 = w[k]/lu;
-    x = (*defaultflag)? y2: y2-y1;
-    i = floor(x);
-    j1 = floor(y1);
-    j2 = floor(y2);
-    alpha1 = i + 1 - x;
-    alpha2 = 1 - alpha1;
-    beta1 = j1 + 1 - y1;
-    beta2 = 1 - beta1;
-    dpw1 = alpha1 * beta1 * dp[k];
-    dpw2 = alpha1 * beta2 * dp[k];
-    dpw3 = alpha2 * beta2 * dp[k];
-    dpw4 = alpha2 * beta1 * dp[k];
-
-    /* by default distribution, we mean names fractions of names that disappeared
-       either because of default or prepayment */
-    for(n=0; n<MIN(My, *N); n++){
-      memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
-    }
-    for(n=0; n<MIN(My, *N); n++){
-      for(m=0; m<MIN(Mx, *N); m++){
-	q[m+(*N)*n] = (1-dp[k]-pp[k]) * q[m+(*N)*n];
-      }
-    }
-    if(*defaultflag){
-      ppw1 = alpha1 * alpha1 * pp[k];
-      ppw2 = alpha1 * alpha2 * pp[k];
-      ppw3 = alpha2 * alpha2 * pp[k];
-      for(n=0; n < MIN(My, *N-j2-1); n++){
-	for(m=0; m < MIN(Mx, *N-i-1); m++){
-	  q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
-	  q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
-	  q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
-	  q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
-
-	  q[i+m+(*N)*(j2+n)] += ppw1 * qtemp[m+(*N)*n];
-	  q[i+m+(*N)*(j2+1+n)] += ppw2 * qtemp[m+(*N)*n];
-	  q[i+1+m+(*N)*(j2+1+n)] += ppw3 * qtemp[m+(*N)*n];
-	  q[i+1+m+(*N)*(j2+n)] += ppw2 * qtemp[m+(*N)*n];
-	}
-      }
-    }else{
-      for(n=0; n < MIN(My, *N-j2-1); n++){
-	for(m=0; m < MIN(Mx, *N-i-1); m++){
-	  q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
-	  q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
-	  q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
-	  q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
-	  q[m+(*N)*(j2+n)] +=  pp[k]*(j2+1-y2) * qtemp[m+(*N)*n];
-	  q[m+(*N)*(j2+1+n)] +=  pp[k]*(y2-j2) * qtemp[m+(*N)*n];
-	}
-      }
-    }
-    Mx += i + 1;
-    My += j2 + 1;
-  }
-  /* correction for weight loss */
-  if(Mx>*N || My>*N){
-    sum = 0;
-    for(m=0; m < MIN(Mx, *N); m++){
-      for(n=0; n < MIN(My, *N); n++){
-	sum += q[m+n*(*N)];
-      }
-    }
-    q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
-  }
-  free(qtemp);
-}
-
-void lossdistrib_prepay_joint_blas(double *dp, double *pp, int *ndp, double *w,
-				   double *S, int *N, int *defaultflag, double *q) {
-  int i, j1, j2, k, m, n;
-  double lu, x, y1, y2, sum;
-  double alpha1, alpha2, beta1, beta2;
-  double dpw1, dpw2, dpw3, dpw4;
-  double ppw1, ppw2, ppw3;
-  double *qtemp;
-  int Mx, My, bound;
-  double pbar;
-  int one = 1;
-
-  lu = 1./(*N-1);
-  qtemp = calloc((*N) * (*N), sizeof(double));
-  q[0] = 1;
-  Mx=1;
-  My=1;
-
-  /* only use one thread, performance is horrible if use multiple threads */
-  openblas_set_num_threads(1);
-  for(k=0; k<(*ndp); k++){
-    y1 = (1-S[k]) * w[k]/lu;
-    y2 = w[k]/lu;
-    x = (*defaultflag)? y2: y2-y1;
-    i = floor(x);
-    j1 = floor(y1);
-    j2 = floor(y2);
-    alpha1 = i + 1 - x;
-    alpha2 = 1 - alpha1;
-    beta1 = j1 + 1 - y1;
-    beta2 = 1 - beta1;
-    dpw1 = alpha1 * beta1 * dp[k];
-    dpw2 = alpha1 * beta2 * dp[k];
-    dpw3 = alpha2 * beta2 * dp[k];
-    dpw4 = alpha2 * beta1 * dp[k];
-
-    /* by default distribution, we mean names fractions of names that disappeared
-       either because of default or prepayment */
-    for(n=0; n<MIN(My, *N); n++){
-      memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
-    }
-    bound = MIN(Mx, *N);
-    pbar = 1-dp[k]-pp[k];
-    for(n=0; n<MIN(My, *N); n++){
-      dscal_(&bound, &pbar, q+(*N)*n, &one);
-    }
-    bound = MIN(Mx, *N-i-1);
-    if(*defaultflag){
-      ppw1 = alpha1 * alpha1 * pp[k];
-      ppw2 = alpha1 * alpha2 * pp[k];
-      ppw3 = alpha2 * alpha2 * pp[k];
-      for(n=0; n < MIN(My, *N-j2-1); n++){
-	daxpy_(&bound, &dpw1, qtemp+(*N)*n, &one, q+i+(*N)*(j1+n), &one);
-	daxpy_(&bound, &dpw2, qtemp+(*N)*n, &one, q+i+(*N)*(j1+1+n), &one);
-	daxpy_(&bound, &dpw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+1+n), &one);
-	daxpy_(&bound, &dpw4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+n), &one);
-
-	daxpy_(&bound, &ppw1, qtemp+(*N)*n, &one, q+i+(*N)*(j2+n), &one);
-	daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+i+(*N)*(j2+1+n), &one);
-	daxpy_(&bound, &ppw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j2+1+n), &one);
-	daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+i+1+(*N)*(j2+n), &one);
-      }
-    }else{
-      ppw1 = pp[k] * (j2+1-y2);
-      ppw2 = pp[k] * (y2-j2);
-      for(n=0; n < MIN(My, *N-j2-1); n++){
-	daxpy_(&bound, &dpw1, qtemp+(*N)*n, &one, q+i+(*N)*(j1+n), &one);
-	daxpy_(&bound, &dpw2, qtemp+(*N)*n, &one, q+i+(*N)*(j1+1+n), &one);
-	daxpy_(&bound, &dpw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+1+n), &one);
-	daxpy_(&bound, &dpw4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+n), &one);
-	daxpy_(&bound, &ppw1, qtemp+(*N)*n, &one, q+(*N)*(j2+n), &one);
-	daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+(*N)*(j2+1+n), &one);
-      }
-    }
-    Mx += i + 1;
-    My += j2 + 1;
-  }
-  /* correction for weight loss */
-  if(Mx>*N || My>*N){
-    sum = 0;
-    for(m=0; m < MIN(Mx, *N); m++){
-      for(n=0; n < MIN(My, *N); n++){
-	sum += q[m+n*(*N)];
-      }
-    }
-    q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
-  }
-  free(qtemp);
-}
-
-double shockprob(double p, double rho, double Z, int give_log){
-  if(rho==1){
-    return((double)(Z<=qnorm(p, 0, 1, 1, 0)));
-  }else{
-    return( pnorm( (qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
-  }
-}
-
-double dqnorm(double x){
-  return 1/dnorm(qnorm(x, 0, 1, 1, 0), 0, 1, 0);
-}
-
-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){
-  /* 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++){
-    q[i] = shockprob(p, rho, Z[i], 0);
-    q[i + nZ] = dshockprob(p, rho, Z[i]);
-  }
-}
-
-double shockseverity(double S, double Z, double rho, double p){
-  if(p==0){
-    return 0;
-  }else{
-    return( exp(shockprob(S * p, rho, Z, 1) - shockprob(p, rho, Z, 1)) );
-  }
-}
-
-double quantile(double* Z, double* w, int nZ, double p0){
-  double cumw;
-  int i;
-  cumw = 0;
-  for(i=0; i<nZ; i++){
-    cumw += w[i];
-    if(cumw > p0){
-      break;
-    }
-  }
-  return( Z[i] );
-}
-
-void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result){
-  double eps = 1e-12;
-  int one = 1;
-  double *q = malloc( 2 * (*nZ) * sizeof(double));
-  double dp, p, phi;
-
-  if(*p0 == 0){
-    *result = 0;
-  }else if(*rho == 1){
-    *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);
-    p = *p0;
-    while(fabs(dp) > eps){
-      phi = 1;
-      while( (p - phi * dp) < 0 || (p - phi * dp) > 1){
-	phi *= 0.8;
-      }
-      p -= phi * dp;
-      shockprobvec2(p, *rho, Z, *nZ, q);
-      dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
-    }
-    *result = p;
-  }
-  free(q);
-}
-
-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++){
-      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++){
-      q[i] = fabs(1 - (1 - *Rtilde) * exp( shockprob(ptilde, *rho, Z[i], 1) -
-                                           shockprob(*pmod, *rho, Z[i], 1)));
-    }
-  }
-}
-
-void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
-                                double *S, int *N, int *defaultflag, double *rho,
-                                double *Z, double *wZ, int *nZ, double *q) {
-  int i, j;
-  double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
-  double* ppshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
-  int N2 = (*N) * (*N);
-  double* qmat = malloc(sizeof(double) * N2 * (*nZ));
-
-  double alpha = 1;
-  double beta = 0;
-  int one = 1;
-
-#pragma omp parallel for private(j)
-  for(i = 0; i < *nZ; i++){
-    for(j = 0; j < *ndp; j++){
-      dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
-      ppshocked[j + (*ndp) * i] = shockprob(pp[j], rho[j], -Z[i], 0);
-    }
-    lossdistrib_prepay_joint_blas(dpshocked + (*ndp) * i, ppshocked + (*ndp) * i, ndp,
-				  w, S + (*ndp) * i, N, defaultflag, qmat + N2 * i);
-  }
-
-  dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);
-
-  free(dpshocked);
-  free(ppshocked);
-  free(qmat);
-}
-
-void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
-                         double *S, int *N, int *defaultflag, double *rho,
-                         double *Z, double *wZ, int *nZ, double *q) {
-  int i, j;
-  double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
-  int N2 = (*N) * (*N);
-  double* qmat = malloc(sizeof(double) * N2 * (*nZ));
-
-  double alpha = 1;
-  double beta = 0;
-  int one = 1;
-
-#pragma omp parallel for private(j)
-  for(i = 0; i < *nZ; i++){
-    for(j = 0; j < *ndp; j++){
-      dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
-    }
-    lossdistrib_joint_blas(dpshocked + (*ndp) * i, ndp, w, S + (*ndp) * i, N,
-			   defaultflag, qmat + N2 * i);
-  }
-
-  dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);
-
-  free(dpshocked);
-  free(qmat);
-}
-
-void BCloss_recov_dist(double *defaultprob, int *dim1, int *dim2,
-		       double *issuerweights, double *recov, double *Z, double *w,
-		       int *n, double *rho, int *N, int *defaultflag,
-		       double *L, double *R) {
-  /*
-    computes the loss and recovery distribution over time with a flat gaussian
-    correlation
-    inputs:
-      defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers
-      and dim2 number of time steps
-      issuerweights: vector of issuer weights (length dim1)
-      recov: vector of recoveries (length dim1)
-      Z: vector of factor values (length n)
-      w: vector of factor weights (length n)
-      rho: correlation beta vector (length dim1)
-      N: number of ticks in the grid
-      defaultflag: if true, computes the default distribution
-    outputs:
-      L: matrix of size (N, dim2)
-      R: matrix of size (N, dim2)
-   */
-  int t, i, j;
-  double g;
-  double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
-  int one = 1;
-  double alpha = 1;
-  double beta = 0;
-
-  gshocked = malloc((*dim1) * (*n) * sizeof(double));
-  Rshocked = malloc((*dim1) * (*n) * sizeof(double));
-  Sshocked = malloc((*dim1) * (*n) * sizeof(double));
-  Lw = malloc((*N) * (*n) * sizeof(double));
-  Rw = malloc((*N) * (*n) * sizeof(double));
-
-
-  for(t=0; t < (*dim2); t++) {
-    memset(Lw, 0, (*N) * (*n) * sizeof(double));
-    memset(Rw, 0, (*N) * (*n) * sizeof(double));
-    #pragma omp parallel for private(j, g)
-    for(i=0; i < *n; i++){
-      for(j=0; j < (*dim1); j++){
-        g = defaultprob[j + (*dim1) * t];
-        gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
-        Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
-        Rshocked[j+(*dim1)*i] = 1 - Sshocked[j+(*dim1)*i];
-      }
-      lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
-		  Lw + i * (*N));
-      lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Rshocked + (*dim1)*i, N, defaultflag,
-		  Rw + i * (*N));
-    }
-    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);
-  }
-  free(gshocked);
-  free(Rshocked);
-  free(Sshocked);
-  free(Lw);
-  free(Rw);
-}
-
-
-void BCloss_dist(double *defaultprob, int *dim1, int *dim2,
-      double *issuerweights, double *recov, double *Z, double *w,
-      int *n, double *rho, int *N, int *defaultflag,
-      double *L) {
-  /*
-    computes the loss distribution over time with a flat gaussian
-    correlation
-    inputs:
-      defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers
-      and dim2 number of time steps
-      issuerweights: vector of issuer weights (length dim1)
-      recov: vector of recoveries (length dim1)
-      Z: vector of factor values (length n)
-      w: vector of factor weights (length n)
-      rho: correlation beta vector (length dim1)
-      N: number of ticks in the grid
-      defaultflag: if true, computes the default distribution
-    outputs:
-      L: matrix of size (N, dim2)
-   */
-  int t, i, j;
-  double g;
-  double *gshocked, *Sshocked, *Lw;
-  int one = 1;
-  double alpha = 1;
-  double beta = 0;
-
-  gshocked = malloc((*dim1) * (*n) * sizeof(double));
-  Sshocked = malloc((*dim1) * (*n) * sizeof(double));
-  Lw = malloc((*N) * (*n) * sizeof(double));
-
-  for(t=0; t < (*dim2); t++) {
-    memset(Lw, 0, (*N) * (*n) * sizeof(double));
-    #pragma omp parallel for private(j, g)
-    for(i=0; i < *n; i++){
-      for(j=0; j < (*dim1); j++){
-        g = defaultprob[j + (*dim1) * t];
-        gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
-        Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
-      }
-      lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
-		       Lw + i * (*N));
-    }
-    dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
-  }
-  free(gshocked);
-  free(Sshocked);
-  free(Lw);
-}
+#include <R.h>
+#include <Rmath.h>
+#include <string.h>
+#include <omp.h>
+#include "lossdistrib.h"
+
+#define MIN(x, y) (((x) < (y)) ? (x) : (y))
+
+extern int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* lda,
+                  double* x, int* incx, double* beta, double* y, int* incy);
+extern double ddot_(int* n, double* dx, int* incx, double* dy, int* incy);
+extern int dscal_(int* n, double* da, double* dx, int* incx);
+extern int daxpy_(int* n, double* da, double* dx, int* incx, double* dy, int* incy);
+extern int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO);
+
+extern void openblas_set_num_threads(int);
+
+void GHquad(int* n, double* Z, double* w) {
+    // Setup for eigenvalue computations
+    char JOBZ   = 'V'; // Compute eigenvalues & vectors
+    int INFO;
+    int i;
+    // Initialize array for workspace
+    double * WORK   = malloc(sizeof(double)*(2*(*n)-2));
+
+    // Initialize array for eigenvectors
+    double * V      = malloc(sizeof(double)*(*n)*(*n));
+
+    for(i = 0; i<(*n)-1; i++){
+      w[i] = sqrt((i+1.)/2);
+    }
+
+    // Run eigen decomposition
+    dstev_(&JOBZ, n, Z, w, V, n, WORK, &INFO);
+
+    for (i=0; i<(*n); i++) {
+      w[i] = V[i*(*n)] * V[i*(*n)];
+      Z[i] *= sqrt(2);
+    }
+
+    // Deallocate temporary arrays
+    free(WORK);
+    free(V);
+}
+
+void lossdistrib(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
+                 double *q) {
+    /* recursive algorithm with first order correction for computing
+       the loss distribution.
+       p vector of default probabilities
+       np length of p
+       w issuer weights
+       S vector of severities (should be same length as p)
+       N number of ticks in the grid
+       defaultflag if true compute the default distribution
+       q the loss distribution */
+
+    int i, j, d1, d2, M;
+    double lu, d, p1, p2, sum;
+    double *qtemp;
+
+    lu = 1./(*N-1);
+    qtemp = calloc(*N, 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);
+        p2 = p[i] - p1;
+        memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
+        for(j=0; j < MIN(M, *N); j++){
+            q[j] = (1-p[i]) * q[j];
+        }
+        for(j=0; j < MIN(M, *N-d2); j++){
+            q[d1+j] += p1 * qtemp[j];
+            q[d2+j] += p2 * qtemp[j];
+        }
+        M+=d2;
+    }
+
+    /* correction for weight loss */
+    if(M > *N){
+        sum = 0;
+        for(j=0; j<MIN(M, *N); j++){
+            sum += q[j];
+        }
+        q[*N-1] += 1-sum;
+    }
+    free(qtemp);
+}
+
+void lossdistrib_blas(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
+                      double *q) {
+    /* recursive algorithm with first order correction for computing
+       the loss distribution.
+       p: vector of default probabilities
+       np: length of p
+       w: issuer weights
+       S: vector of severities (should be same length as p)
+       N: number of ticks in the grid
+       defaultflag: if true compute the default distribution
+       q: the loss distribution */
+
+    int i, j, d1, d2, M;
+    double lu, d, p1, p2, sum;
+    double *qtemp;
+    int bound;
+    double pbar;
+    int one = 1;
+    openblas_set_num_threads(1);
+    lu = 1./(*N-1);
+    qtemp = calloc(*N, 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);
+        p2 = p[i] - p1;
+        memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
+        pbar = 1-p[i];
+        bound = MIN(M, *N);
+        dscal_(&bound, &pbar, q, &one);
+        bound = MIN(M, *N-d2);
+        daxpy_(&bound, &p1, qtemp, &one, q+d1, &one);
+        daxpy_(&bound, &p2, qtemp, &one, q+d2, &one);
+        M += d2;
+    }
+    /* correction for weight loss */
+    if(M > *N){
+        sum = 0;
+        for(j=0; j<MIN(M, *N); j++){
+            sum += q[j];
+        }
+        q[*N-1] += 1-sum;
+    }
+    free(qtemp);
+}
+
+void lossdistrib_Z(double *p, int *np, double *w, double *S, int *N, int *defaultflag,
+                   double *rho, double *Z, int *nZ, double *q){
+    int i, j;
+    double* pshocked = malloc(sizeof(double) * (*np) * (*nZ));
+
+#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);
+        }
+        lossdistrib_blas(pshocked + (*np) * i, np, w, S + (*np) * i, N,
+                         defaultflag, q + (*N) * i);
+    }
+    free(pshocked);
+}
+
+void lossdistrib_truncated(double *p, int *np, double *w, double *S, int *N,
+                           int *T, int *defaultflag, double *q) {
+    /* recursive algorithm with first order correction for computing
+       the loss distribution.
+       input:
+       p vector of default probabilities
+       np length of p
+       S vector of severities (should be same length as p)
+       N number of ticks in the grid
+       T where to truncate the distribution
+       defaultflag if true computes the default distribution
+       output:
+       q the loss distribution */
+
+    int i, j, d1, d2, M;
+    double lu, d, p1, p2;
+    double *q1, *q2;
+
+    lu = 1./(*N-1);
+    q1 = calloc( *T, sizeof(double));
+    q2 = calloc( *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);
+        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];
+        };
+        M += d2;
+    }
+    free(q1);
+    free(q2);
+}
+
+void lossdistrib_joint(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
+    /* recursive algorithm with first order correction
+       computes jointly the loss and recovery distribution
+       p vector of default probabilities
+       np length of p
+       w vector of issuer weights (length np)
+       S vector of severities (should be same length as p)
+       N number of ticks in the grid
+       defaultflag if true computes the default distribution
+       q the joint probability distribution */
+
+    int i, j, k, m, n;
+    int Mx, My;
+    double lu, x, y, sum;
+    double alpha1, alpha2, beta1, beta2;
+    double w1, w2, w3, w4;
+    double *qtemp;
+
+    lu = 1./(*N-1);
+    qtemp = calloc( (*N) * (*N), sizeof(double));
+    q[0] = 1;
+    Mx=1;
+    My=1;
+    for(k=0; k<(*np); k++){
+        x = (*defaultflag)? w[k] /lu : S[k] * w[k] / lu;
+        y = (1-S[k]) * w[k] / lu;
+        i = floor(x);
+        j = floor(y);
+        alpha1 = i + 1 - x;
+        alpha2 = 1 - alpha1;
+        beta1 = j + 1 - y;
+        beta2 = 1 - beta1;
+        w1 = alpha1 * beta1;
+        w2 = alpha1 * beta2;
+        w3 = alpha2 * beta2;
+        w4 = alpha2 * beta1;
+
+        for(n=0; n<MIN(My, *N); n++){
+            memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+        }
+        for(n=0; n<MIN(My, *N); n++){
+            for(m=0; m<MIN(Mx, *N); m++){
+                q[m+(*N)*n] = (1-p[k])* q[m+(*N)*n];
+            }
+        }
+        for(n=0; n < MIN(My, *N-j-1); n++){
+            for(m=0; m < MIN(Mx, *N-i-1); m++){
+                q[i+m+(*N)*(j+n)] += w1 * p[k] * qtemp[m+(*N)*n];
+                q[i+m+(*N)*(j+1+n)] += w2 * p[k] * qtemp[m+(*N)*n];
+                q[i+1+m+(*N)*(j+1+n)] += w3 * p[k] * qtemp[m+(*N)*n];
+                q[i+1+m+(*N)*(j+n)] += w4 * p[k] *qtemp[m+(*N)*n];
+            }
+        }
+        Mx += i+1;
+        My += j+1;
+    }
+    /* correction for weight loss */
+    if(Mx>*N || My>*N){
+        sum = 0;
+        for(m=0; m < MIN(Mx, *N); m++){
+            for(n=0; n < MIN(My, *N); n++){
+                sum += q[m+n*(*N)];
+            }
+        }
+        q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
+    }
+    free(qtemp);
+}
+
+void lossdistrib_joint_blas(double *p, int *np, double *w, double *S, int *N, int *defaultflag, double *q) {
+    /* recursive algorithm with first order correction
+       computes jointly the loss and recovery distribution
+       p vector of default probabilities
+       np length of p
+       w vector of issuer weights (length np)
+       S vector of severities (should be same length as p)
+       N number of ticks in the grid
+       defaultflag if true computes the default distribution
+       q the joint probability distribution */
+
+    int i, j, k, m, n;
+    int Mx, My;
+    double lu, x, y, sum, pbar;
+    double alpha1, alpha2, beta1, beta2;
+    double w1, w2, w3, w4;
+    double *qtemp;
+    int bound;
+    int one = 1;
+
+    /* only use one thread, performance is horrible if use multiple threads */
+    openblas_set_num_threads(1);
+
+    lu = 1./(*N-1);
+    qtemp = calloc( (*N) * (*N), sizeof(double));
+    q[0] = 1;
+    Mx=1;
+    My=1;
+    for(k=0; k<(*np); k++){
+        x = (*defaultflag)? w[k] /lu : S[k] * w[k] / lu;
+        y = (1-S[k]) * w[k] / lu;
+        i = floor(x);
+        j = floor(y);
+        alpha1 = i + 1 - x;
+        alpha2 = 1 - alpha1;
+        beta1 = j + 1 - y;
+        beta2 = 1 - beta1;
+        w1 = alpha1 * beta1 * p[k];
+        w2 = alpha1 * beta2 * p[k];
+        w3 = alpha2 * beta2 * p[k];
+        w4 = alpha2 * beta1 * p[k];
+
+        for(n=0; n<MIN(My, *N); n++){
+            memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+        }
+
+        bound = MIN(Mx, *N);
+        pbar = 1-p[k];
+        for(n=0; n<MIN(My, *N); n++){
+            dscal_(&bound, &pbar, q+(*N)*n, &one);
+        }
+        bound = MIN(Mx, *N-i-1);
+        for(n=0; n < MIN(My, *N-j-1); n++){
+            daxpy_(&bound, &w1, qtemp+(*N)*n, &one, q+i+(*N)*(j+n), &one);
+            daxpy_(&bound, &w2, qtemp+(*N)*n, &one, q+i+(*N)*(j+1+n), &one);
+            daxpy_(&bound, &w3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j+1+n), &one);
+            daxpy_(&bound, &w4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j+n), &one);
+        }
+        Mx += i+1;
+        My += j+1;
+    }
+    /* correction for weight loss */
+    if(Mx>*N || My>*N){
+        sum = 0;
+        for(m=0; m < MIN(Mx, *N); m++){
+            for(n=0; n < MIN(My, *N); n++){
+                sum += q[m+n*(*N)];
+            }
+        }
+        q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
+    }
+    free(qtemp);
+}
+
+void recovdist(double *dp, double *pp, int *n, double *w, double *S, int *N, double *q) {
+    /* recursive algorithm with first order correction for computing
+       the recovery distribution in case of prepayment.
+       dp vector of default probabilities
+       pp vector of prepay probabilities
+       n length of p
+       S vector of severities (should be same length as p)
+       w vector of weights
+       N number of ticks in the grid
+       q the loss distribution */
+
+    int i, j, d1l, d1u, d2l, d2u;
+    int M;
+    double lu, d1, d2, dp1, dp2, pp1, pp2, sum;
+    double *qtemp;
+
+    lu = 1./(*N - 1);
+    qtemp = calloc( (*N), sizeof(double));
+    q[0] = 1;
+    M=1;
+    for(i=0; i<(*n); i++){
+        d1 = w[i] * (1-S[i]) /lu;
+        d2 = w[i]/lu;
+        d1l = floor(d1);
+        d1u = d1l + 1;
+        d2l = floor(d2);
+        d2u = d2l + 1;
+        dp1 = dp[i] * (d1u - d1);
+        dp2 = dp[i] - dp1;
+        pp1 = pp[i] * (d2u - d2);
+        pp2 = pp[i] - pp1;
+        memcpy(qtemp, q, MIN(M, *N) * sizeof(double));
+        for(j = 0; j< MIN(M, *N); j++){
+            q[j] = (1-dp[i]-pp[i]) * q[j];
+        }
+        for(j=0; j < MIN(M, *N-d2u); j++){
+            q[d1l+j] += dp1 * qtemp[j];
+            q[d1u+j] += dp2 * qtemp[j];
+            q[d2l+j] += pp1 * qtemp[j];
+            q[d2u+j] += pp2 * qtemp[j];
+        };
+        M += d2u;
+    }
+    /* correction for weight loss */
+    if(M>*N){
+        sum = 0;
+        for(j=0; j<MIN(M, *N); j++){
+            sum += q[j];
+        }
+        q[*N-1] += 1-sum;
+    }
+    free(qtemp);
+}
+
+void lossdistrib_prepay_joint(double *dp, double *pp, int *ndp, double *w,
+                              double *S, int *N, int *defaultflag, double *q) {
+    int i, j1, j2, k, m, n;
+    double lu, x, y1, y2, sum;
+    double alpha1, alpha2, beta1, beta2;
+    double dpw1, dpw2, dpw3, dpw4;
+    double ppw1, ppw2, ppw3;
+    double *qtemp;
+    int Mx, My;
+
+    lu = 1./(*N-1);
+    qtemp = calloc((*N) * (*N), sizeof(double));
+    q[0] = 1;
+    Mx=1;
+    My=1;
+
+    for(k=0; k<(*ndp); k++){
+        y1 = (1-S[k]) * w[k]/lu;
+        y2 = w[k]/lu;
+        x = (*defaultflag)? y2: y2-y1;
+        i = floor(x);
+        j1 = floor(y1);
+        j2 = floor(y2);
+        alpha1 = i + 1 - x;
+        alpha2 = 1 - alpha1;
+        beta1 = j1 + 1 - y1;
+        beta2 = 1 - beta1;
+        dpw1 = alpha1 * beta1 * dp[k];
+        dpw2 = alpha1 * beta2 * dp[k];
+        dpw3 = alpha2 * beta2 * dp[k];
+        dpw4 = alpha2 * beta1 * dp[k];
+
+        /* by default distribution, we mean names fractions of names that disappeared
+           either because of default or prepayment */
+        for(n=0; n<MIN(My, *N); n++){
+            memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+        }
+        for(n=0; n<MIN(My, *N); n++){
+            for(m=0; m<MIN(Mx, *N); m++){
+                q[m+(*N)*n] = (1-dp[k]-pp[k]) * q[m+(*N)*n];
+            }
+        }
+        if(*defaultflag){
+            ppw1 = alpha1 * alpha1 * pp[k];
+            ppw2 = alpha1 * alpha2 * pp[k];
+            ppw3 = alpha2 * alpha2 * pp[k];
+            for(n=0; n < MIN(My, *N-j2-1); n++){
+                for(m=0; m < MIN(Mx, *N-i-1); m++){
+                    q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
+                    q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
+                    q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
+                    q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
+
+                    q[i+m+(*N)*(j2+n)] += ppw1 * qtemp[m+(*N)*n];
+                    q[i+m+(*N)*(j2+1+n)] += ppw2 * qtemp[m+(*N)*n];
+                    q[i+1+m+(*N)*(j2+1+n)] += ppw3 * qtemp[m+(*N)*n];
+                    q[i+1+m+(*N)*(j2+n)] += ppw2 * qtemp[m+(*N)*n];
+                }
+            }
+        }else{
+            for(n=0; n < MIN(My, *N-j2-1); n++){
+                for(m=0; m < MIN(Mx, *N-i-1); m++){
+                    q[i+m+(*N)*(j1+n)] += dpw1 * qtemp[m+(*N)*n];
+                    q[i+m+(*N)*(j1+1+n)] += dpw2 * qtemp[m+(*N)*n];
+                    q[i+1+m+(*N)*(j1+1+n)] += dpw3 * qtemp[m+(*N)*n];
+                    q[i+1+m+(*N)*(j1+n)] += dpw4 * qtemp[m+(*N)*n];
+                    q[m+(*N)*(j2+n)] +=  pp[k]*(j2+1-y2) * qtemp[m+(*N)*n];
+                    q[m+(*N)*(j2+1+n)] +=  pp[k]*(y2-j2) * qtemp[m+(*N)*n];
+                }
+            }
+        }
+        Mx += i + 1;
+        My += j2 + 1;
+    }
+    /* correction for weight loss */
+    if(Mx>*N || My>*N){
+        sum = 0;
+        for(m=0; m < MIN(Mx, *N); m++){
+            for(n=0; n < MIN(My, *N); n++){
+                sum += q[m+n*(*N)];
+            }
+        }
+        q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
+    }
+    free(qtemp);
+}
+
+void lossdistrib_prepay_joint_blas(double *dp, double *pp, int *ndp, double *w,
+                                   double *S, int *N, int *defaultflag, double *q) {
+    int i, j1, j2, k, m, n;
+    double lu, x, y1, y2, sum;
+    double alpha1, alpha2, beta1, beta2;
+    double dpw1, dpw2, dpw3, dpw4;
+    double ppw1, ppw2, ppw3;
+    double *qtemp;
+    int Mx, My, bound;
+    double pbar;
+    int one = 1;
+
+    lu = 1./(*N-1);
+    qtemp = calloc((*N) * (*N), sizeof(double));
+    q[0] = 1;
+    Mx=1;
+    My=1;
+
+    /* only use one thread, performance is horrible if use multiple threads */
+    openblas_set_num_threads(1);
+    for(k=0; k<(*ndp); k++){
+        y1 = (1-S[k]) * w[k]/lu;
+        y2 = w[k]/lu;
+        x = (*defaultflag)? y2: y2-y1;
+        i = floor(x);
+        j1 = floor(y1);
+        j2 = floor(y2);
+        alpha1 = i + 1 - x;
+        alpha2 = 1 - alpha1;
+        beta1 = j1 + 1 - y1;
+        beta2 = 1 - beta1;
+        dpw1 = alpha1 * beta1 * dp[k];
+        dpw2 = alpha1 * beta2 * dp[k];
+        dpw3 = alpha2 * beta2 * dp[k];
+        dpw4 = alpha2 * beta1 * dp[k];
+
+        /* by default distribution, we mean names fractions of names that disappeared
+           either because of default or prepayment */
+        for(n=0; n<MIN(My, *N); n++){
+            memcpy(qtemp+n*(*N), q+n*(*N), MIN(Mx, *N) * sizeof(double));
+        }
+        bound = MIN(Mx, *N);
+        pbar = 1-dp[k]-pp[k];
+        for(n=0; n<MIN(My, *N); n++){
+            dscal_(&bound, &pbar, q+(*N)*n, &one);
+        }
+        bound = MIN(Mx, *N-i-1);
+        if(*defaultflag){
+            ppw1 = alpha1 * alpha1 * pp[k];
+            ppw2 = alpha1 * alpha2 * pp[k];
+            ppw3 = alpha2 * alpha2 * pp[k];
+            for(n=0; n < MIN(My, *N-j2-1); n++){
+                daxpy_(&bound, &dpw1, qtemp+(*N)*n, &one, q+i+(*N)*(j1+n), &one);
+                daxpy_(&bound, &dpw2, qtemp+(*N)*n, &one, q+i+(*N)*(j1+1+n), &one);
+                daxpy_(&bound, &dpw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+1+n), &one);
+                daxpy_(&bound, &dpw4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+n), &one);
+
+                daxpy_(&bound, &ppw1, qtemp+(*N)*n, &one, q+i+(*N)*(j2+n), &one);
+                daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+i+(*N)*(j2+1+n), &one);
+                daxpy_(&bound, &ppw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j2+1+n), &one);
+                daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+i+1+(*N)*(j2+n), &one);
+            }
+        }else{
+            ppw1 = pp[k] * (j2+1-y2);
+            ppw2 = pp[k] * (y2-j2);
+            for(n=0; n < MIN(My, *N-j2-1); n++){
+                daxpy_(&bound, &dpw1, qtemp+(*N)*n, &one, q+i+(*N)*(j1+n), &one);
+                daxpy_(&bound, &dpw2, qtemp+(*N)*n, &one, q+i+(*N)*(j1+1+n), &one);
+                daxpy_(&bound, &dpw3, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+1+n), &one);
+                daxpy_(&bound, &dpw4, qtemp+(*N)*n, &one, q+i+1+(*N)*(j1+n), &one);
+                daxpy_(&bound, &ppw1, qtemp+(*N)*n, &one, q+(*N)*(j2+n), &one);
+                daxpy_(&bound, &ppw2, qtemp+(*N)*n, &one, q+(*N)*(j2+1+n), &one);
+            }
+        }
+        Mx += i + 1;
+        My += j2 + 1;
+    }
+    /* correction for weight loss */
+    if(Mx>*N || My>*N){
+        sum = 0;
+        for(m=0; m < MIN(Mx, *N); m++){
+            for(n=0; n < MIN(My, *N); n++){
+                sum += q[m+n*(*N)];
+            }
+        }
+        q[MIN(*N, Mx)*MIN(My,*N)-1] += 1 - sum;
+    }
+    free(qtemp);
+}
+
+double shockprob(double p, double rho, double Z, int give_log){
+    if(rho==1){
+        return((double)(Z<=qnorm(p, 0, 1, 1, 0)));
+    }else{
+        return( pnorm( (qnorm(p, 0, 1, 1, 0) - sqrt(rho) * Z)/sqrt(1 - rho), 0, 1, 1, give_log));
+    }
+}
+
+double dqnorm(double x){
+    return 1/dnorm(qnorm(x, 0, 1, 1, 0), 0, 1, 0);
+}
+
+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){
+    /* 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++){
+        q[i] = shockprob(p, rho, Z[i], 0);
+        q[i + nZ] = dshockprob(p, rho, Z[i]);
+    }
+}
+
+double shockseverity(double S, double Z, double rho, double p){
+    if(p==0){
+        return 0;
+    }else{
+        return( exp(shockprob(S * p, rho, Z, 1) - shockprob(p, rho, Z, 1)) );
+    }
+}
+
+double quantile(double* Z, double* w, int nZ, double p0){
+    double cumw;
+    int i;
+    cumw = 0;
+    for(i=0; i<nZ; i++){
+        cumw += w[i];
+        if(cumw > p0){
+            break;
+        }
+    }
+    return( Z[i] );
+}
+
+void fitprob(double* Z, double* w, int* nZ, double* rho, double* p0, double* result){
+    double eps = 1e-12;
+    int one = 1;
+    double *q = malloc( 2 * (*nZ) * sizeof(double));
+    double dp, p, phi;
+
+    if(*p0 == 0){
+        *result = 0;
+    }else if(*rho == 1){
+        *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);
+        p = *p0;
+        while(fabs(dp) > eps){
+            phi = 1;
+            while( (p - phi * dp) < 0 || (p - phi * dp) > 1){
+                phi *= 0.8;
+            }
+            p -= phi * dp;
+            shockprobvec2(p, *rho, Z, *nZ, q);
+            dp = (ddot_(nZ, q, &one, w, &one) - *p0)/ddot_(nZ, q + *nZ, &one, w, &one);
+        }
+        *result = p;
+    }
+    free(q);
+}
+
+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++){
+            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++){
+            q[i] = fabs(1 - (1 - *Rtilde) * exp( shockprob(ptilde, *rho, Z[i], 1) -
+                                                 shockprob(*pmod, *rho, Z[i], 1)));
+        }
+    }
+}
+
+void lossdistrib_prepay_joint_Z(double *dp, double *pp, int *ndp, double *w,
+                                double *S, int *N, int *defaultflag, double *rho,
+                                double *Z, double *wZ, int *nZ, double *q) {
+    int i, j;
+    double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
+    double* ppshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
+    int N2 = (*N) * (*N);
+    double* qmat = malloc(sizeof(double) * N2 * (*nZ));
+
+    double alpha = 1;
+    double beta = 0;
+    int one = 1;
+
+    #pragma omp parallel for private(j)
+    for(i = 0; i < *nZ; i++){
+        for(j = 0; j < *ndp; j++){
+            dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
+            ppshocked[j + (*ndp) * i] = shockprob(pp[j], rho[j], -Z[i], 0);
+        }
+        lossdistrib_prepay_joint_blas(dpshocked + (*ndp) * i, ppshocked + (*ndp) * i, ndp,
+                                      w, S + (*ndp) * i, N, defaultflag, qmat + N2 * i);
+    }
+
+    dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);
+
+    free(dpshocked);
+    free(ppshocked);
+    free(qmat);
+}
+
+void lossdistrib_joint_Z(double *dp, int *ndp, double *w,
+                         double *S, int *N, int *defaultflag, double *rho,
+                         double *Z, double *wZ, int *nZ, double *q) {
+    int i, j;
+    double* dpshocked = malloc(sizeof(double) * (*ndp) * (*nZ));
+    int N2 = (*N) * (*N);
+    double* qmat = malloc(sizeof(double) * N2 * (*nZ));
+
+    double alpha = 1;
+    double beta = 0;
+    int one = 1;
+
+#pragma omp parallel for private(j)
+    for(i = 0; i < *nZ; i++){
+        for(j = 0; j < *ndp; j++){
+            dpshocked[j + (*ndp) * i] = shockprob(dp[j], rho[j], Z[i], 0);
+        }
+        lossdistrib_joint_blas(dpshocked + (*ndp) * i, ndp, w, S + (*ndp) * i, N,
+                               defaultflag, qmat + N2 * i);
+    }
+
+    dgemv_("n", &N2, nZ, &alpha, qmat, &N2, wZ, &one, &beta, q, &one);
+
+    free(dpshocked);
+    free(qmat);
+}
+
+void BCloss_recov_dist(double *defaultprob, int *dim1, int *dim2,
+                       double *issuerweights, double *recov, double *Z, double *w,
+                       int *n, double *rho, int *N, int *defaultflag,
+                       double *L, double *R) {
+    /*
+      computes the loss and recovery distribution over time with a flat gaussian
+      correlation
+      inputs:
+      defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers
+      and dim2 number of time steps
+      issuerweights: vector of issuer weights (length dim1)
+      recov: vector of recoveries (length dim1)
+      Z: vector of factor values (length n)
+      w: vector of factor weights (length n)
+      rho: correlation beta vector (length dim1)
+      N: number of ticks in the grid
+      defaultflag: if true, computes the default distribution
+      outputs:
+      L: matrix of size (N, dim2)
+      R: matrix of size (N, dim2)
+    */
+    int t, i, j;
+    double g;
+    double *gshocked, *Rshocked, *Sshocked, *Lw, *Rw;
+    int one = 1;
+    double alpha = 1;
+    double beta = 0;
+
+    gshocked = malloc((*dim1) * (*n) * sizeof(double));
+    Rshocked = malloc((*dim1) * (*n) * sizeof(double));
+    Sshocked = malloc((*dim1) * (*n) * sizeof(double));
+    Lw = malloc((*N) * (*n) * sizeof(double));
+    Rw = malloc((*N) * (*n) * sizeof(double));
+
+
+    for(t=0; t < (*dim2); t++) {
+        memset(Lw, 0, (*N) * (*n) * sizeof(double));
+        memset(Rw, 0, (*N) * (*n) * sizeof(double));
+        #pragma omp parallel for private(j, g)
+        for(i=0; i < *n; i++){
+            for(j=0; j < (*dim1); j++){
+                g = defaultprob[j + (*dim1) * t];
+                gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
+                Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
+                Rshocked[j+(*dim1)*i] = 1 - Sshocked[j+(*dim1)*i];
+            }
+            lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
+                             Lw + i * (*N));
+            lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Rshocked + (*dim1)*i, N, defaultflag,
+                             Rw + i * (*N));
+        }
+        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);
+    }
+    free(gshocked);
+    free(Rshocked);
+    free(Sshocked);
+    free(Lw);
+    free(Rw);
+}
+
+
+void BCloss_dist(double *defaultprob, int *dim1, int *dim2,
+                 double *issuerweights, double *recov, double *Z, double *w,
+                 int *n, double *rho, int *N, int *defaultflag,
+                 double *L) {
+    /*
+      computes the loss distribution over time with a flat gaussian
+      correlation
+      inputs:
+      defaultprob: matrix of size dim1 x dim2. dim1 is the number of issuers
+      and dim2 number of time steps
+      issuerweights: vector of issuer weights (length dim1)
+      recov: vector of recoveries (length dim1)
+      Z: vector of factor values (length n)
+      w: vector of factor weights (length n)
+      rho: correlation beta vector (length dim1)
+      N: number of ticks in the grid
+      defaultflag: if true, computes the default distribution
+      outputs:
+      L: matrix of size (N, dim2)
+    */
+    int t, i, j;
+    double g;
+    double *gshocked, *Sshocked, *Lw;
+    int one = 1;
+    double alpha = 1;
+    double beta = 0;
+
+    gshocked = malloc((*dim1) * (*n) * sizeof(double));
+    Sshocked = malloc((*dim1) * (*n) * sizeof(double));
+    Lw = malloc((*N) * (*n) * sizeof(double));
+
+    for(t=0; t < (*dim2); t++) {
+        memset(Lw, 0, (*N) * (*n) * sizeof(double));
+        #pragma omp parallel for private(j, g)
+        for(i=0; i < *n; i++){
+            for(j=0; j < (*dim1); j++){
+                g = defaultprob[j + (*dim1) * t];
+                gshocked[j+(*dim1)*i] = shockprob(g, rho[j], Z[i], 0);
+                Sshocked[j+(*dim1)*i] = shockseverity(1-recov[j], Z[i], rho[j], g);
+            }
+            lossdistrib_blas(gshocked + (*dim1) * i, dim1, issuerweights, Sshocked + (*dim1)*i, N, defaultflag,
+                             Lw + i * (*N));
+        }
+        dgemv_("n", N, n, &alpha, Lw, N, w, &one, &beta, L + t * (*N), &one);
+    }
+    free(gshocked);
+    free(Sshocked);
+    free(Lw);
+}
-- 
cgit v1.2.3-70-g09d2