diff options
Diffstat (limited to 'lib/src')
| -rw-r--r-- | lib/src/network.cpp | 417 | ||||
| -rw-r--r-- | lib/src/problem.cpp | 217 |
2 files changed, 301 insertions, 333 deletions
diff --git a/lib/src/network.cpp b/lib/src/network.cpp index 1d9b3f9..b081a3c 100644 --- a/lib/src/network.cpp +++ b/lib/src/network.cpp @@ -1,13 +1,6 @@ #include "network.hpp" - -class nanException: public std::exception -{ - virtual const char* what() const throw() - { - return "The linear problem returned NaN."; - } -} nanex; +#include <iostream> class nofuseException: public std::exception { @@ -17,210 +10,6 @@ class nofuseException: public std::exception } } nofuseex; -problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), voltcurmat(vcmat), c(c), axis(axis) { - b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); - for (unsigned i = 0; i < G.edges.size(); i++) { - graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; - graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; - - if (G.edges[i].crossings[axis]) { - bool ind; - if (axis == 1) { - ind = v0.y < v1.y ? 0 : 1; - } else { - ind = v0.x < v1.x ? 0 : 1; - } - - ((double *)b->x)[G.edges[i].v[ind]] += 1.0; - ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0; - } - } - - unsigned nnz = G.vertices.size() + G.edges.size(); - - cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c); - - for (unsigned i = 0; i < G.vertices.size(); i++) { - ((CHOL_INT *)t->i)[i] = i; - ((CHOL_INT *)t->j)[i] = i; - ((double *)t->x)[i] = 0.0; - } - - unsigned terms = G.vertices.size(); - - std::unordered_map<std::array<unsigned, 2>, unsigned> known_edges; - - for (unsigned i = 0; i < G.edges.size(); i++) { - unsigned v0 = G.edges[i].v[0]; - unsigned v1 = G.edges[i].v[1]; - - ((double *)t->x)[v0]++; - ((double *)t->x)[v1]++; - - unsigned s0 = v0 < v1 ? v0 : v1; - unsigned s1 = v0 < v1 ? v1 : v0; - - auto it = known_edges.find({s0, s1}); - - if (it == known_edges.end()) { - ((CHOL_INT *)t->i)[terms] = s1; - ((CHOL_INT *)t->j)[terms] = s0; - ((double *)t->x)[terms] = -1.0; - - known_edges[{s0, s1}] = terms; - terms++; - } else { - ((double *)t->x)[it->second] -= 1.0; - } - } - - ((double *)t->x)[0]++; - - t->nnz = terms; - - cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c); - CHOL_F(free_triplet)(&t, c); - factor = CHOL_F(analyze)(laplacian, c); - CHOL_F(factorize)(laplacian, factor, c); - CHOL_F(free_sparse)(&laplacian, c); -} - -problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) { - cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c); - - t->nnz = 2 * G.edges.size(); - - for (unsigned i = 0; i < G.edges.size(); i++) { - ((CHOL_INT *)t->i)[2 * i] = i; - ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0]; - ((double *)t->x)[2 * i] = 1.0; - - ((CHOL_INT *)t->i)[2 * i + 1] = i; - ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1]; - ((double *)t->x)[2 * i + 1] = -1.0; - } - - voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c); - - CHOL_F(free_triplet)(&t, c); -} - -problem::problem(const problem& other) : G(other.G), c(other.c), axis(other.axis) { - b = CHOL_F(copy_dense)(other.b, c); - factor = CHOL_F(copy_factor)(other.factor, c); - voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c); -} - -problem::~problem() { - CHOL_F(free_dense)(&b, c); - CHOL_F(free_factor)(&factor, c); - CHOL_F(free_sparse)(&voltcurmat, c); -} - -current_info problem::solve(std::vector<bool>& fuses) { - cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); - - if (((double *)x->x)[0] != ((double *)x->x)[0]) { - throw nanex; - } - - cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c); - - double alpha[2] = {1, 0}; - double beta[2] = {0, 0}; - CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); - - std::vector<double> currents(G.edges.size()); - - double total_current = 0; - - for (int i = 0; i < G.edges.size(); i++) { - if (fuses[i]) { - currents[i] = 0; - } else { - currents[i] = ((double *)y->x)[i]; - - graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; - graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; - - if (G.edges[i].crossings[axis]) { - bool comp; - if (axis == 1) { - comp = v0.y > v1.y; - } else { - comp = v0.x > v1.x; - } - - if (comp) { - currents[i] += 1.0; - total_current += currents[i]; - } else { - currents[i] -= 1.0; - total_current -= currents[i]; - } - } - } - } - - CHOL_F(free_dense)(&x, c); - CHOL_F(free_dense)(&y, c); - - return {total_current, currents}; -} - -void problem::break_edge(unsigned e, bool unbreak) { - unsigned v0 = G.edges[e].v[0]; - unsigned v1 = G.edges[e].v[1]; - - unsigned n = factor->n; - - cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); - - unsigned s1, s2; - s1 = v0 < v1 ? v0 : v1; - s2 = v0 < v1 ? v1 : v0; - - CHOL_INT *pp = (CHOL_INT *)update_mat->p; - CHOL_INT *ii = (CHOL_INT *)update_mat->i; - double *xx = (double *)update_mat->x; - - for (unsigned i = 0; i <= s1; i++) { - pp[i] = 0; - } - - for (unsigned 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, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c); - - CHOL_F(updown)(unbreak, perm_update_mat, factor, c); - - CHOL_F(free_sparse)(&perm_update_mat, c); - CHOL_F(free_sparse)(&update_mat, c); - - graph::coordinate r0 = G.vertices[v0].r; - graph::coordinate r1 = G.vertices[v1].r; - - if (G.edges[e].crossings[axis]) { - bool ind; - if (axis == 1) { - ind = r0.y < r1.y ? unbreak : !unbreak; - } else { - ind = r0.x < r1.x ? unbreak : !unbreak; - } - - ((double *)b->x)[G.edges[e].v[ind]] -= 1.0; - ((double *)b->x)[G.edges[e].v[!ind]] += 1.0; - } -} - network::network(const graph& G) : G(G), fuses(G.edges.size()), thresholds(G.edges.size()) {} network::network(const network& n) : G(n.G), fuses(n.fuses), thresholds(n.thresholds) {} @@ -244,16 +33,19 @@ void network::set_thresholds(double beta, std::mt19937& rng) { } } -fuse_network::fuse_network(const graph& G, cholmod_common *c) : network(G), ohm(G, 1, c) { -} - -void fuse_network::fracture(hooks& m, double cutoff) { +void network::fracture(hooks& m, bool one_axis) { m.pre_fracture(*this); while (true) { - current_info ci = ohm.solve(fuses); + current_info c = this->get_current_info(); + + double min_cond = 1.0 / G.edges.size(); - if (ci.conductivity < cutoff * G.vertices.size()) { + if (c.conductivity[0] < min_cond && c.conductivity[1] < min_cond) { + break; + } + + if (one_axis && (c.conductivity[0] < min_cond || c.conductivity[1] < min_cond)) { break; } @@ -261,8 +53,9 @@ void fuse_network::fracture(hooks& m, double cutoff) { long double max_val = std::numeric_limits<long double>::lowest(); for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && fabs(ci.currents[i]) > cutoff) { - long double val = logl(fabs(ci.currents[i])) - thresholds[i]; + if (!fuses[i] && c.currents[i] > min_cond) { + long double val = logl(c.currents[i]) - thresholds[i]; + if (val > max_val) { max_val = val; max_pos = i; @@ -274,172 +67,130 @@ void fuse_network::fracture(hooks& m, double cutoff) { throw nofuseex; } - fuses[max_pos] = true; - ohm.break_edge(max_pos); - - m.bond_broken(*this, ci, max_pos); + this->break_edge(max_pos); + m.bond_broken(*this, c, max_pos); } m.post_fracture(*this); } -elastic_network::elastic_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) { -} - -elastic_network::elastic_network(const elastic_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) { -} - -void elastic_network::fracture(hooks& m, double weight, double cutoff) { - this->weight = weight; - m.pre_fracture(*this); - - while (true) { - current_info ctot = this->get_current_info(); - - if (ctot.conductivity == 0) { - break; - } - - unsigned max_pos = UINT_MAX; - long double max_val = std::numeric_limits<long double>::lowest(); - - for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && ctot.currents[i] > 1.0 / G.edges.size()) { - long double val = logl(ctot.currents[i]) - thresholds[i]; - if (val > max_val) { - max_val = val; - max_pos = i; - } - } - } - - if (max_pos == UINT_MAX) { - throw nofuseex; - } - this->break_edge(max_pos); +fuse_network::fuse_network(const graph& G, cholmod_common* c, double weight) : + network(G), ohm_x(G, 0, c), ohm_y(G, 1, c) {} - m.bond_broken(*this, ctot, max_pos); - } +fuse_network::fuse_network(const fuse_network& n) : + network(n), ohm_x(n.ohm_x), ohm_y(n.ohm_y) { +} - m.post_fracture(*this); +void fuse_network::break_edge(unsigned e, bool unbreak) { + fuses[e] = !unbreak; + ohm_x.break_edge(e, unbreak); + ohm_y.break_edge(e, unbreak); } -current_info elastic_network::get_current_info() { - current_info cx = hook_x.solve(fuses); - current_info cy = hook_y.solve(fuses); +current_info fuse_network::get_current_info() { + current_info cx = ohm_x.solve(fuses); + current_info cy = ohm_y.solve(fuses); - bool done_x = cx.conductivity < 1.0 / G.edges.size(); - bool done_y = cy.conductivity < 1.0 / G.edges.size(); + bool done_x = cx.conductivity[0] < 1.0 / G.edges.size(); + bool done_y = cy.conductivity[1] < 1.0 / G.edges.size(); current_info ctot; ctot.currents.resize(G.edges.size()); + ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]}; - if (done_x && done_y) { - ctot.conductivity = 0; - } else if (done_x) { - ctot.conductivity = cy.conductivity; + if (done_x && !done_y) { for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity; + ctot.currents[i] = fabs(weight * cy.currents[i] / cy.conductivity[1]); } - } else if (done_y) { - ctot.conductivity = cx.conductivity; + } else if (done_y && !done_x) { for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity; + ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0]); } - } else { - ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); + } else if (!done_x && !done_y) { for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2)); + ctot.currents[i] = fabs((1 - weight) * cx.currents[i] / cx.conductivity[0] + + weight * cy.currents[i] / cy.conductivity[1]); } } return ctot; } + +elastic_network::elastic_network(const graph& G, cholmod_common* c, double weight) : + network(G), weight(weight), + hook_x(G, 0, c), hook_y(G, 1, c) {} + +elastic_network::elastic_network(const elastic_network& n) : + network(n), weight(n.weight), hook_x(n.hook_x), hook_y(n.hook_y) {} + void elastic_network::break_edge(unsigned e, bool unbreak) { fuses[e] = !unbreak; hook_x.break_edge(e, unbreak); hook_y.break_edge(e, unbreak); } -percolation_network::percolation_network(const graph& G, cholmod_common *c) : network(G), hook_x(G, 1, c), hook_y(G, 0, c) { -} +current_info elastic_network::get_current_info() { + current_info cx = hook_x.solve(fuses); + current_info cy = hook_y.solve(fuses); -percolation_network::percolation_network(const percolation_network& n) : network(n), hook_x(n.hook_x), hook_y(n.hook_y) { -} + bool done_x = cx.conductivity[0] < 1.0 / G.edges.size(); + bool done_y = cy.conductivity[1] < 1.0 / G.edges.size(); -void percolation_network::fracture(hooks& m, double weight, double cutoff) { - this->weight = weight; - m.pre_fracture(*this); - - while (true) { - current_info ctot = this->get_current_info(); + current_info ctot; + ctot.currents.resize(G.edges.size()); + ctot.conductivity[0] = cx.conductivity[0]; + ctot.conductivity[1] = cy.conductivity[1]; - if (ctot.conductivity == 0) { - break; + if (done_x && !done_y) { + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity[1]; } - - unsigned max_pos = UINT_MAX; - long double max_val = std::numeric_limits<long double>::lowest(); - + } else if (done_y && !done_x) { for (unsigned i = 0; i < G.edges.size(); i++) { - if (!fuses[i] && ctot.currents[i] > 1.0 / G.edges.size()) { - long double val = - thresholds[i]; - if (val > max_val) { - max_val = val; - max_pos = i; - } - } + ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity[0]; } - - if (max_pos == UINT_MAX) { - throw nofuseex; + } else if (!done_x && !done_y) { + for (unsigned i = 0; i < G.edges.size(); i++) { + ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity[0], 2) + + pow(weight * cy.currents[i] / cy.conductivity[1], 2)); } - - this->break_edge(max_pos); - - m.bond_broken(*this, ctot, max_pos); } - m.post_fracture(*this); + return ctot; } -current_info percolation_network::get_current_info() { - current_info cx = hook_x.solve(fuses); - current_info cy = hook_y.solve(fuses); - bool done_x = cx.conductivity < 1.0 / G.edges.size(); - bool done_y = cy.conductivity < 1.0 / G.edges.size(); +percolation_network::percolation_network(const graph& G, cholmod_common* c) : + network(G), px(G, 0, c), py(G, 1, c) {} + +percolation_network::percolation_network(const percolation_network& n) : + network(n), px(n.px), py(n.py) {} +current_info percolation_network::get_current_info() { current_info ctot; - ctot.currents.resize(G.edges.size()); + ctot.currents.resize(G.edges.size(), 0); - if (done_x && done_y) { - ctot.conductivity = 0; - } else if (done_x) { - ctot.conductivity = cy.conductivity; - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = weight * fabs(cy.currents[i]) / cy.conductivity; - } - } else if (done_y) { - ctot.conductivity = cx.conductivity; - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = (1 - weight) * fabs(cx.currents[i]) / cx.conductivity; - } - } else { - ctot.conductivity = sqrt(pow((1 - weight) * cx.conductivity, 2) + pow(weight * cy.conductivity, 2)); - for (unsigned i = 0; i < G.edges.size(); i++) { - ctot.currents[i] = sqrt(pow((1 - weight) * cx.currents[i] / cx.conductivity, 2) + pow(weight * cy.currents[i] / cy.conductivity, 2)); + current_info cx = px.solve(fuses); + current_info cy = py.solve(fuses); + + ctot.conductivity = {cx.conductivity[0], cy.conductivity[1]}; + + double min_cur = 1.0 / G.edges.size(); + + for (unsigned i = 0; i < G.edges.size(); i++) { + if (fabs(cx.currents[i]) > min_cur || fabs(cy.currents[i]) > min_cur) { + ctot.currents[i] = 1.0; } } - + return ctot; } void percolation_network::break_edge(unsigned e, bool unbreak) { fuses[e] = !unbreak; - hook_x.break_edge(e, unbreak); - hook_y.break_edge(e, unbreak); + px.break_edge(e, unbreak); + py.break_edge(e, unbreak); } diff --git a/lib/src/problem.cpp b/lib/src/problem.cpp new file mode 100644 index 0000000..f66a35c --- /dev/null +++ b/lib/src/problem.cpp @@ -0,0 +1,217 @@ + +#include "problem.hpp" + +class nanException: public std::exception +{ + virtual const char* what() const throw() + { + return "The linear problem returned NaN."; + } +} nanex; + +problem::problem(const graph& G, unsigned axis, cholmod_sparse *vcmat, cholmod_common *c) : G(G), axis(axis), voltcurmat(vcmat), c(c) { + b = CHOL_F(zeros)(G.vertices.size(), 1, CHOLMOD_REAL, c); + for (unsigned i = 0; i < G.edges.size(); i++) { + graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; + graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; + + if (G.edges[i].crossings[axis]) { + bool ind; + if (axis == 1) { + ind = v0.y < v1.y ? 0 : 1; + } else { + ind = v0.x < v1.x ? 0 : 1; + } + + ((double *)b->x)[G.edges[i].v[ind]] += 1.0; + ((double *)b->x)[G.edges[i].v[!ind]] -= 1.0; + } + } + + unsigned nnz = G.vertices.size() + G.edges.size(); + + cholmod_triplet *t = CHOL_F(allocate_triplet)(G.vertices.size(), G.vertices.size(), nnz, 1, CHOLMOD_REAL, c); + + for (unsigned i = 0; i < G.vertices.size(); i++) { + ((CHOL_INT *)t->i)[i] = i; + ((CHOL_INT *)t->j)[i] = i; + ((double *)t->x)[i] = 0.0; + } + + unsigned terms = G.vertices.size(); + + std::unordered_map<std::array<unsigned, 2>, unsigned> known_edges; + + for (unsigned i = 0; i < G.edges.size(); i++) { + unsigned v0 = G.edges[i].v[0]; + unsigned v1 = G.edges[i].v[1]; + + ((double *)t->x)[v0]++; + ((double *)t->x)[v1]++; + + unsigned s0 = v0 < v1 ? v0 : v1; + unsigned s1 = v0 < v1 ? v1 : v0; + + auto it = known_edges.find({s0, s1}); + + if (it == known_edges.end()) { + ((CHOL_INT *)t->i)[terms] = s1; + ((CHOL_INT *)t->j)[terms] = s0; + ((double *)t->x)[terms] = -1.0; + + known_edges[{s0, s1}] = terms; + terms++; + } else { + ((double *)t->x)[it->second] -= 1.0; + } + } + + ((double *)t->x)[0]++; + + t->nnz = terms; + + cholmod_sparse *laplacian = CHOL_F(triplet_to_sparse)(t, terms, c); + CHOL_F(free_triplet)(&t, c); + factor = CHOL_F(analyze)(laplacian, c); + CHOL_F(factorize)(laplacian, factor, c); + CHOL_F(free_sparse)(&laplacian, c); +} + +problem::problem(const graph& G, unsigned axis, cholmod_common *c) : problem(G, axis, NULL, c) { + cholmod_triplet *t = CHOL_F(allocate_triplet)(G.edges.size(), G.vertices.size(), 2 * G.edges.size(), 0, CHOLMOD_REAL, c); + + t->nnz = 2 * G.edges.size(); + + for (unsigned i = 0; i < G.edges.size(); i++) { + ((CHOL_INT *)t->i)[2 * i] = i; + ((CHOL_INT *)t->j)[2 * i] = G.edges[i].v[0]; + ((double *)t->x)[2 * i] = 1.0; + + ((CHOL_INT *)t->i)[2 * i + 1] = i; + ((CHOL_INT *)t->j)[2 * i + 1] = G.edges[i].v[1]; + ((double *)t->x)[2 * i + 1] = -1.0; + } + + voltcurmat = CHOL_F(triplet_to_sparse)(t, 2 * G.edges.size(), c); + + CHOL_F(free_triplet)(&t, c); +} + +problem::problem(const problem& other) : G(other.G), axis(other.axis), c(other.c) { + b = CHOL_F(copy_dense)(other.b, c); + factor = CHOL_F(copy_factor)(other.factor, c); + voltcurmat = CHOL_F(copy_sparse)(other.voltcurmat, c); +} + +problem::~problem() { + CHOL_F(free_dense)(&b, c); + CHOL_F(free_factor)(&factor, c); + CHOL_F(free_sparse)(&voltcurmat, c); +} + +current_info problem::solve(std::vector<bool>& fuses) { + cholmod_dense *x = CHOL_F(solve)(CHOLMOD_A, factor, b, c); + + if (((double *)x->x)[0] != ((double *)x->x)[0]) { + throw nanex; + } + + cholmod_dense *y = CHOL_F(allocate_dense)(G.edges.size(), 1, G.edges.size(), CHOLMOD_REAL, c); + + double alpha[2] = {1, 0}; + double beta[2] = {0, 0}; + CHOL_F(sdmult)(voltcurmat, 0, alpha, beta, x, y, c); + + std::vector<double> currents(G.edges.size()); + + std::array<double, 2> total_current = {0, 0}; + + for (int i = 0; i < G.edges.size(); i++) { + if (fuses[i]) { + currents[i] = 0; + } else { + currents[i] = ((double *)y->x)[i]; + + graph::coordinate v0 = G.vertices[G.edges[i].v[0]].r; + graph::coordinate v1 = G.vertices[G.edges[i].v[1]].r; + + if (G.edges[i].crossings[axis]) { + bool comp; + if (axis == 1) { + comp = v0.y > v1.y; + } else { + comp = v0.x > v1.x; + } + + if (comp) { + currents[i] += 1.0; + total_current[axis] += currents[i]; + } else { + currents[i] -= 1.0; + total_current[axis] -= currents[i]; + } + } + } + } + + total_current[!axis] = 0.0; + + CHOL_F(free_dense)(&x, c); + CHOL_F(free_dense)(&y, c); + + return {total_current, currents}; +} + +void problem::break_edge(unsigned e, bool unbreak) { + unsigned v0 = G.edges[e].v[0]; + unsigned v1 = G.edges[e].v[1]; + + unsigned n = factor->n; + + cholmod_sparse *update_mat = CHOL_F(allocate_sparse)(n, n, 2, true, true, 0, CHOLMOD_REAL, c); + + unsigned s1, s2; + s1 = v0 < v1 ? v0 : v1; + s2 = v0 < v1 ? v1 : v0; + + CHOL_INT *pp = (CHOL_INT *)update_mat->p; + CHOL_INT *ii = (CHOL_INT *)update_mat->i; + double *xx = (double *)update_mat->x; + + for (unsigned i = 0; i <= s1; i++) { + pp[i] = 0; + } + + for (unsigned 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, (CHOL_INT *)factor->Perm, factor->n, NULL, -1, true, true, c); + + CHOL_F(updown)(unbreak, perm_update_mat, factor, c); + + CHOL_F(free_sparse)(&perm_update_mat, c); + CHOL_F(free_sparse)(&update_mat, c); + + graph::coordinate r0 = G.vertices[v0].r; + graph::coordinate r1 = G.vertices[v1].r; + + if (G.edges[e].crossings[axis]) { + bool ind; + if (axis == 1) { + ind = r0.y < r1.y ? unbreak : !unbreak; + } else { + ind = r0.x < r1.x ? unbreak : !unbreak; + } + + ((double *)b->x)[G.edges[e].v[ind]] -= 1.0; + ((double *)b->x)[G.edges[e].v[!ind]] += 1.0; + } +} + |