diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Makevars | 2 | ||||
| -rw-r--r-- | src/lossdistrib.c | 30 | ||||
| -rw-r--r-- | src/lossdistrib.h | 2 |
3 files changed, 33 insertions, 1 deletions
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);
|