From 2bb0740b68fdb62d45adc00204b3990ca42ade77 Mon Sep 17 00:00:00 2001 From: Jaron Kent-Dobias Date: Mon, 22 Aug 2016 10:11:14 -0400 Subject: started repo again without all the data files gunking the works --- src/gen_laplacian.c | 231 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 231 insertions(+) create mode 100644 src/gen_laplacian.c (limited to 'src/gen_laplacian.c') diff --git a/src/gen_laplacian.c b/src/gen_laplacian.c new file mode 100644 index 0000000..6f4263a --- /dev/null +++ b/src/gen_laplacian.c @@ -0,0 +1,231 @@ + +#include "fracture.h" + +cholmod_sparse *gen_adjacency(const finst *instance, bool dual, bool breakv, + unsigned int pad, cholmod_common *c) { + unsigned int ne; + if (dual) + ne = ((int)instance->network->num_edges - + (int)instance->num_remaining_edges); + else { + ne = instance->num_remaining_edges; + if (!breakv && instance->network->boundary != TORUS_BOUND) { + ne += instance->network->bound_inds[2]; + } + } + + unsigned int nnz = 2 * ne; + + unsigned int nv; + if (dual) + nv = instance->network->num_dual_verts; + else { + if (breakv) nv = instance->network->num_verts_break; + else nv = instance->network->num_verts; + if (instance->network->boundary != TORUS_BOUND && !breakv) { + nv += 2; + } + } + + cholmod_triplet *t = + CHOL_F(allocate_triplet)(nv + pad, nv + pad, nnz, 0, CHOLMOD_REAL, c); + CHOL_INT *ri = (CHOL_INT *)t->i; + CHOL_INT *ci = (CHOL_INT *)t->j; + double *ai = (double *)t->x; + + t->nnz = nnz; + + unsigned int *etv; + if (dual) + etv = instance->network->dual_edges_to_verts; + else { + if (breakv) etv = instance->network->edges_to_verts_break; + else etv = instance->network->edges_to_verts; + } + + unsigned int count = 0; + + for (unsigned int i = 0; i < instance->network->num_edges; i++) { + if ((instance->fuses[i] && dual) || (!instance->fuses[i] && !dual)) { + unsigned int v1 = etv[2 * i]; + unsigned int v2 = etv[2 * i + 1]; + + ri[2 * count] = v1; + ci[2 * count] = v2; + ai[2 * count] = 1; + + ri[2 * count + 1] = v2; + ci[2 * count + 1] = v1; + ai[2 * count + 1] = 1; + + count++; + } + } + + if (!breakv && instance->network->boundary != TORUS_BOUND && !dual) { + for (unsigned int i = 0; i < 2; i++) { + for (unsigned int j = 0; j < instance->network->bound_inds[i+1] - instance->network->bound_inds[i]; j++) { + ri[2*count] = instance->network->num_verts + i; + ci[2*count] = instance->network->bound_verts[instance->network->bound_inds[i] + j]; + ai[2*count] = 1; + ri[2*count+1] = instance->network->bound_verts[instance->network->bound_inds[i] + j]; + ci[2*count+1] = instance->network->num_verts + i; + ai[2*count+1] = 1; + count++; + } + } + } + + + cholmod_sparse *s = CHOL_F(triplet_to_sparse)(t, nnz, c); + CHOL_F(free_triplet)(&t, c); + + return s; +} + +cholmod_sparse *gen_laplacian(const finst *instance, cholmod_common *c, + bool symmetric) { + const fnet *network = instance->network; + unsigned int num_verts = network->num_verts_break; + double *vert_coords = network->vert_coords; + unsigned int num_bounds = network->num_bounds; + double inf = instance->inf; + bool voltage_bound = instance->voltage_bound; + unsigned int *bound_inds = network->bound_inds; + unsigned int *bound_verts = network->bound_verts; + + unsigned int num_gedges; + if (voltage_bound) { + num_gedges = 0; + } else { + num_gedges = network->bound_inds[num_bounds] + (num_bounds - 2) * 2; + } + if (network->boundary == TORUS_BOUND) num_gedges = bound_inds[1]; + unsigned int num_gverts = network->break_dim; + + unsigned int nnz = num_gverts + 2 * num_gedges; + + cholmod_triplet *temp_m = + CHOL_F(allocate_triplet)(num_gverts, num_gverts, nnz, 0, CHOLMOD_REAL, c); + + CHOL_INT *rowind = (CHOL_INT *)temp_m->i; + CHOL_INT *colind = (CHOL_INT *)temp_m->j; + double *acoo = (double *)temp_m->x; + + temp_m->nnz = nnz; + + for (unsigned int i = 0; i < num_gverts; 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); + CHOL_INT *ia = (CHOL_INT *)adjacency->p; + CHOL_INT *ja = (CHOL_INT *)adjacency->i; + double *aa = (double *)adjacency->x; + + for (unsigned int i = 0; i < num_verts; i++) { + for (unsigned int j = 0; j < ia[i + 1] - ia[i]; j++) { + acoo[i] += aa[ia[i] + j]; + } + } + + if (network->boundary == TORUS_BOUND) { + for (unsigned int i = 0; i < network->bound_inds[1]; i++) { + unsigned int vv = network->bound_verts[i]; + rowind[num_gverts + 2*i] = vv; + colind[num_gverts + 2*i] = network->num_verts + i; + acoo[num_gverts + 2*i] = -1; + rowind[num_gverts + 2*i+1] = network->num_verts + i; + colind[num_gverts + 2*i+1] = vv; + acoo[num_gverts + 2*i+1] = -1; + acoo[vv]++; + acoo[network->num_verts + i]++; + } + acoo[bound_verts[0]]++; + } + + if (network->boundary != TORUS_BOUND) { + if (voltage_bound) { + for (unsigned int i = 0; i < 2; i++) { + for (unsigned int j = 0; j < bound_inds[i + 1] - bound_inds[i]; j++) { + acoo[bound_verts[bound_inds[i] + j]]++; + } + } + } else { + acoo[0]++; + for (unsigned int i = 0; i < num_bounds; i++) { + 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; + } + + unsigned int start = num_gverts; + + for (unsigned int i = 0; i < num_bounds; i++) { + for (unsigned int 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; + } + } + + for (unsigned int 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; + } + } + } + + + for (unsigned int i = 0; i < num_gverts; i++) { + if (acoo[i] == 0) + acoo[i] = 1; + } + + 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)); + + 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); + + 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