#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);
}