fork download

copy

/* to get the segfault compile with 
 
   gcc ex.c -lgsl -lgscblas -O2
 
*/   
 
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <signal.h>
#include <math.h>
 
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_eigen.h>
 
#define N 16
 
#define NTOT 76 
#define VAL 30
 
#define SC_EIG 2
 
#define EPS 0.000001
 
#define G6LEN(n)  (((n)*((n)-1)/2+5)/6+1) 
#define BIAS6 63
#define TOPBIT6 32
 
static gsl_matrix *adj;
static gsl_matrix *partitioned_adj;
 
static gsl_matrix *adj_ext;
 
static gsl_vector_complex *eval_ext;
static gsl_matrix_complex *evec_ext;
static gsl_eigen_nonsymmv_workspace *w_ext;
 
static char gcode[G6LEN(N)+3]; 
 
const int eigenvalues[NTOT] = {30, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 
    2, 2, 2, 2, 2, 2, 2, 2, 2, 2, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, -8, 
    -8, -8, -8, -8};
 
 
char *adjtog6(gsl_matrix *adj, unsigned n) {
 
    unsigned i,j;
    int  k;
 
    char *p,x;
 
    p = gcode;
    *p++ = BIAS6+n;
 
    k = 6;
    x = 0;
 
    for (j = 1; j < n; ++j) {
        for (i = 0; i < j; ++i) {
            x <<= 1;
            if (gsl_matrix_get(adj, j,i) == 1) x |= 1;
 
            if (--k == 0) {
                *p++ = BIAS6 + x;
                k = 6;
                x = 0;
            }
        }
    }
 
    if (k != 6) *p++ = BIAS6 + (x << k);
 
    *p++ = '\n';
    *p = '\0';
 
    return gcode;
}
 
static gsl_matrix *partitioned_am(gsl_matrix *adj, gsl_matrix *partitioned_adj, unsigned n) {
 
    unsigned i,j;
 
    unsigned total_edges = 0; 
 
    gsl_matrix_set_zero(partitioned_adj);
 
    for (i= 0; i < n; i++) {
        unsigned deg = 0;
        for (j = 0; j < n; j++) {
            unsigned el = gsl_matrix_get(adj,i,j);
            if (el)  {
                gsl_matrix_set(partitioned_adj, i,j,1);
                deg+=1;
            }   
        }
        gsl_matrix_set(partitioned_adj,i, n, VAL-deg);
        gsl_matrix_set(partitioned_adj, n,i, (double)(VAL-deg)/(NTOT-n));
        total_edges += deg;
    }
    gsl_matrix_set(partitioned_adj,n, n, (double) 2*(NTOT*VAL/2 + total_edges/2 - n*VAL)/(NTOT-n));
 
    return partitioned_adj;
}
 
static int cmp(const void* elem1, const void* elem2) {
 
    const double *a = elem1, *b = elem2; 
 
    return *a > *b ? -1 : *a < *b ? 1 : 0;
}
 
 
static double *spectrum(gsl_matrix *m, gsl_vector_complex *eval, gsl_matrix_complex *evec, 
        gsl_eigen_nonsymmv_workspace *w, unsigned size) {
 
    static double eigs[N+2];
    unsigned i;
 
    gsl_eigen_nonsymmv (m, eval, evec, w);    
 
    for(i = 0; i < size ; i++) {
        gsl_complex eval_i = gsl_vector_complex_get (eval, i);
        eigs[i] = GSL_REAL(eval_i);
    }
 
    qsort(eigs,size, sizeof(double),cmp); 
 
    return eigs;
}
 
static unsigned does_interlace(double eigs_sc[N+2], unsigned n) {
 
    unsigned i;
 
    for (i = 0; i < n; i++) {
 
        double expr = eigenvalues[i]-eigs_sc[i];
 
        if (expr <= EPS && fabs(expr) >= EPS) { 
            return 0;    
        }
 
        expr = eigs_sc[i] - eigenvalues[NTOT-n+i];
 
        if (expr <= EPS && fabs(expr) >= EPS) {
            return 0;
        } 
    }
 
    return 1; 
}
 
 
static unsigned isInterlaced(gsl_matrix *adj) {
 
    gsl_matrix *m = partitioned_am(adj,partitioned_adj,N+1);
 
    double *eigs = spectrum(m, eval_ext, evec_ext, w_ext, N+2);
 
    return does_interlace(eigs,N+2);
}
 
