From 1e1fdfc2e3892667bccaf317a01defd8832041c7 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 16 Jan 2017 01:31:10 -0500 Subject: fixed voltage and torus conditions, current and free boundaries and broken right now --- src/gen_laplacian.c | 253 +++++++++++++++------------------------------------- 1 file changed, 72 insertions(+), 181 deletions(-) (limited to 'src/gen_laplacian.c') diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c index 411734c..1cc93a4 100644 --- a/src/gen_laplacian.c +++ b/src/gen_laplacian.c @@ -1,118 +1,79 @@ #include "fracture.h" -cholmod_sparse *gen_adjacency(const net_t *instance, bool dual, bool breakv, - uint_t pad, cholmod_common *c) { - uint_t ne = instance->graph->ne; - if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) { - ne += instance->graph->bound_inds[2]; +cholmod_sparse *gen_adjacency(const net_t *net, bool dual, bool use_gp, bool symmetric, cholmod_common *c) { + const graph_t *g = net->graph; + uint_t nv; + uint_t ne; + uint_t nre; + uint_t *ev; + + if (use_gp) { + nv = net->dim; + ne = net->nep; + nre = (int_t)net->nep - (int_t)net->num_broken; + ev = net->evp; + } else if (dual) { + nv = g->dnv; + ne = g->ne; + nre = net->num_broken; + ev = g->dev; + } else { + nv = g->nv; + ne = g->ne; + nre = (int_t)g->ne - (int_t)net->num_broken; + ev = g->ev; } - uint_t nnz = 2 * ne; + uint_t nnz = nre; - uint_t nv; - if (dual) - nv = instance->graph->dnv; - else { - if (breakv) nv = instance->graph->nv_break; - else nv = instance->graph->nv; - if (instance->graph->boundary != TORUS_BOUND && !breakv) { - nv += 2; - } - } + cholmod_triplet *t = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c); - cholmod_triplet *t = - CHOL_F(allocate_triplet)(nv + pad, nv + pad, nnz, 0, CHOLMOD_REAL, c); int_t *ri = (int_t *)t->i; int_t *ci = (int_t *)t->j; double *ai = (double *)t->x; t->nnz = nnz; - uint_t *etv; - if (dual) - etv = instance->graph->dev; - else { - if (breakv) etv = instance->graph->ev_break; - else etv = instance->graph->ev; - } - - uint_t count = 0; - - for (uint_t i = 0; i < instance->graph->ne; i++) { - if ((instance->fuses[i] && dual) || (!instance->fuses[i] && !dual)) { - uint_t v1 = etv[2 * i]; - uint_t v2 = etv[2 * i + 1]; + uint_t a = 0; - ri[2 * count] = v1; - ci[2 * count] = v2; - ai[2 * count] = 1; + for (uint_t i = 0; i < ne; i++) { + if ((net->fuses[i] && dual) || (!net->fuses[i] && !dual)) { + uint_t v1 = ev[2 * i]; + uint_t v2 = ev[2 * i + 1]; - ri[2 * count + 1] = v2; - ci[2 * count + 1] = v1; - ai[2 * count + 1] = 1; + uint_t s1 = v1 < v2 ? v1 : v2; + uint_t s2 = v1 < v2 ? v2 : v1; - count++; - } else { - uint_t v1 = etv[2 * i]; - uint_t v2 = etv[2 * i + 1]; + ri[a] = s2; + ci[a] = s1; + ai[a] = 1; - ri[2 * count] = v1; - ci[2 * count] = v2; - ai[2 * count] = 0; - - ri[2 * count + 1] = v2; - ci[2 * count + 1] = v1; - ai[2 * count + 1] = 0; - - count++; + a++; } } - if (!breakv && instance->graph->boundary != TORUS_BOUND && !dual) { - for (uint_t i = 0; i < 2; i++) { - for (uint_t j = 0; j < instance->graph->bound_inds[i+1] - instance->graph->bound_inds[i]; j++) { - ri[2*count] = instance->graph->nv + i; - ci[2*count] = instance->graph->bound_verts[instance->graph->bound_inds[i] + j]; - ai[2*count] = 1; - ri[2*count+1] = instance->graph->bound_verts[instance->graph->bound_inds[i] + j]; - ci[2*count+1] = instance->graph->nv + i; - ai[2*count+1] = 1; - count++; - } - } - } - - cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c); CHOL_F(free_triplet)(&t, c); + if (!symmetric) { + cholmod_sparse *tmp_s = CHOL_F(copy)(s, 0, 1, c); + CHOL_F(free_sparse)(&s, c); + s = tmp_s; + } + return s; } -cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c, - bool symmetric) { - const graph_t *network = instance->graph; - uint_t num_verts = network->nv_break; - uint_t num_bounds = network->num_bounds; - double inf = instance->inf; - bool voltage_bound = instance->voltage_bound; - uint_t *bound_inds = network->bound_inds; - uint_t *bound_verts = network->bound_verts; - - uint_t num_gedges; - if (voltage_bound) { - num_gedges = network->bound_inds[num_bounds]; - } else { - num_gedges = network->bound_inds[num_bounds] + (num_bounds - 2) * 2; - } - if (network->boundary == TORUS_BOUND) num_gedges = bound_inds[1]; - uint_t num_gverts = network->break_dim; +cholmod_sparse *gen_laplacian(const net_t *net, cholmod_common *c) { + const graph_t *g = net->graph; + uint_t nv = net->dim; + uint_t ne = net->nep; + uint_t *ev = net->evp; - uint_t nnz = num_gverts + 2 * num_gedges; + uint_t nnz = nv; - cholmod_triplet *temp_m = - CHOL_F(allocate_triplet)(num_gverts, num_gverts, nnz, 0, CHOLMOD_REAL, c); + cholmod_triplet *temp_m = CHOL_F(allocate_triplet)(nv, nv, nnz, 1, CHOLMOD_REAL, c); int_t *rowind = (int_t *)temp_m->i; int_t *colind = (int_t *)temp_m->j; @@ -120,126 +81,56 @@ cholmod_sparse *gen_laplacian(const net_t *instance, cholmod_common *c, temp_m->nnz = nnz; - for (uint_t i = 0; i < num_gverts; i++) { + for (uint_t i = 0; i < nv; i++) { rowind[i] = i; colind[i] = i; acoo[i] = 0; } - cholmod_sparse *adjacency = gen_adjacency(instance, false, true, (int)num_gverts - (int)num_verts, c); - int_t *ia = (int_t *)adjacency->p; - double *aa = (double *)adjacency->x; + cholmod_sparse *adjacency = gen_adjacency(net, false, true, true, c); - for (uint_t i = 0; i < num_verts; i++) { - for (uint_t j = 0; j < ia[i + 1] - ia[i]; j++) { - acoo[i] += aa[ia[i] + j]; - } - } + for (uint_t i = 0; i < ne; i++) { + if (!net->fuses[i]){ + uint_t v1 = ev[2 * i]; + uint_t v2 = ev[2 * i + 1]; - if (network->boundary == TORUS_BOUND) { - for (uint_t i = 0; i < network->bound_inds[1]; i++) { - uint_t vv = network->bound_verts[i]; - rowind[num_gverts + 2*i] = vv; - colind[num_gverts + 2*i] = network->nv + i; - acoo[num_gverts + 2*i] = -1; - rowind[num_gverts + 2*i+1] = network->nv + i; - colind[num_gverts + 2*i+1] = vv; - acoo[num_gverts + 2*i+1] = -1; - acoo[vv]++; - acoo[network->nv + i]++; + acoo[v1]++; + acoo[v2]++; } - acoo[bound_verts[0]]++; } - if (network->boundary != TORUS_BOUND) { - if (voltage_bound) { - for (uint_t i = 0; i < num_bounds; i++) { - for (uint_t j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) { - acoo[bound_verts[bound_inds[i] + j]]++; - - rowind[nnz - 1 - 2 * (bound_inds[i] + j)] = bound_verts[bound_inds[i] + j]; - colind[nnz - 1 - 2 * (bound_inds[i] + j)] = num_gverts - 1; - acoo[nnz - 1 - 2 * (bound_inds[i] + j)] = -1; - rowind[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = num_gverts - 1; - colind[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = bound_verts[bound_inds[i] + j]; - acoo[nnz - 1 - 2 * (bound_inds[i] + j) - 1] = -1; - acoo[num_gverts - 1]++; - acoo[bound_verts[bound_inds[i] + j]]++; - } - } - } else { - for (uint_t i = 0; i < num_bounds; i++) { - if (network->boundary != EMBEDDED_BOUND) acoo[0]++; - - rowind[num_verts + i] = num_verts + i; - colind[num_verts + i] = num_verts + i; - if (i < 2) - acoo[num_verts + i] = (network->bound_inds[i + 1] - - network->bound_inds[i] + network->num_bounds - 2) * - inf; - else - acoo[num_verts + i] = - (network->bound_inds[i + 1] - network->bound_inds[i] + 2) * inf; - } + if (net->voltage_bound && g->boundary != TORUS_BOUND) { + for (uint_t i = 0; i < net->dim; i++) { + uint_t v = g->nbi[i]; + for (uint_t j = 0; j < g->vei[v+1] - g->vei[v]; j++) { + uint_t e = g->ve[g->vei[v] + j]; + uint_t v0 = g->ev[2 * e]; + uint_t v1 = g->ev[2 * e + 1]; - uint_t start = num_gverts; - - for (uint_t i = 0; i < 2; i++) { - for (uint_t j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) { - rowind[start + bound_inds[i] + j] = bound_verts[bound_inds[i] + j]; - colind[start + bound_inds[i] + j] = num_verts + i; - acoo[start + bound_inds[i] + j] = -inf; - acoo[bound_verts[bound_inds[i] + j]] += inf; - colind[start + bound_inds[num_bounds] + bound_inds[i] + j] = - bound_verts[bound_inds[i] + j]; - rowind[start + bound_inds[num_bounds] + bound_inds[i] + j] = - num_verts + i; - acoo[start + bound_inds[num_bounds] + bound_inds[i] + j] = -inf; + if (g->bq[v0] || g->bq[v1]) { + acoo[i]++; } } - - for (uint_t i = 0; i < num_bounds - 2; i++) { - rowind[nnz - 2 * i - 1] = num_verts; - colind[nnz - 2 * i - 1] = num_verts + 2 + i; - acoo[nnz - 2 * i - 1] = -inf; - rowind[nnz - 2 * i - 2] = num_verts + 1; - colind[nnz - 2 * i - 2] = num_verts + 2 + i; - acoo[nnz - 2 * i - 2] = -inf; - colind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts; - rowind[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = num_verts + 2 + i; - acoo[nnz - 2 * i - 1 - 2 * (num_bounds - 2)] = -inf; - colind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 1; - rowind[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = num_verts + 2 + i; - acoo[nnz - 2 * i - 2 - 2 * (num_bounds - 2)] = -inf; - } } + } else { + acoo[0]++; } - - for (uint_t i = 0; i < num_gverts; i++) { - if (instance->marks[i] != instance->marks[network->nv]) - acoo[i]++; + for (uint_t i = 0; i < nv; i++) { + if (acoo[i] == 0) acoo[i]++; } - assert(CHOL_F(check_triplet)(temp_m, c)); + //assert(CHOL_F(check_triplet)(temp_m, c)); cholmod_sparse *t_out = CHOL_F(triplet_to_sparse)(temp_m, nnz, c); - assert(CHOL_F(check_sparse)(t_out, c)); + //assert(CHOL_F(check_sparse)(t_out, c)); double alpha[2] = {1, 0}; double beta[2] = {-1, 0}; - cholmod_sparse *t_laplacian = - CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c); - - cholmod_sparse *laplacian; - if (symmetric) - laplacian = CHOL_F(copy)(t_laplacian, 1, 1, c); - else - laplacian = CHOL_F(copy)(t_laplacian, 0, 1, c); + cholmod_sparse *laplacian = CHOL_F(add)(t_out, adjacency, alpha, beta, 1, 1, c); CHOL_F(free_sparse)(&t_out, c); CHOL_F(free_sparse)(&adjacency, c); - CHOL_F(free_sparse)(&t_laplacian, c); CHOL_F(free_triplet)(&temp_m, c); return laplacian; -- cgit v1.2.3-70-g09d2