From 6eac91bf7bcf3515e661d74b49e8179938bd1488 Mon Sep 17 00:00:00 2001
From: Guillaume Horel <guillaume.horel@serenitascapital.com>
Date: Wed, 7 May 2014 10:50:13 -0400
Subject: add GHquad function, get rid of statmod dependency

---
 src/Makevars      |  2 +-
 src/lossdistrib.c | 30 ++++++++++++++++++++++++++++++
 src/lossdistrib.h |  2 ++
 3 files changed, 33 insertions(+), 1 deletion(-)

(limited to 'src')

diff --git a/src/Makevars b/src/Makevars
index f6790ff..0ef4741 100644
--- a/src/Makevars
+++ b/src/Makevars
@@ -1,2 +1,2 @@
 PKG_CFLAGS=$(SHLIB_OPENMP_CFLAGS)
-PKG_LIBS=$(SHLIB_OPENMP_CFLAGS)
+PKG_LIBS=$(SHLIB_OPENMP_CFLAGS) -llapack
diff --git a/src/lossdistrib.c b/src/lossdistrib.c
index 77ea335..f12ab2f 100644
--- a/src/lossdistrib.c
+++ b/src/lossdistrib.c
@@ -11,8 +11,38 @@ extern int dgemv_(char* trans, int *m, int *n, double* alpha, double* A, int* ld
 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
diff --git a/src/lossdistrib.h b/src/lossdistrib.h
index d4bd974..785b4e8 100644
--- a/src/lossdistrib.h
+++ b/src/lossdistrib.h
@@ -44,3 +44,5 @@ void BClossdist(double *SurvProb, int *dim1, int *dim2, double *issuerweights,
                 int *defaultflag, double *L, double *R);
 
 double quantile(double* Z, double* w, int nZ, double p0);
+
+void GHquad(int *n, double* Z, double* w);
-- 
cgit v1.2.3-70-g09d2