diff options
| author | Guillaume Horel <guillaume.horel@serenitascapital.com> | 2014-05-07 10:50:13 -0400 |
|---|---|---|
| committer | Guillaume Horel <guillaume.horel@serenitascapital.com> | 2014-05-07 10:50:13 -0400 |
| commit | 6eac91bf7bcf3515e661d74b49e8179938bd1488 (patch) | |
| tree | 76a7ff22c3d6f5be749c2efd9fb01138447c384b /src/lossdistrib.c | |
| parent | 216fe79edf63926a28a4a435b78ebf03fc96b1d1 (diff) | |
| download | lossdistrib-6eac91bf7bcf3515e661d74b49e8179938bd1488.tar.gz | |
add GHquad function, get rid of statmod dependency
Diffstat (limited to 'src/lossdistrib.c')
| -rw-r--r-- | src/lossdistrib.c | 30 |
1 files changed, 30 insertions, 0 deletions
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
|