static void expand(gsl_matrix *adj, unsigned n) {
 
    unsigned i,j;
 
    gsl_matrix *adj_ext = gsl_matrix_alloc(n+1,n+1);
    gsl_matrix_set_zero(adj_ext);
 
    for (i = 0; i < n; i++) {
        for (j = i+1; j < n; j++) {
            double val = gsl_matrix_get(adj,i,j);
            if (val) {
                gsl_matrix_set(adj_ext, i,j, 1.0);
                gsl_matrix_set(adj_ext, j,i, 1.0);
            }
        }
    }
 
    for (i = 0; i < (1U <<n); i++) {
 
        for (j = 0; j < n; j++) {
            if (i & (1U << j) ) {
                gsl_matrix_set(adj_ext, n, j,1);
                gsl_matrix_set(adj_ext, j, n,1);
            }
        }
 
        if (isInterlaced(adj_ext)) {
            printf("%s", adjtog6(adj_ext, n+1));
        }
 
        for (j = 0; j < n; j++) {
            if (i & (1 << j) ) {
                gsl_matrix_set(adj_ext, n, j,0);
                gsl_matrix_set(adj_ext, j, n,0);
            }
        }
 
    }
    gsl_matrix_free(adj_ext);
}
 
static void stringtoadj(char *s) {
 
	char *p;
	int i,j,k,x = 0;
 
    gsl_matrix_set_zero(adj);
 
    p = s + 1;
 
    k = 1;
    for (j = 1; j < N; ++j) {
        for (i = 0; i < j; ++i) {
            if (--k == 0) {
		        k = 6;
		        x = *(p++) - BIAS6;
            }
 
            if (x & TOPBIT6) {
                gsl_matrix_set(adj,i,j,1);
                gsl_matrix_set(adj,j,i,1);
            }
            x <<= 1;
        }
    }
}
 
int main(int argc, char *argv[]) {
 
    adj = gsl_matrix_calloc(N, N);
    partitioned_adj = gsl_matrix_calloc(N+2, N+2);
 
    adj_ext = gsl_matrix_alloc(N+1, N+1);
 
    eval_ext = gsl_vector_complex_alloc(N+2);
    evec_ext = gsl_matrix_complex_alloc (N+2, N+2);
    w_ext = gsl_eigen_nonsymmv_alloc(N+2);
 
    if (!w_ext || !evec_ext || !eval_ext ||!adj_ext) {
        puts("Errrr");
        return 1;
    }    
 
    if (!adj || !partitioned_adj) { 
        puts("Some error somewhere");
        return 1;
    }
 
    stringtoadj("O??O@BI@?Ej@hBAt}YFwa");
    expand(adj,N);
 
    return 0;
 
}

/home/SaNhCt/ccDWWfFn.o: In function `adjtog6':
prog.c:(.text+0xa5): undefined reference to `gsl_matrix_get'
/home/SaNhCt/ccDWWfFn.o: In function `main':
prog.c:(.text.startup+0x19): undefined reference to `gsl_matrix_calloc'
prog.c:(.text.startup+0x29): undefined reference to `gsl_matrix_calloc'
prog.c:(.text.startup+0x39): undefined reference to `gsl_matrix_alloc'
prog.c:(.text.startup+0x4a): undefined reference to `gsl_vector_complex_alloc'
prog.c:(.text.startup+0x5a): undefined reference to `gsl_matrix_complex_alloc'
prog.c:(.text.startup+0x6b): undefined reference to `gsl_eigen_nonsymmv_alloc'
prog.c:(.text.startup+0xce): undefined reference to `gsl_matrix_set_zero'
prog.c:(.text.startup+0x125): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x141): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x16a): undefined reference to `gsl_matrix_alloc'
prog.c:(.text.startup+0x174): undefined reference to `gsl_matrix_set_zero'
prog.c:(.text.startup+0x19f): undefined reference to `gsl_matrix_get'
prog.c:(.text.startup+0x1c4): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x1da): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x225): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x23a): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x255): undefined reference to `gsl_matrix_set_zero'
prog.c:(.text.startup+0x287): undefined reference to `gsl_matrix_get'
prog.c:(.text.startup+0x2be): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x2f5): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x312): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x353): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x371): undefined reference to `gsl_eigen_nonsymmv'
prog.c:(.text.startup+0x38a): undefined reference to `gsl_vector_complex_get'
prog.c:(.text.startup+0x485): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x49a): undefined reference to `gsl_matrix_set'
prog.c:(.text.startup+0x4bd): undefined reference to `gsl_matrix_free'
collect2: error: ld returned 1 exit status

Standard output is empty