1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
|
#include "stdlib.h"
#include "math.h"
extern int dstev_(char* JOBZ, int* n, double* D, double* E, double* Z, int* ldz, double* WORK, int* INFO);
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);
}
|
