diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/break_edge.c | 4 | ||||
| -rw-r--r-- | src/corr_test.c | 4 | ||||
| -rw-r--r-- | src/factor_update.c | 38 | ||||
| -rw-r--r-- | src/fracture.c | 4 | ||||
| -rw-r--r-- | src/fracture.h | 21 | ||||
| -rw-r--r-- | src/get_current.c | 98 | ||||
| -rw-r--r-- | src/net.c | 130 | ||||
| -rw-r--r-- | src/net_conductivity.c (renamed from src/get_conductivity.c) | 2 | ||||
| -rw-r--r-- | src/net_copy.c | 35 | ||||
| -rw-r--r-- | src/net_create.c | 69 | ||||
| -rw-r--r-- | src/net_currents.c | 35 | ||||
| -rw-r--r-- | src/net_fracture.c | 6 | ||||
| -rw-r--r-- | src/net_free.c | 14 | ||||
| -rw-r--r-- | src/net_notch.c | 29 | ||||
| -rw-r--r-- | src/net_voltages.c | 22 | ||||
| -rw-r--r-- | src/rand.c | 20 | ||||
| -rw-r--r-- | src/update_factor.c | 44 |
17 files changed, 266 insertions, 309 deletions
diff --git a/src/break_edge.c b/src/break_edge.c index daed1c5..54eaf34 100644 --- a/src/break_edge.c +++ b/src/break_edge.c @@ -7,7 +7,7 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) { unsigned int v1 = instance->graph->ev_break[2 * edge]; unsigned int v2 = instance->graph->ev_break[2 * edge + 1]; - if (instance->factor != NULL) update_factor(instance->factor, v1, v2, c); + if (instance->factor != NULL) factor_update(instance->factor, v1, v2, c); if (instance->graph->boundary != TORUS_BOUND) { unsigned int w1 = instance->graph->ev[2 * edge]; @@ -93,7 +93,7 @@ bool break_edge(net_t *instance, unsigned int edge, cholmod_common *c) { lap_x[lap_p[tw2] + i] = 1; } - set_connected(instance->adjacency, instance->dual_marks, dw1, instance->dual_marks[v1], -1, 0); + set_connected(instance->adjacency, instance->dual_marks, dw1, instance->dual_marks[dw2], -1, 0); } diff --git a/src/corr_test.c b/src/corr_test.c index 4805fb0..3d56c07 100644 --- a/src/corr_test.c +++ b/src/corr_test.c @@ -7,8 +7,8 @@ int main() { unsigned int width = 64; - graph_t *network = ini_square_network(width, true, false, &c); - net_t *instance = net_create(network, 1e-14, 0.001, 0, true, &c); + graph_t *network = graph_create(VORONOI_LATTICE, TORUS_BOUND, 64, false, &c); + net_t *instance = net_create(network, 1e-14, 0.5, 0, true, &c); net_fracture(instance, &c, 1e-10); double *corr = get_corr(instance, NULL, &c); diff --git a/src/factor_update.c b/src/factor_update.c new file mode 100644 index 0000000..f27971b --- /dev/null +++ b/src/factor_update.c @@ -0,0 +1,38 @@ + +#include "fracture.h" + +void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c) { + uint_t n = factor->n; + + cholmod_sparse *update_mat = + CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); + + uint_t s1, s2; + s1 = v1 < v2 ? v1 : v2; + s2 = v1 > v2 ? v1 : v2; + + int_t *pp = (int_t *)update_mat->p; + int_t *ii = (int_t *)update_mat->i; + double *xx = (double *)update_mat->x; + + for (uint_t i = 0; i <= s1; i++) { + pp[i] = 0; + } + + for (uint_t i = s1 + 1; i <= n; i++) { + pp[i] = 2; + } + + ii[0] = s1; + ii[1] = s2; + xx[0] = 1; + xx[1] = -1; + + cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( + update_mat, factor->Perm, factor->n, NULL, -1, true, true, c); + + CHOL_F(updown)(false, perm_update_mat, factor, c); + + CHOL_F(free_sparse)(&perm_update_mat, c); + CHOL_F(free_sparse)(&update_mat, c); +} diff --git a/src/fracture.c b/src/fracture.c index f36cfdb..f38568f 100644 --- a/src/fracture.c +++ b/src/fracture.c @@ -320,7 +320,7 @@ int main(int argc, char *argv[]) { } if (save_stress_field || save_voltage_field) { - double *tmp_voltages = get_voltage(net, &c); + double *tmp_voltages = net_voltages(net, &c); if (save_voltage_field) { for (uint_t j = 0; j < g->nv; j++) { voltage_field[3 * voltage_pos] = g->vx[2 * j]; @@ -330,7 +330,7 @@ int main(int argc, char *argv[]) { } } if (save_stress_field) { - double *tmp_currents = get_current_v(net, tmp_voltages, &c); + double *tmp_currents = net_currents(net, tmp_voltages, &c); for (uint_t j = 0; j < g->ne; j++) { stress_field[3 * stress_pos] = g->ex[2 * j]; stress_field[3 * stress_pos + 1] = g->ex[2 * j + 1]; diff --git a/src/fracture.h b/src/fracture.h index 2401d25..5e1b3e4 100644 --- a/src/fracture.h +++ b/src/fracture.h @@ -29,6 +29,8 @@ #define CINT_MAX INT_MAX #define CHOL_F(x) cholmod_##x +#define GSL_RAND_GEN gsl_rng_mt19937 + typedef enum lattice_t { VORONOI_LATTICE, SQUARE_LATTICE @@ -111,16 +113,13 @@ int dual_edge_to_verts(unsigned int width, bool periodic, unsigned int edge, double dual_vert_to_coord(unsigned int width, bool periodic, unsigned int vert, bool index); -bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2, - cholmod_common *c); - -data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff); - -double *get_current(const net_t *instance, cholmod_common *c); -double *get_current_v(const net_t *instance, double *voltages, cholmod_common *c); -double *get_voltage(const net_t *instance, cholmod_common *c); +void factor_update(cholmod_factor *factor, uint_t v1, uint_t v2, cholmod_common *c); void net_notch(net_t *net, double notch_len, cholmod_common *c); +data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff); +double *net_voltages(const net_t *net, cholmod_common *c); +double *net_currents(const net_t *net, const double *voltages, cholmod_common *c); +double net_conductivity(const net_t *net, const double *voltages); void update_boundary(net_t *instance, const double *avg_field); @@ -174,10 +173,12 @@ data_t *data_create(uint_t num_edges); void data_free(data_t *data); void data_update(data_t *data, uint_t last_broke, double strength, double conductivity); -double get_conductivity(net_t *inst, double *current, cholmod_common *c); - graph_t *graph_create(lattice_t lattice, bound_t bound, uint_t L, bool dual, cholmod_common *c); uint_t find_cycles(uint_t num_edges, const bool *fuses, const uint_t *ev, const uint_t *vei, const uint_t *ve, int **cycles); bool set_connected(const cholmod_sparse *laplacian, uint_t *marks, int vertex, int label, int stop_at, int exclude); + +unsigned long int rand_seed(); + +double rand_dist_pow(const gsl_rng *r, double beta); diff --git a/src/get_current.c b/src/get_current.c deleted file mode 100644 index a83c399..0000000 --- a/src/get_current.c +++ /dev/null @@ -1,98 +0,0 @@ - -#include "fracture.h" - -double *get_voltage(const net_t *instance, cholmod_common *c) { - cholmod_dense *b = instance->boundary_cond; - cholmod_factor *factor = instance->factor; - - cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); - - if (((double *)x->x)[0] != ((double *)x->x)[0]) { - for (uint_t i = 0; i < instance->graph->ne; i++) { - printf("%d ", instance->fuses[i]); - } - printf("\n"); - printf("GET_VOLTAGE: value is NaN\n"); - exit(EXIT_FAILURE); - } - - double *field = (double *)x->x; - x->x = NULL; - - CHOL_F(free_dense)(&x, c); - - - return field; -} - -double *get_current(const net_t *instance, cholmod_common *c) { - unsigned int num_edges = instance->graph->ne; - unsigned int num_gverts = instance->graph->break_dim; - cholmod_sparse *voltcurmat = instance->graph->voltcurmat; - - double *voltages = get_voltage(instance, c); - if (voltages == NULL) { - return NULL; - } - cholmod_dense *x = CHOL_F(allocate_dense)( - num_gverts, 1, num_gverts, CHOLMOD_REAL, c); - double *tmp_x = x->x; - x->x = voltages; - - cholmod_dense *y = - CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c); - - double alpha[2] = {1, 0}; - double beta[2] = {0, 0}; - CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); - - double *field = (double *)malloc(num_edges * sizeof(double)); - - for (int i = 0; i < num_edges; i++) { - if (instance->fuses[i]) - field[i] = 0; - else - field[i] = ((double *)y->x)[i]; - } - - x->x = tmp_x; - free(voltages); - CHOL_F(free_dense)(&x, c); - CHOL_F(free_dense)(&y, c); - - return field; -} - - -double *get_current_v(const net_t *instance, double *voltages, cholmod_common *c) { - unsigned int num_edges = instance->graph->ne; - unsigned int num_gverts = instance->graph->break_dim; - cholmod_sparse *voltcurmat = instance->graph->voltcurmat; - - cholmod_dense *x = CHOL_F(allocate_dense)( - num_gverts, 1, num_gverts, CHOLMOD_REAL, c); - double *tmp_x = x->x; - x->x = voltages; - - cholmod_dense *y = - CHOL_F(allocate_dense)(num_edges, 1, num_edges, CHOLMOD_REAL, c); - - double alpha[2] = {1, 0}; - double beta[2] = {0, 0}; - CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); - - double *field = (double *)malloc(num_edges * sizeof(double)); - - for (int i = 0; i < num_edges; i++) { - if (instance->fuses[i]) - field[i] = 0; - else - field[i] = ((double *)y->x)[i]; - } - - x->x = tmp_x; - CHOL_F(free_dense)(&x, c); - CHOL_F(free_dense)(&y, c); - - return field; -} diff --git a/src/net.c b/src/net.c new file mode 100644 index 0000000..9f6965b --- /dev/null +++ b/src/net.c @@ -0,0 +1,130 @@ + +#include "fracture.h" + +double *get_thres(uint_t ne, double beta) { + double *thres = (double *)malloc(ne * sizeof(double)); + assert(thres != NULL); + + gsl_rng *r = gsl_rng_alloc(GSL_RAND_GEN); + gsl_rng_set(r, rand_seed()); + + for (uint_t i = 0; i < ne; i++) { + thres[i] = rand_dist_pow(r, beta); + } + + gsl_rng_free(r); + + return thres; +} + +void net_notch(net_t *net, double notch_len, cholmod_common *c) { + for (uint_t i = 0; i < net->graph->ne; i++) { + uint_t v1, v2; + double v1x, v1y, v2x, v2y, dy; + bool crosses_center, not_wrapping, correct_length; + + v1 = net->graph->ev[2 * i]; + v2 = net->graph->ev[2 * i + 1]; + + v1x = net->graph->vx[2 * v1]; + v1y = net->graph->vx[2 * v1 + 1]; + v2x = net->graph->vx[2 * v2]; + v2y = net->graph->vx[2 * v2 + 1]; + + dy = v1y - v2y; + + crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5); + not_wrapping = fabs(dy) < 0.5; + //correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len; + correct_length = v1x < notch_len && v2x < notch_len; + + if (crosses_center && not_wrapping && correct_length) { + break_edge(net, i, c); + } + } +} + +net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) { + net_t *net = (net_t *)calloc(1, sizeof(net_t)); + assert(net != NULL); + + net->graph = g; + net->fuses = (bool *)calloc(g->ne, sizeof(bool)); + assert(net->fuses != NULL); + net->thres = get_thres(g->ne, beta); + net->inf = inf; + + net->voltage_bound = vb; + net->boundary_cond = bound_set(g, vb, notch_len, c); + + if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c); + else net->adjacency = gen_adjacency(net, true, false, 0, c); + + net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t)); + assert(net->marks != NULL); + + net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t)); + assert(net->dual_marks != NULL); + + for (uint_t i = 0; i < (net->graph->break_dim); i++) { + net->marks[i] = 1; + } + for (uint_t i = 0; i < (net->graph->dnv); i++) { + net->dual_marks[i] = i+1; + } + net->num_components = 1; + + net_notch(net, notch_len, c); + + { + cholmod_sparse *laplacian = gen_laplacian(net, c, true); + net->factor = CHOL_F(analyze)(laplacian, c); + CHOL_F(factorize)(laplacian, net->factor, c); + CHOL_F(free_sparse)(&laplacian, c); + } + + return net; +} + +net_t *net_copy(const net_t *net, cholmod_common *c) { + net_t *net_copy = (net_t *)calloc(1, sizeof(net_t)); + assert(net_copy != NULL); + memcpy(net_copy, net, sizeof(net_t)); + + size_t fuses_size = (net->graph)->ne * sizeof(bool); + net_copy->fuses = (bool *)malloc(fuses_size); + assert(net_copy->fuses != NULL); + memcpy(net_copy->fuses, net->fuses, fuses_size); + + size_t thres_size = (net->graph)->ne * sizeof(double); + net_copy->thres = (double *)malloc(thres_size); + assert(net_copy->thres != NULL); + memcpy(net_copy->thres, net->thres, thres_size); + + size_t marks_size = (net->graph->break_dim) * sizeof(uint_t); + net_copy->marks = (uint_t *)malloc(marks_size); + assert(net_copy->marks != NULL); + memcpy(net_copy->marks, net->marks, marks_size); + + size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t); + net_copy->dual_marks = (uint_t *)malloc(dual_marks_size); + assert(net_copy->dual_marks != NULL); + memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size); + + net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c); + net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c); + net_copy->factor = CHOL_F(copy_factor)(net->factor, c); + + return net_copy; +} + +void net_free(net_t *net, cholmod_common *c) { + free(net->fuses); + free(net->thres); + CHOL_F(free_dense)(&(net->boundary_cond), c); + CHOL_F(free_sparse)(&(net->adjacency), c); + CHOL_F(free_factor)(&(net->factor), c); + free(net->marks); + free(net->dual_marks); + free(net); +} diff --git a/src/get_conductivity.c b/src/net_conductivity.c index 23b7056..e9325bb 100644 --- a/src/get_conductivity.c +++ b/src/net_conductivity.c @@ -1,7 +1,7 @@ #include "fracture.h" -double get_conductivity(net_t *net, double *voltages, cholmod_common *c) { +double net_conductivity(const net_t *net, const double *voltages) { if (net->voltage_bound) { // the voltage drop across the network is fixed to one with voltage // boundary conditions, so the conductivity is the total current flowing diff --git a/src/net_copy.c b/src/net_copy.c deleted file mode 100644 index dcb4080..0000000 --- a/src/net_copy.c +++ /dev/null @@ -1,35 +0,0 @@ - -#include "fracture.h" - -net_t *net_copy(const net_t *net, cholmod_common *c) { - net_t *net_copy = (net_t *)calloc(1, sizeof(net_t)); - assert(net_copy != NULL); - memcpy(net_copy, net, sizeof(net_t)); - - size_t fuses_size = (net->graph)->ne * sizeof(bool); - net_copy->fuses = (bool *)malloc(fuses_size); - assert(net_copy->fuses != NULL); - memcpy(net_copy->fuses, net->fuses, fuses_size); - - size_t thres_size = (net->graph)->ne * sizeof(double); - net_copy->thres = (double *)malloc(thres_size); - assert(net_copy->thres != NULL); - memcpy(net_copy->thres, net->thres, thres_size); - - size_t marks_size = (net->graph->break_dim) * sizeof(uint_t); - net_copy->marks = (uint_t *)malloc(marks_size); - assert(net_copy->marks != NULL); - memcpy(net_copy->marks, net->marks, marks_size); - - size_t dual_marks_size = (net->graph->dnv) * sizeof(uint_t); - net_copy->dual_marks = (uint_t *)malloc(dual_marks_size); - assert(net_copy->dual_marks != NULL); - memcpy(net_copy->dual_marks, net->dual_marks, dual_marks_size); - - net_copy->adjacency = CHOL_F(copy_sparse)(net->adjacency, c); - net_copy->boundary_cond = CHOL_F(copy_dense)(net->boundary_cond, c); - net_copy->factor = CHOL_F(copy_factor)(net->factor, c); - - return net_copy; -} - diff --git a/src/net_create.c b/src/net_create.c deleted file mode 100644 index fedcb3a..0000000 --- a/src/net_create.c +++ /dev/null @@ -1,69 +0,0 @@ - -#include "fracture.h" - -double *get_thres(uint_t ne, double beta) { - assert(beta > 0); - - double *thres = (double *)malloc(ne * sizeof(double)); - assert(thres != NULL); - - gsl_set_error_handler_off(); - - gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937); - { - FILE *rf = fopen("/dev/urandom", "r"); - unsigned long int seed; - fread(&seed, sizeof(unsigned long int), 1, rf); - fclose(rf); - gsl_rng_set(r, seed); - } - - for (uint_t i = 0; i < ne; i++) { - while ((thres[i] = exp(log(gsl_ran_flat(r, 0, 1)) / beta)) == 0.0); - } - - gsl_rng_free(r); - - return thres; -} - -net_t *net_create(const graph_t *g, double inf, double beta, double notch_len, bool vb, cholmod_common *c) { - net_t *net = (net_t *)calloc(1, sizeof(net_t)); - assert(net != NULL); - - net->graph = g; - net->fuses = (bool *)calloc(g->ne, sizeof(bool)); - assert(net->fuses != NULL); - net->thres = get_thres(g->ne, beta); - net->inf = inf; - - net->voltage_bound = vb; - net->boundary_cond = bound_set(g, vb, notch_len, c); - - if (g->boundary != TORUS_BOUND) net->adjacency = gen_adjacency(net, false, false, 0, c); - else net->adjacency = gen_adjacency(net, true, false, 0, c); - - net->marks = (uint_t *)malloc((net->graph->break_dim) * sizeof(uint_t)); - net->dual_marks = (uint_t *)malloc((net->graph->dnv) * sizeof(uint_t)); - assert(net->marks != NULL); - - for (uint_t i = 0; i < (net->graph->break_dim); i++) { - net->marks[i] = 1; - } - for (uint_t i = 0; i < (net->graph->dnv); i++) { - net->dual_marks[i] = i+1; - } - net->num_components = 1; - - net_notch(net, notch_len, c); - - { - cholmod_sparse *laplacian = gen_laplacian(net, c, true); - net->factor = CHOL_F(analyze)(laplacian, c); - CHOL_F(factorize)(laplacian, net->factor, c); - CHOL_F(free_sparse)(&laplacian, c); - } - - return net; -} - diff --git a/src/net_currents.c b/src/net_currents.c new file mode 100644 index 0000000..431818f --- /dev/null +++ b/src/net_currents.c @@ -0,0 +1,35 @@ + +#include "fracture.h" + +double *net_currents(const net_t *net, const double *voltages, cholmod_common *c) { + uint_t ne = net->graph->ne; + uint_t dim = net->graph->break_dim; + cholmod_sparse *voltcurmat = net->graph->voltcurmat; + + cholmod_dense *x = CHOL_F(allocate_dense)(dim, 1, dim, CHOLMOD_REAL, c); + + double *tmp_x = x->x; + x->x = (void *)voltages; + + cholmod_dense *y = CHOL_F(allocate_dense)(ne, 1, ne, CHOLMOD_REAL, c); + + double alpha[2] = {1, 0}; + double beta[2] = {0, 0}; + CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + + double *currents = (double *)malloc(ne * sizeof(double)); + + for (int i = 0; i < ne; i++) { + if (net->fuses[i]) { + currents[i] = 0; + } else { + currents[i] = ((double *)y->x)[i]; + } + } + + x->x = tmp_x; + CHOL_F(free_dense)(&x, c); + CHOL_F(free_dense)(&y, c); + + return currents; +} diff --git a/src/net_fracture.c b/src/net_fracture.c index 48e81f9..54f37dd 100644 --- a/src/net_fracture.c +++ b/src/net_fracture.c @@ -31,10 +31,10 @@ data_t *net_fracture(net_t *net, cholmod_common *c, double cutoff) { data_t *data = data_create(net->graph->ne); while (true) { - double *voltages = get_voltage(net, c); - double *currents = get_current_v(net, voltages, c); + double *voltages = net_voltages(net, c); + double *currents = net_currents(net, voltages, c); - double conductivity = get_conductivity(net, voltages, c); + double conductivity = net_conductivity(net, voltages); if (conductivity < 1e-12 && net->graph->boundary == TORUS_BOUND) { free(voltages); diff --git a/src/net_free.c b/src/net_free.c deleted file mode 100644 index 8b5af50..0000000 --- a/src/net_free.c +++ /dev/null @@ -1,14 +0,0 @@ - -#include "fracture.h" - -void net_free(net_t *net, cholmod_common *c) { - free(net->fuses); - free(net->thres); - CHOL_F(free_dense)(&(net->boundary_cond), c); - CHOL_F(free_sparse)(&(net->adjacency), c); - CHOL_F(free_factor)(&(net->factor), c); - free(net->marks); - free(net->dual_marks); - free(net); -} - diff --git a/src/net_notch.c b/src/net_notch.c deleted file mode 100644 index ccbb387..0000000 --- a/src/net_notch.c +++ /dev/null @@ -1,29 +0,0 @@ - -#include "fracture.h" - -void net_notch(net_t *net, double notch_len, cholmod_common *c) { - for (uint_t i = 0; i < net->graph->ne; i++) { - uint_t v1, v2; - double v1x, v1y, v2x, v2y, dy; - bool crosses_center, not_wrapping, correct_length; - - v1 = net->graph->ev[2 * i]; - v2 = net->graph->ev[2 * i + 1]; - - v1x = net->graph->vx[2 * v1]; - v1y = net->graph->vx[2 * v1 + 1]; - v2x = net->graph->vx[2 * v2]; - v2y = net->graph->vx[2 * v2 + 1]; - - dy = v1y - v2y; - - crosses_center = (v1y >= 0.5 && v2y <= 0.5) || (v1y <= 0.5 && v2y >= 0.5); - not_wrapping = fabs(dy) < 0.5; - //correct_length = v1x + dx / dy * (v1y - 0.5) <= notch_len; - correct_length = v1x < notch_len && v2x < notch_len; - - if (crosses_center && not_wrapping && correct_length) { - break_edge(net, i, c); - } - } -} diff --git a/src/net_voltages.c b/src/net_voltages.c new file mode 100644 index 0000000..dedf5b2 --- /dev/null +++ b/src/net_voltages.c @@ -0,0 +1,22 @@ + +#include "fracture.h" + +double *net_voltages(const net_t *net, cholmod_common *c) { + cholmod_dense *b = net->boundary_cond; + cholmod_factor *factor = net->factor; + + cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); + + if (((double *)x->x)[0] != ((double *)x->x)[0]) { + printf("GET_VOLTAGE: value is NaN\n"); + exit(EXIT_FAILURE); + } + + double *voltages = (double *)x->x; + x->x = NULL; + + CHOL_F(free_dense)(&x, c); + + return voltages; +} + diff --git a/src/rand.c b/src/rand.c new file mode 100644 index 0000000..1a6a3d4 --- /dev/null +++ b/src/rand.c @@ -0,0 +1,20 @@ + +#include "fracture.h" + +unsigned long int rand_seed() { + FILE *f = fopen("/dev/urandom", "r"); + unsigned long int seed; + fread(&seed, sizeof(unsigned long int), 1, f); + fclose(f); + return seed; +} + +double rand_dist_pow(const gsl_rng *r, double beta) { + double x = 0; + + // underflow means that for very small beta x is sometimes identically zero, + // which causes problems + while ((x = exp(log(gsl_rng_uniform_pos(r)) / beta)) == 0.0); + + return x; +} diff --git a/src/update_factor.c b/src/update_factor.c deleted file mode 100644 index debfe7b..0000000 --- a/src/update_factor.c +++ /dev/null @@ -1,44 +0,0 @@ - -#include "fracture.h" - -bool update_factor(cholmod_factor *factor, unsigned int v1, unsigned int v2, - cholmod_common *c) { - int n = factor->n; - assert(v1 < n); - assert(v2 < n); - - cholmod_sparse *update_mat = - CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); - - unsigned int v3, v4; - v3 = v1 < v2 ? v1 : v2; - v4 = v1 > v2 ? v1 : v2; - - for (int i = 0; i < n; i++) { - if (i <= v3) - ((int_t *)update_mat->p)[i] = 0; - else - ((int_t *)update_mat->p)[i] = 2; - } - ((int_t *)update_mat->p)[n] = 2; - ((int_t *)update_mat->i)[0] = v3; - ((int_t *)update_mat->i)[1] = v4; - ((double *)update_mat->x)[0] = 1; - ((double *)update_mat->x)[1] = -1; - - // assert(CHOL_F(check_sparse)(update_mat, c)); - - cholmod_sparse *perm_update_mat = CHOL_F(submatrix)( - update_mat, factor->Perm, factor->n, NULL, -1, true, true, c); - - // assert(CHOL_F(check_sparse)(perm_update_mat, c)); - - CHOL_F(updown)(false, perm_update_mat, factor, c); - - // assert(CHOL_F(check_factor)(factor, c)); - - CHOL_F(free_sparse)(&perm_update_mat, c); - CHOL_F(free_sparse)(&update_mat, c); - - return true; -} |