started repo again without all the data files gunking the works

1 files changed, 231 insertions, 0 deletions

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;
+}