-#include "general.h"
-#include "SparseMatrix.h"
-#include "Multilevel.h"
+/* vim:set shiftwidth=4 ts=8: */
+
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+#include "Multilevel.h"
+#include "memory.h"
+#include "logic.h"
+#include "assert.h"
+#include "arith.h"
+
+Multilevel_control *Multilevel_control_new()
+{
+ Multilevel_control *ctrl;
+
+ ctrl = GNEW(Multilevel_control);
+ ctrl->minsize = 4;
+ ctrl->min_coarsen_factor = 0.75;
+ ctrl->maxlevel = 1 << 30;
+ ctrl->randomize = TRUE;
+ ctrl->coarsen_scheme =
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST;
+ ctrl->coarsen_scheme = COARSEN_HYBRID;
+ return ctrl;
+}
-Multilevel_control Multilevel_control_new(){
- Multilevel_control ctrl;
+void Multilevel_control_delete(Multilevel_control * ctrl)
+{
+ free(ctrl);
+}
- ctrl = MALLOC(sizeof(struct Multilevel_control_struct));
- ctrl->minsize = 4;
- ctrl->min_coarsen_factor = 0.75;
- ctrl->maxlevel = 1<<30;
- ctrl->randomize = TRUE;
- ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST;
- ctrl->coarsen_scheme = COARSEN_HYBRID;
- return ctrl;
+static Multilevel *Multilevel_init(SparseMatrix * A, real * node_weights)
+{
+ Multilevel *grid;
+ if (!A)
+ return NULL;
+ assert(A->m == A->n);
+ grid = GNEW(Multilevel);
+ grid->level = 0;
+ grid->n = A->n;
+ grid->A = A;
+ grid->P = NULL;
+ grid->R = NULL;
+ grid->node_weights = node_weights;
+ grid->next = NULL;
+ grid->prev = NULL;
+ grid->delete_top_level_A = FALSE;
+ return grid;
}
-void Multilevel_control_delete(Multilevel_control ctrl){
- FREE(ctrl);
+void Multilevel_delete(Multilevel * grid)
+{
+ if (!grid)
+ return;
+ if (grid->A) {
+ if (grid->level == 0) {
+ if (grid->delete_top_level_A)
+ SparseMatrix_delete(grid->A);
+ } else {
+ SparseMatrix_delete(grid->A);
+ }
+ }
+ SparseMatrix_delete(grid->P);
+ SparseMatrix_delete(grid->R);
+ if (grid->node_weights && grid->level > 0)
+ free(grid->node_weights);
+ Multilevel_delete(grid->next);
+ free(grid);
}
-static Multilevel Multilevel_init(SparseMatrix A, real *node_weights){
- Multilevel grid;
- if (!A) return NULL;
- assert(A->m == A->n);
- grid = MALLOC(sizeof(struct Multilevel_struct));
- grid->level = 0;
- grid->n = A->n;
- grid->A = A;
- grid->P = NULL;
- grid->R = NULL;
- grid->node_weights = node_weights;
- grid->next = NULL;
- grid->prev = NULL;
- grid->delete_top_level_A = FALSE;
- return grid;
+static int irand(int n){
+ /* 0, 1, ..., n-1 */
+ assert(n > 1);
+ /*return (int) MIN(floor(drand()*n),n-1);*/
+ return rand()%n;
}
-void Multilevel_delete(Multilevel grid){
- if (!grid) return;
- if (grid->A){
- if (grid->level == 0) {
- if (grid->delete_top_level_A) SparseMatrix_delete(grid->A);
- } else {
- SparseMatrix_delete(grid->A);
+static int *random_permutation(int n)
+{
+ int *p;
+ int i, j, pp, len;
+ if (n <= 0)
+ return NULL;
+ p = N_GNEW(n, int);
+ for (i = 0; i < n; i++)
+ p[i] = i;
+
+ len = n;
+ while (len > 1) {
+ j = irand(len);
+ pp = p[len - 1];
+ p[len - 1] = p[j];
+ p[j] = pp;
+ len--;
}
- }
- SparseMatrix_delete(grid->P);
- SparseMatrix_delete(grid->R);
- if (grid->node_weights && grid->level > 0) FREE(grid->node_weights);
- Multilevel_delete(grid->next);
- FREE(grid);
+ return p;
}
-static void maximal_independent_vertex_set(SparseMatrix A, int randomize, int **vset, int *nvset, int *nzc){
- int i, ii, j, *ia, *ja, m, n, *p = NULL;
- assert(A);
- assert(SparseMatrix_known_strucural_symmetric(A));
- ia = A->ia;
- ja = A->ja;
- m = A->m;
- n = A->n;
- assert(n == m);
- *vset = MALLOC(sizeof(int)*m);
- for (i = 0; i < m; i++) (*vset)[i] = MAX_IND_VTX_SET_0;
- *nvset = 0;
- *nzc = 0;
-
- if (!randomize){
- for (i = 0; i < m; i++){
- if ((*vset)[i] == MAX_IND_VTX_SET_0){
- (*vset)[i] = (*nvset)++;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- (*vset)[ja[j]] = MAX_IND_VTX_SET_U;
- (*nzc)++;
- }
- }
- }
- } else {
- p = random_permutation(m);
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- if ((*vset)[i] == MAX_IND_VTX_SET_0){
- (*vset)[i] = (*nvset)++;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- (*vset)[ja[j]] = MAX_IND_VTX_SET_U;
- (*nzc)++;
- }
- }
+static void maximal_independent_vertex_set(SparseMatrix * A, int randomize,
+ int **vset, int *nvset,
+ int *nzc)
+{
+ int i, ii, j, *ia, *ja, m, n, *p = NULL;
+ assert(A);
+ assert(SparseMatrix_known_strucural_symmetric(A));
+ ia = A->ia;
+ ja = A->ja;
+ m = A->m;
+ n = A->n;
+ assert(n == m);
+ *vset = N_GNEW(m, int);
+ for (i = 0; i < m; i++)
+ (*vset)[i] = MAX_IND_VTX_SET_0;
+ *nvset = 0;
+ *nzc = 0;
+
+ if (!randomize) {
+ for (i = 0; i < m; i++) {
+ if ((*vset)[i] == MAX_IND_VTX_SET_0) {
+ (*vset)[i] = (*nvset)++;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ (*vset)[ja[j]] = MAX_IND_VTX_SET_U;
+ (*nzc)++;
+ }
+ }
+ }
+ } else {
+ p = random_permutation(m);
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ if ((*vset)[i] == MAX_IND_VTX_SET_0) {
+ (*vset)[i] = (*nvset)++;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ (*vset)[ja[j]] = MAX_IND_VTX_SET_U;
+ (*nzc)++;
+ }
+ }
+ }
+ free(p);
}
- FREE(p);
- }
- (*nzc) += *nvset;
+ (*nzc) += *nvset;
}
-static void maximal_independent_edge_set(SparseMatrix A, int randomize, int **matching, int *nmatch){
- int i, ii, j, *ia, *ja, m, n, *p = NULL;
- assert(A);
- assert(SparseMatrix_known_strucural_symmetric(A));
- ia = A->ia;
- ja = A->ja;
- m = A->m;
- n = A->n;
- assert(n == m);
- *matching = MALLOC(sizeof(int)*m);
- for (i = 0; i < m; i++) (*matching)[i] = i;
- *nmatch = n;
-
- if (!randomize){
- for (i = 0; i < m; i++){
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
- (*matching)[ja[j]] = i;
- (*matching)[i] = ja[j];
- (*nmatch)--;
- }
- }
- }
- } else {
- p = random_permutation(m);
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
- (*matching)[ja[j]] = i;
- (*matching)[i] = ja[j];
- (*nmatch)--;
- }
- }
+static void maximal_independent_edge_set(SparseMatrix * A, int randomize,
+ int **matching, int *nmatch)
+{
+ int i, ii, j, *ia, *ja, m, n, *p = NULL;
+ assert(A);
+ assert(SparseMatrix_known_strucural_symmetric(A));
+ ia = A->ia;
+ ja = A->ja;
+ m = A->m;
+ n = A->n;
+ assert(n == m);
+ *matching = N_GNEW(m, int);
+ for (i = 0; i < m; i++)
+ (*matching)[i] = i;
+ *nmatch = n;
+
+ if (!randomize) {
+ for (i = 0; i < m; i++) {
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i) {
+ (*matching)[ja[j]] = i;
+ (*matching)[i] = ja[j];
+ (*nmatch)--;
+ }
+ }
+ }
+ } else {
+ p = random_permutation(m);
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i) {
+ (*matching)[ja[j]] = i;
+ (*matching)[i] = ja[j];
+ (*nmatch)--;
+ }
+ }
+ }
+ free(p);
}
- FREE(p);
- }
}
-static void maximal_independent_edge_set_heavest_edge_pernode(SparseMatrix A, int randomize, int **matching, int *nmatch){
- int i, ii, j, *ia, *ja, m, n, *p = NULL;
- real *a, amax = 0;
- int first = TRUE, jamax = 0;
-
- assert(A);
- assert(SparseMatrix_known_strucural_symmetric(A));
- ia = A->ia;
- ja = A->ja;
- m = A->m;
- n = A->n;
- assert(n == m);
- *matching = MALLOC(sizeof(int)*m);
- for (i = 0; i < m; i++) (*matching)[i] = i;
- *nmatch = n;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
- assert(A->type == MATRIX_TYPE_REAL);
-
- a = (real*) A->a;
- if (!randomize){
- for (i = 0; i < m; i++){
- first = TRUE;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
- if (first) {
- amax = a[j];
- jamax = ja[j];
- first = FALSE;
- } else {
- if (a[j] > amax){
- amax = a[j];
- jamax = ja[j];
- }
- }
- }
- }
- if (!first){
- (*matching)[jamax] = i;
- (*matching)[i] = jamax;
- (*nmatch)--;
- }
- }
- } else {
- p = random_permutation(m);
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- if ((*matching)[i] != i) continue;
- first = TRUE;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
- if (first) {
- amax = a[j];
- jamax = ja[j];
- first = FALSE;
- } else {
- if (a[j] > amax){
- amax = a[j];
- jamax = ja[j];
- }
- }
- }
- }
- if (!first){
- (*matching)[jamax] = i;
- (*matching)[i] = jamax;
- (*nmatch)--;
- }
+static void maximal_independent_edge_set_heavest_edge_pernode(SparseMatrix
+ * A,
+ int
+ randomize,
+ int
+ **matching,
+ int *nmatch)
+{
+ int i, ii, j, *ia, *ja, m, n, *p = NULL;
+ real *a, amax = 0;
+ int first = TRUE, jamax = 0;
+
+ assert(A);
+ assert(SparseMatrix_known_strucural_symmetric(A));
+ ia = A->ia;
+ ja = A->ja;
+ m = A->m;
+ n = A->n;
+ assert(n == m);
+ *matching = N_GNEW(m, int);
+ for (i = 0; i < m; i++)
+ (*matching)[i] = i;
+ *nmatch = n;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+ assert(A->type == MATRIX_TYPE_REAL);
+
+ a = (real *) A->a;
+ if (!randomize) {
+ for (i = 0; i < m; i++) {
+ first = TRUE;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i) {
+ if (first) {
+ amax = a[j];
+ jamax = ja[j];
+ first = FALSE;
+ } else {
+ if (a[j] > amax) {
+ amax = a[j];
+ jamax = ja[j];
+ }
+ }
+ }
+ }
+ if (!first) {
+ (*matching)[jamax] = i;
+ (*matching)[i] = jamax;
+ (*nmatch)--;
+ }
+ }
+ } else {
+ p = random_permutation(m);
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ if ((*matching)[i] != i)
+ continue;
+ first = TRUE;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i) {
+ if (first) {
+ amax = a[j];
+ jamax = ja[j];
+ first = FALSE;
+ } else {
+ if (a[j] > amax) {
+ amax = a[j];
+ jamax = ja[j];
+ }
+ }
+ }
+ }
+ if (!first) {
+ (*matching)[jamax] = i;
+ (*matching)[i] = jamax;
+ (*nmatch)--;
+ }
+ }
+ free(p);
}
- FREE(p);
- }
}
#define node_degree(i) (ia[(i)+1] - ia[(i)])
-static void maximal_independent_edge_set_heavest_edge_pernode_leaves_first(SparseMatrix A, int randomize, int **cluster, int **clusterp, int *ncluster){
- int i, ii, j, *ia, *ja, m, n, *p = NULL, q;
- real *a, amax = 0;
- int first = TRUE, jamax = 0;
- int *matched, nz, ncmax = 0, nz0, nzz,k ;
- enum {UNMATCHED = -2, MATCHED = -1};
-
- assert(A);
- assert(SparseMatrix_known_strucural_symmetric(A));
- ia = A->ia;
- ja = A->ja;
- m = A->m;
- n = A->n;
- assert(n == m);
- *cluster = MALLOC(sizeof(int)*m);
- *clusterp = MALLOC(sizeof(int)*(m+1));
- matched = MALLOC(sizeof(int)*m);
-
- for (i = 0; i < m; i++) matched[i] = i;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
- assert(A->type == MATRIX_TYPE_REAL);
-
- *ncluster = 0;
- (*clusterp)[0] = 0;
- nz = 0;
- a = (real*) A->a;
- if (!randomize){
- for (i = 0; i < m; i++){
- if (matched[i] == MATCHED || node_degree(i) != 1) continue;
- q = ja[ia[i]];
- assert(matched[q] != MATCHED);
- matched[q] = MATCHED;
- (*cluster)[nz++] = q;
- for (j = ia[q]; j < ia[q+1]; j++){
- if (q == ja[j]) continue;
- if (node_degree(ja[j]) == 1){
- matched[ja[j]] = MATCHED;
- (*cluster)[nz++] = ja[j];
- }
- }
- ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]);
- nz0 = (*clusterp)[*ncluster];
- if (nz - nz0 <= MAX_CLUSTER_SIZE){
- (*clusterp)[++(*ncluster)] = nz;
- } else {
- (*clusterp)[++(*ncluster)] = ++nz0;
- nzz = nz0;
- for (k = nz0; k < nz && nzz < nz; k++){
- nzz += MAX_CLUSTER_SIZE - 1;
- nzz = MIN(nz, nzz);
- (*clusterp)[++(*ncluster)] = nzz;
- }
- }
+static void
+maximal_independent_edge_set_heavest_edge_pernode_leaves_first(SparseMatrix
+ * A,
+ int
+ randomize,
+ int
+ **cluster,
+ int
+ **clusterp,
+ int
+ *ncluster)
+{
+ int i, ii, j, *ia, *ja, m, n, *p = NULL, q;
+ real *a, amax = 0;
+ int first = TRUE, jamax = 0;
+ int *matched, nz, ncmax = 0, nz0, nzz, k;
+ enum { UNMATCHED = -2, MATCHED = -1 };
+
+ assert(A);
+ assert(SparseMatrix_known_strucural_symmetric(A));
+ ia = A->ia;
+ ja = A->ja;
+ m = A->m;
+ n = A->n;
+ assert(n == m);
+ *cluster = N_GNEW(m, int);
+ *clusterp = N_GNEW(m + 1, int);
+ matched = N_GNEW(m, int);
+
+ for (i = 0; i < m; i++)
+ matched[i] = i;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+ assert(A->type == MATRIX_TYPE_REAL);
+
+ *ncluster = 0;
+ (*clusterp)[0] = 0;
+ nz = 0;
+ a = (real *) A->a;
+ if (!randomize) {
+ for (i = 0; i < m; i++) {
+ if (matched[i] == MATCHED || node_degree(i) != 1)
+ continue;
+ q = ja[ia[i]];
+ assert(matched[q] != MATCHED);
+ matched[q] = MATCHED;
+ (*cluster)[nz++] = q;
+ for (j = ia[q]; j < ia[q + 1]; j++) {
+ if (q == ja[j])
+ continue;
+ if (node_degree(ja[j]) == 1) {
+ matched[ja[j]] = MATCHED;
+ (*cluster)[nz++] = ja[j];
+ }
+ }
+ ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]);
+ nz0 = (*clusterp)[*ncluster];
+ if (nz - nz0 <= MAX_CLUSTER_SIZE) {
+ (*clusterp)[++(*ncluster)] = nz;
+ } else {
+ (*clusterp)[++(*ncluster)] = ++nz0;
+ nzz = nz0;
+ for (k = nz0; k < nz && nzz < nz; k++) {
+ nzz += MAX_CLUSTER_SIZE - 1;
+ nzz = MIN(nz, nzz);
+ (*clusterp)[++(*ncluster)] = nzz;
+ }
+ }
- }
- #ifdef DEBUG_print
- if (Verbose)
- fprintf(stderr, "%d leaves and parents for %d clusters, largest cluster = %d\n",nz, *ncluster, ncmax);
+ }
+#ifdef DEBUG_print
+ if (Verbose)
+ fprintf(stderr,
+ "%d leaves and parents for %d clusters, largest cluster = %d\n",
+ nz, *ncluster, ncmax);
#endif
- for (i = 0; i < m; i++){
- first = TRUE;
- if (matched[i] == MATCHED) continue;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
- if (first) {
- amax = a[j];
- jamax = ja[j];
- first = FALSE;
- } else {
- if (a[j] > amax){
- amax = a[j];
- jamax = ja[j];
- }
- }
- }
- }
- if (!first){
- matched[jamax] = MATCHED;
- matched[i] = MATCHED;
- (*cluster)[nz++] = i;
- (*cluster)[nz++] = jamax;
- (*clusterp)[++(*ncluster)] = nz;
- }
- }
+ for (i = 0; i < m; i++) {
+ first = TRUE;
+ if (matched[i] == MATCHED)
+ continue;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if (matched[ja[j]] != MATCHED && matched[i] != MATCHED) {
+ if (first) {
+ amax = a[j];
+ jamax = ja[j];
+ first = FALSE;
+ } else {
+ if (a[j] > amax) {
+ amax = a[j];
+ jamax = ja[j];
+ }
+ }
+ }
+ }
+ if (!first) {
+ matched[jamax] = MATCHED;
+ matched[i] = MATCHED;
+ (*cluster)[nz++] = i;
+ (*cluster)[nz++] = jamax;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
+ }
- /* dan yi dian, wu ban */
- for (i = 0; i < m; i++){
- if (matched[i] == i){
- (*cluster)[nz++] = i;
- (*clusterp)[++(*ncluster)] = nz;
- }
- }
- assert(nz == n);
-
- } else {
- p = random_permutation(m);
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- if (matched[i] == MATCHED || node_degree(i) != 1) continue;
- q = ja[ia[i]];
- assert(matched[q] != MATCHED);
- matched[q] = MATCHED;
- (*cluster)[nz++] = q;
- for (j = ia[q]; j < ia[q+1]; j++){
- if (q == ja[j]) continue;
- if (node_degree(ja[j]) == 1){
- matched[ja[j]] = MATCHED;
- (*cluster)[nz++] = ja[j];
- }
- }
- ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]);
- nz0 = (*clusterp)[*ncluster];
- if (nz - nz0 <= MAX_CLUSTER_SIZE){
- (*clusterp)[++(*ncluster)] = nz;
- } else {
- (*clusterp)[++(*ncluster)] = ++nz0;
- nzz = nz0;
- for (k = nz0; k < nz && nzz < nz; k++){
- nzz += MAX_CLUSTER_SIZE - 1;
- nzz = MIN(nz, nzz);
- (*clusterp)[++(*ncluster)] = nzz;
- }
- }
- }
+ /* dan yi dian, wu ban */
+ for (i = 0; i < m; i++) {
+ if (matched[i] == i) {
+ (*cluster)[nz++] = i;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
+ }
+ assert(nz == n);
+
+ } else {
+ p = random_permutation(m);
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ if (matched[i] == MATCHED || node_degree(i) != 1)
+ continue;
+ q = ja[ia[i]];
+ assert(matched[q] != MATCHED);
+ matched[q] = MATCHED;
+ (*cluster)[nz++] = q;
+ for (j = ia[q]; j < ia[q + 1]; j++) {
+ if (q == ja[j])
+ continue;
+ if (node_degree(ja[j]) == 1) {
+ matched[ja[j]] = MATCHED;
+ (*cluster)[nz++] = ja[j];
+ }
+ }
+ ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]);
+ nz0 = (*clusterp)[*ncluster];
+ if (nz - nz0 <= MAX_CLUSTER_SIZE) {
+ (*clusterp)[++(*ncluster)] = nz;
+ } else {
+ (*clusterp)[++(*ncluster)] = ++nz0;
+ nzz = nz0;
+ for (k = nz0; k < nz && nzz < nz; k++) {
+ nzz += MAX_CLUSTER_SIZE - 1;
+ nzz = MIN(nz, nzz);
+ (*clusterp)[++(*ncluster)] = nzz;
+ }
+ }
+ }
- #ifdef DEBUG_print
- if (Verbose)
- fprintf(stderr, "%d leaves and parents for %d clusters, largest cluster = %d\n",nz, *ncluster, ncmax);
+#ifdef DEBUG_print
+ if (Verbose)
+ fprintf(stderr,
+ "%d leaves and parents for %d clusters, largest cluster = %d\n",
+ nz, *ncluster, ncmax);
#endif
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- first = TRUE;
- if (matched[i] == MATCHED) continue;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
- if (first) {
- amax = a[j];
- jamax = ja[j];
- first = FALSE;
- } else {
- if (a[j] > amax){
- amax = a[j];
- jamax = ja[j];
- }
- }
- }
- }
- if (!first){
- matched[jamax] = MATCHED;
- matched[i] = MATCHED;
- (*cluster)[nz++] = i;
- (*cluster)[nz++] = jamax;
- (*clusterp)[++(*ncluster)] = nz;
- }
- }
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ first = TRUE;
+ if (matched[i] == MATCHED)
+ continue;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if (matched[ja[j]] != MATCHED && matched[i] != MATCHED) {
+ if (first) {
+ amax = a[j];
+ jamax = ja[j];
+ first = FALSE;
+ } else {
+ if (a[j] > amax) {
+ amax = a[j];
+ jamax = ja[j];
+ }
+ }
+ }
+ }
+ if (!first) {
+ matched[jamax] = MATCHED;
+ matched[i] = MATCHED;
+ (*cluster)[nz++] = i;
+ (*cluster)[nz++] = jamax;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
+ }
- /* dan yi dian, wu ban */
- for (i = 0; i < m; i++){
- if (matched[i] == i){
- (*cluster)[nz++] = i;
- (*clusterp)[++(*ncluster)] = nz;
- }
- }
+ /* dan yi dian, wu ban */
+ for (i = 0; i < m; i++) {
+ if (matched[i] == i) {
+ (*cluster)[nz++] = i;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
+ }
- FREE(p);
- }
+ free(p);
+ }
- FREE(matched);
+ free(matched);
}
-static void maximal_independent_edge_set_heavest_edge_pernode_supernodes_first(SparseMatrix A, int randomize, int **cluster, int **clusterp, int *ncluster){
- int i, ii, j, *ia, *ja, m, n, *p = NULL;
- real *a, amax = 0;
- int first = TRUE, jamax = 0;
- int *matched, nz, nz0;
- enum {UNMATCHED = -2, MATCHED = -1};
- int nsuper, *super = NULL, *superp = NULL;
-
- assert(A);
- assert(SparseMatrix_known_strucural_symmetric(A));
- ia = A->ia;
- ja = A->ja;
- m = A->m;
- n = A->n;
- assert(n == m);
- *cluster = MALLOC(sizeof(int)*m);
- *clusterp = MALLOC(sizeof(int)*(m+1));
- matched = MALLOC(sizeof(int)*m);
-
- for (i = 0; i < m; i++) matched[i] = i;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
- assert(A->type == MATRIX_TYPE_REAL);
-
- SparseMatrix_decompose_to_supervariables(A, &nsuper, &super, &superp);
-
- *ncluster = 0;
- (*clusterp)[0] = 0;
- nz = 0;
- a = (real*) A->a;
-
- for (i = 0; i < nsuper; i++){
- if (superp[i+1] - superp[i] <= 1) continue;
- nz0 = (*clusterp)[*ncluster];
- for (j = superp[i]; j < superp[i+1]; j++){
- matched[super[j]] = MATCHED;
- (*cluster)[nz++] = super[j];
- if (nz - nz0 >= MAX_CLUSTER_SIZE){
- (*clusterp)[++(*ncluster)] = nz;
- nz0 = nz;
- }
- }
- if (nz > nz0) (*clusterp)[++(*ncluster)] = nz;
- }
-
- if (!randomize){
- for (i = 0; i < m; i++){
- first = TRUE;
- if (matched[i] == MATCHED) continue;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
- if (first) {
- amax = a[j];
- jamax = ja[j];
- first = FALSE;
- } else {
- if (a[j] > amax){
- amax = a[j];
- jamax = ja[j];
- }
- }
- }
- }
- if (!first){
- matched[jamax] = MATCHED;
- matched[i] = MATCHED;
- (*cluster)[nz++] = i;
- (*cluster)[nz++] = jamax;
- (*clusterp)[++(*ncluster)] = nz;
- }
- }
+static void
+maximal_independent_edge_set_heavest_edge_pernode_supernodes_first
+(SparseMatrix * A, int randomize, int **cluster, int **clusterp, int *ncluster)
+{
+ int i, ii, j, *ia, *ja, m, n, *p = NULL;
+ real *a, amax = 0;
+ int first = TRUE, jamax = 0;
+ int *matched, nz, nz0;
+ enum { UNMATCHED = -2, MATCHED = -1 };
+ int nsuper, *super = NULL, *superp = NULL;
- /* dan yi dian, wu ban */
- for (i = 0; i < m; i++){
- if (matched[i] == i){
- (*cluster)[nz++] = i;
- (*clusterp)[++(*ncluster)] = nz;
- }
- }
- assert(nz == n);
-
- } else {
- p = random_permutation(m);
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- first = TRUE;
- if (matched[i] == MATCHED) continue;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
- if (first) {
- amax = a[j];
- jamax = ja[j];
- first = FALSE;
- } else {
- if (a[j] > amax){
- amax = a[j];
- jamax = ja[j];
- }
- }
- }
- }
- if (!first){
- matched[jamax] = MATCHED;
- matched[i] = MATCHED;
- (*cluster)[nz++] = i;
- (*cluster)[nz++] = jamax;
- (*clusterp)[++(*ncluster)] = nz;
- }
+ assert(A);
+ assert(SparseMatrix_known_strucural_symmetric(A));
+ ia = A->ia;
+ ja = A->ja;
+ m = A->m;
+ n = A->n;
+ assert(n == m);
+ *cluster = N_GNEW(m, int);
+ *clusterp = N_GNEW(m + 1, int);
+ matched = N_GNEW(m, int);
+
+ for (i = 0; i < m; i++)
+ matched[i] = i;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+ assert(A->type == MATRIX_TYPE_REAL);
+
+ SparseMatrix_decompose_to_supervariables(A, &nsuper, &super, &superp);
+
+ *ncluster = 0;
+ (*clusterp)[0] = 0;
+ nz = 0;
+ a = (real *) A->a;
+
+ for (i = 0; i < nsuper; i++) {
+ if (superp[i + 1] - superp[i] <= 1)
+ continue;
+ nz0 = (*clusterp)[*ncluster];
+ for (j = superp[i]; j < superp[i + 1]; j++) {
+ matched[super[j]] = MATCHED;
+ (*cluster)[nz++] = super[j];
+ if (nz - nz0 >= MAX_CLUSTER_SIZE) {
+ (*clusterp)[++(*ncluster)] = nz;
+ nz0 = nz;
+ }
+ }
+ if (nz > nz0)
+ (*clusterp)[++(*ncluster)] = nz;
}
- /* dan yi dian, wu ban */
- for (i = 0; i < m; i++){
- if (matched[i] == i){
- (*cluster)[nz++] = i;
- (*clusterp)[++(*ncluster)] = nz;
- }
- }
- FREE(p);
+ if (!randomize) {
+ for (i = 0; i < m; i++) {
+ first = TRUE;
+ if (matched[i] == MATCHED)
+ continue;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if (matched[ja[j]] != MATCHED && matched[i] != MATCHED) {
+ if (first) {
+ amax = a[j];
+ jamax = ja[j];
+ first = FALSE;
+ } else {
+ if (a[j] > amax) {
+ amax = a[j];
+ jamax = ja[j];
+ }
+ }
+ }
+ }
+ if (!first) {
+ matched[jamax] = MATCHED;
+ matched[i] = MATCHED;
+ (*cluster)[nz++] = i;
+ (*cluster)[nz++] = jamax;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
+ }
+
+ /* dan yi dian, wu ban */
+ for (i = 0; i < m; i++) {
+ if (matched[i] == i) {
+ (*cluster)[nz++] = i;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
+ }
+ assert(nz == n);
- }
+ } else {
+ p = random_permutation(m);
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ first = TRUE;
+ if (matched[i] == MATCHED)
+ continue;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if (matched[ja[j]] != MATCHED && matched[i] != MATCHED) {
+ if (first) {
+ amax = a[j];
+ jamax = ja[j];
+ first = FALSE;
+ } else {
+ if (a[j] > amax) {
+ amax = a[j];
+ jamax = ja[j];
+ }
+ }
+ }
+ }
+ if (!first) {
+ matched[jamax] = MATCHED;
+ matched[i] = MATCHED;
+ (*cluster)[nz++] = i;
+ (*cluster)[nz++] = jamax;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
+ }
- FREE(super);
+ /* dan yi dian, wu ban */
+ for (i = 0; i < m; i++) {
+ if (matched[i] == i) {
+ (*cluster)[nz++] = i;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
+ }
+ free(p);
- FREE(superp);
+ }
- FREE(matched);
-}
+ free(super);
+
+ free(superp);
-static int scomp(const void *s1, const void *s2){
- real *ss1, *ss2;
- ss1 = (real*) s1;
- ss2 = (real*) s2;
-
- if ((ss1)[1] > (ss2)[1]){
- return -1;
- } else if ((ss1)[1] < (ss2)[1]){
- return 1;
- }
- return 0;
+ free(matched);
}
-static void maximal_independent_edge_set_heavest_cluster_pernode_leaves_first(SparseMatrix A, int csize,
- int randomize, int **cluster, int **clusterp, int *ncluster){
- int i, ii, j, *ia, *ja, m, n, *p = NULL, q, iv, ncmax;
- real *a;
- int *matched, nz, nz0, nzz,k, nv;
- enum {UNMATCHED = -2, MATCHED = -1};
- real *vlist;
-
- assert(A);
- assert(SparseMatrix_known_strucural_symmetric(A));
- ia = A->ia;
- ja = A->ja;
- m = A->m;
- n = A->n;
- assert(n == m);
- *cluster = MALLOC(sizeof(int)*m);
- *clusterp = MALLOC(sizeof(int)*(m+1));
- matched = MALLOC(sizeof(int)*m);
- vlist = MALLOC(sizeof(real)*2*m);
-
- for (i = 0; i < m; i++) matched[i] = i;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
- assert(A->type == MATRIX_TYPE_REAL);
-
- *ncluster = 0;
- (*clusterp)[0] = 0;
- nz = 0;
- a = (real*) A->a;
-
- p = random_permutation(m);
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- if (matched[i] == MATCHED || node_degree(i) != 1) continue;
- q = ja[ia[i]];
- assert(matched[q] != MATCHED);
- matched[q] = MATCHED;
- (*cluster)[nz++] = q;
- for (j = ia[q]; j < ia[q+1]; j++){
- if (q == ja[j]) continue;
- if (node_degree(ja[j]) == 1){
- matched[ja[j]] = MATCHED;
- (*cluster)[nz++] = ja[j];
- }
- }
- ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]);
- nz0 = (*clusterp)[*ncluster];
- if (nz - nz0 <= MAX_CLUSTER_SIZE){
- (*clusterp)[++(*ncluster)] = nz;
- } else {
- (*clusterp)[++(*ncluster)] = ++nz0;
- nzz = nz0;
- for (k = nz0; k < nz && nzz < nz; k++){
- nzz += MAX_CLUSTER_SIZE - 1;
- nzz = MIN(nz, nzz);
- (*clusterp)[++(*ncluster)] = nzz;
- }
+static int scomp(const void *s1, const void *s2)
+{
+ real *ss1, *ss2;
+ ss1 = (real *) s1;
+ ss2 = (real *) s2;
+
+ if ((ss1)[1] > (ss2)[1]) {
+ return -1;
+ } else if ((ss1)[1] < (ss2)[1]) {
+ return 1;
}
- }
-
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- if (matched[i] == MATCHED) continue;
- nv = 0;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if (matched[ja[j]] != MATCHED && matched[i] != MATCHED){
- vlist[2*nv] = ja[j];
- vlist[2*nv+1] = a[j];
- nv++;
- }
+ return 0;
+}
+
+static void
+maximal_independent_edge_set_heavest_cluster_pernode_leaves_first
+(SparseMatrix * A, int csize, int randomize, int **cluster, int **clusterp,
+int *ncluster)
+{
+ int i, ii, j, *ia, *ja, m, n, *p = NULL, q, iv, ncmax;
+ real *a;
+ int *matched, nz, nz0, nzz, k, nv;
+ enum { UNMATCHED = -2, MATCHED = -1 };
+ real *vlist;
+
+ assert(A);
+ assert(SparseMatrix_known_strucural_symmetric(A));
+ ia = A->ia;
+ ja = A->ja;
+ m = A->m;
+ n = A->n;
+ assert(n == m);
+ *cluster = N_GNEW(m, int);
+ *clusterp = N_GNEW(m + 1, int);
+ matched = N_GNEW(m, int);
+ vlist = N_GNEW(2 * m, real);
+
+ for (i = 0; i < m; i++)
+ matched[i] = i;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+ assert(A->type == MATRIX_TYPE_REAL);
+
+ *ncluster = 0;
+ (*clusterp)[0] = 0;
+ nz = 0;
+ a = (real *) A->a;
+
+ p = random_permutation(m);
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ if (matched[i] == MATCHED || node_degree(i) != 1)
+ continue;
+ q = ja[ia[i]];
+ assert(matched[q] != MATCHED);
+ matched[q] = MATCHED;
+ (*cluster)[nz++] = q;
+ for (j = ia[q]; j < ia[q + 1]; j++) {
+ if (q == ja[j])
+ continue;
+ if (node_degree(ja[j]) == 1) {
+ matched[ja[j]] = MATCHED;
+ (*cluster)[nz++] = ja[j];
+ }
+ }
+ ncmax = MAX(ncmax, nz - (*clusterp)[*ncluster]);
+ nz0 = (*clusterp)[*ncluster];
+ if (nz - nz0 <= MAX_CLUSTER_SIZE) {
+ (*clusterp)[++(*ncluster)] = nz;
+ } else {
+ (*clusterp)[++(*ncluster)] = ++nz0;
+ nzz = nz0;
+ for (k = nz0; k < nz && nzz < nz; k++) {
+ nzz += MAX_CLUSTER_SIZE - 1;
+ nzz = MIN(nz, nzz);
+ (*clusterp)[++(*ncluster)] = nzz;
+ }
+ }
}
- if (nv > 0){
- qsort(vlist, nv, sizeof(real)*2, scomp);
- for (j = 0; j < MIN(csize - 1, nv); j++){
- iv = (int) vlist[2*j];
- matched[iv] = MATCHED;
- (*cluster)[nz++] = iv;
- }
- matched[i] = MATCHED;
- (*cluster)[nz++] = i;
- (*clusterp)[++(*ncluster)] = nz;
+
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ if (matched[i] == MATCHED)
+ continue;
+ nv = 0;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if (matched[ja[j]] != MATCHED && matched[i] != MATCHED) {
+ vlist[2 * nv] = ja[j];
+ vlist[2 * nv + 1] = a[j];
+ nv++;
+ }
+ }
+ if (nv > 0) {
+ qsort(vlist, nv, sizeof(real) * 2, scomp);
+ for (j = 0; j < MIN(csize - 1, nv); j++) {
+ iv = (int) vlist[2 * j];
+ matched[iv] = MATCHED;
+ (*cluster)[nz++] = iv;
+ }
+ matched[i] = MATCHED;
+ (*cluster)[nz++] = i;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
}
- }
-
- /* dan yi dian, wu ban */
- for (i = 0; i < m; i++){
- if (matched[i] == i){
- (*cluster)[nz++] = i;
- (*clusterp)[++(*ncluster)] = nz;
+
+ /* dan yi dian, wu ban */
+ for (i = 0; i < m; i++) {
+ if (matched[i] == i) {
+ (*cluster)[nz++] = i;
+ (*clusterp)[++(*ncluster)] = nz;
+ }
}
- }
- FREE(p);
+ free(p);
- FREE(matched);
+ free(matched);
}
-static void maximal_independent_edge_set_heavest_edge_pernode_scaled(SparseMatrix A, int randomize, int **matching, int *nmatch){
- int i, ii, j, *ia, *ja, m, n, *p = NULL;
- real *a, amax = 0;
- int first = TRUE, jamax = 0;
-
- assert(A);
- assert(SparseMatrix_known_strucural_symmetric(A));
- ia = A->ia;
- ja = A->ja;
- m = A->m;
- n = A->n;
- assert(n == m);
- *matching = MALLOC(sizeof(int)*m);
- for (i = 0; i < m; i++) (*matching)[i] = i;
- *nmatch = n;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
- assert(A->type == MATRIX_TYPE_REAL);
-
- a = (real*) A->a;
- if (!randomize){
- for (i = 0; i < m; i++){
- first = TRUE;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
- if (first) {
- amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]);
- jamax = ja[j];
- first = FALSE;
- } else {
- if (a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]) > amax){
- amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]);
- jamax = ja[j];
- }
- }
- }
- }
- if (!first){
- (*matching)[jamax] = i;
- (*matching)[i] = jamax;
- (*nmatch)--;
- }
- }
- } else {
- p = random_permutation(m);
- for (ii = 0; ii < m; ii++){
- i = p[ii];
- if ((*matching)[i] != i) continue;
- first = TRUE;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i){
- if (first) {
- amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]);
- jamax = ja[j];
- first = FALSE;
- } else {
- if (a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]) > amax){
- amax = a[j]/(ia[i+1]-ia[i])/(ia[ja[j]+1]-ia[ja[j]]);
- jamax = ja[j];
- }
- }
- }
- }
- if (!first){
- (*matching)[jamax] = i;
- (*matching)[i] = jamax;
- (*nmatch)--;
- }
+static void
+maximal_independent_edge_set_heavest_edge_pernode_scaled(SparseMatrix * A,
+ int randomize,
+ int **matching,
+ int *nmatch)
+{
+ int i, ii, j, *ia, *ja, m, n, *p = NULL;
+ real *a, amax = 0;
+ int first = TRUE, jamax = 0;
+
+ assert(A);
+ assert(SparseMatrix_known_strucural_symmetric(A));
+ ia = A->ia;
+ ja = A->ja;
+ m = A->m;
+ n = A->n;
+ assert(n == m);
+ *matching = N_GNEW(m, int);
+ for (i = 0; i < m; i++)
+ (*matching)[i] = i;
+ *nmatch = n;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+ assert(A->type == MATRIX_TYPE_REAL);
+
+ a = (real *) A->a;
+ if (!randomize) {
+ for (i = 0; i < m; i++) {
+ first = TRUE;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i) {
+ if (first) {
+ amax =
+ a[j] / (ia[i + 1] - ia[i]) / (ia[ja[j] + 1] -
+ ia[ja[j]]);
+ jamax = ja[j];
+ first = FALSE;
+ } else {
+ if (a[j] / (ia[i + 1] - ia[i]) /
+ (ia[ja[j] + 1] - ia[ja[j]]) > amax) {
+ amax =
+ a[j] / (ia[i + 1] -
+ ia[i]) / (ia[ja[j] + 1] -
+ ia[ja[j]]);
+ jamax = ja[j];
+ }
+ }
+ }
+ }
+ if (!first) {
+ (*matching)[jamax] = i;
+ (*matching)[i] = jamax;
+ (*nmatch)--;
+ }
+ }
+ } else {
+ p = random_permutation(m);
+ for (ii = 0; ii < m; ii++) {
+ i = p[ii];
+ if ((*matching)[i] != i)
+ continue;
+ first = TRUE;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ if ((*matching)[ja[j]] == ja[j] && (*matching)[i] == i) {
+ if (first) {
+ amax =
+ a[j] / (ia[i + 1] - ia[i]) / (ia[ja[j] + 1] -
+ ia[ja[j]]);
+ jamax = ja[j];
+ first = FALSE;
+ } else {
+ if (a[j] / (ia[i + 1] - ia[i]) /
+ (ia[ja[j] + 1] - ia[ja[j]]) > amax) {
+ amax =
+ a[j] / (ia[i + 1] -
+ ia[i]) / (ia[ja[j] + 1] -
+ ia[ja[j]]);
+ jamax = ja[j];
+ }
+ }
+ }
+ }
+ if (!first) {
+ (*matching)[jamax] = i;
+ (*matching)[i] = jamax;
+ (*nmatch)--;
+ }
+ }
+ free(p);
}
- FREE(p);
- }
}
-static void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, real *node_wgt, real **cnode_wgt,
- SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used){
- int *matching = NULL, nmatch, nc, nzc, n, i;
- int *irn = NULL, *jcn = NULL, *ia = NULL, *ja = NULL;
- real *val = NULL;
- SparseMatrix B = NULL;
- int *vset = NULL, nvset, ncov, j;
- int *cluster, *clusterp, ncluster;
-
- assert(A->m == A->n);
- *cA = NULL;
- *P = NULL;
- *R = NULL;
- n = A->m;
-
- *coarsen_scheme_used = ctrl->coarsen_scheme;
-
- switch (ctrl->coarsen_scheme){
- case COARSEN_HYBRID:
+static void Multilevel_coarsen(SparseMatrix * A, SparseMatrix * *cA,
+ real * node_wgt, real ** cnode_wgt,
+ SparseMatrix * *P, SparseMatrix * *R,
+ Multilevel_control * ctrl,
+ int *coarsen_scheme_used)
+{
+ int *matching = NULL, nmatch, nc, nzc, n, i;
+ int *irn = NULL, *jcn = NULL, *ia = NULL, *ja = NULL;
+ real *val = NULL;
+ SparseMatrix *B = NULL;
+ int *vset = NULL, nvset, ncov, j;
+ int *cluster, *clusterp, ncluster;
+
+ assert(A->m == A->n);
+ *cA = NULL;
+ *P = NULL;
+ *R = NULL;
+ n = A->m;
+
+ *coarsen_scheme_used = ctrl->coarsen_scheme;
+
+ switch (ctrl->coarsen_scheme) {
+ case COARSEN_HYBRID:
#ifdef DEBUG_PRINT
- if (Verbose)
- fprintf(stderr, "hybrid scheme, try COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST first\n");
+ if (Verbose)
+ fprintf(stderr,
+ "hybrid scheme, try COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST first\n");
#endif
- *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST;
- Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
+ *coarsen_scheme_used = ctrl->coarsen_scheme =
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST;
+ Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl,
+ coarsen_scheme_used);
- if (!(*cA)) {
+ if (!(*cA)) {
#ifdef DEBUG_PRINT
- if (Verbose)
- fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST\n");
+ if (Verbose)
+ fprintf(stderr,
+ "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST\n");
#endif
- *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST;
- Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
- }
+ *coarsen_scheme_used = ctrl->coarsen_scheme =
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST;
+ Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl,
+ coarsen_scheme_used);
+ }
- if (!(*cA)) {
+ if (!(*cA)) {
#ifdef DEBUG_PRINT
- if (Verbose)
- fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST\n");
+ if (Verbose)
+ fprintf(stderr,
+ "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST\n");
#endif
- *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST;
- Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
- }
+ *coarsen_scheme_used = ctrl->coarsen_scheme =
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST;
+ Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl,
+ coarsen_scheme_used);
+ }
- if (!(*cA)) {
+ if (!(*cA)) {
#ifdef DEBUG_PRINT
- if (Verbose)
- fprintf(stderr, "switching to COARSEN_INDEPENDENT_VERTEX_SET\n");
+ if (Verbose)
+ fprintf(stderr,
+ "switching to COARSEN_INDEPENDENT_VERTEX_SET\n");
#endif
- *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_VERTEX_SET;
- Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
- }
+ *coarsen_scheme_used = ctrl->coarsen_scheme =
+ COARSEN_INDEPENDENT_VERTEX_SET;
+ Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl,
+ coarsen_scheme_used);
+ }
- if (!(*cA)) {
+ if (!(*cA)) {
#ifdef DEBUG_PRINT
- if (Verbose)
- fprintf(stderr, "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE\n");
+ if (Verbose)
+ fprintf(stderr,
+ "switching to COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE\n");
#endif
- *coarsen_scheme_used = ctrl->coarsen_scheme = COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE;
- Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl, coarsen_scheme_used);
- }
- ctrl->coarsen_scheme = COARSEN_HYBRID;
- break;
- case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST:
- case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST:
- case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST:
- if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST) {
- maximal_independent_edge_set_heavest_edge_pernode_leaves_first(A, ctrl->randomize, &cluster, &clusterp, &ncluster);
- } else if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST) {
- maximal_independent_edge_set_heavest_edge_pernode_supernodes_first(A, ctrl->randomize, &cluster, &clusterp, &ncluster);
- } else {
- maximal_independent_edge_set_heavest_cluster_pernode_leaves_first(A, 4, ctrl->randomize, &cluster, &clusterp, &ncluster);
- }
- assert(ncluster <= n);
- nc = ncluster;
- if (nc > ctrl->min_coarsen_factor*n || nc < ctrl->minsize) {
+ *coarsen_scheme_used = ctrl->coarsen_scheme =
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE;
+ Multilevel_coarsen(A, cA, node_wgt, cnode_wgt, P, R, ctrl,
+ coarsen_scheme_used);
+ }
+ ctrl->coarsen_scheme = COARSEN_HYBRID;
+ break;
+ case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST:
+ case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST:
+ case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST:
+ if (ctrl->coarsen_scheme ==
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST)
+ {
+ maximal_independent_edge_set_heavest_edge_pernode_leaves_first
+ (A, ctrl->randomize, &cluster, &clusterp, &ncluster);
+ } else if (ctrl->coarsen_scheme ==
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST)
+ {
+ maximal_independent_edge_set_heavest_edge_pernode_supernodes_first
+ (A, ctrl->randomize, &cluster, &clusterp, &ncluster);
+ } else {
+ maximal_independent_edge_set_heavest_cluster_pernode_leaves_first
+ (A, 4, ctrl->randomize, &cluster, &clusterp, &ncluster);
+ }
+ assert(ncluster <= n);
+ nc = ncluster;
+ if (nc > ctrl->min_coarsen_factor * n || nc < ctrl->minsize) {
#ifdef DEBUG_PRINT
- if (Verbose)
- fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
+ if (Verbose)
+ fprintf(stderr,
+ "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",
+ nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
#endif
- goto RETURN;
- }
- irn = MALLOC(sizeof(int)*n);
- jcn = MALLOC(sizeof(int)*n);
- val = MALLOC(sizeof(real)*n);
- nzc = 0;
- for (i = 0; i < ncluster; i++){
- for (j = clusterp[i]; j < clusterp[i+1]; j++){
- assert(clusterp[i+1] > clusterp[i]);
- irn[nzc] = cluster[j];
- jcn[nzc] = i;
- val[nzc++] = 1.;
- }
- }
- assert(nzc == n);
- *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL);
- *R = SparseMatrix_transpose(*P);
- B = SparseMatrix_multiply(*R, A);
- *cA = SparseMatrix_multiply(B, *P);
- SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
- *R = SparseMatrix_divide_row_by_degree(*R);
- SparseMatrix_set_symmetric(*cA);
- SparseMatrix_set_pattern_symmetric(*cA);
- *cA = SparseMatrix_remove_diagonal(*cA);
- break;
- case COARSEN_INDEPENDENT_EDGE_SET:
- maximal_independent_edge_set(A, ctrl->randomize, &matching, &nmatch);
- case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE:
- if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE)
- maximal_independent_edge_set_heavest_edge_pernode(A, ctrl->randomize, &matching, &nmatch);
- case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED:
- if (ctrl->coarsen_scheme == COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED)
- maximal_independent_edge_set_heavest_edge_pernode_scaled(A, ctrl->randomize, &matching, &nmatch);
- nc = nmatch;
- if (nc > ctrl->min_coarsen_factor*n || nc < ctrl->minsize) {
+ goto RETURN;
+ }
+ irn = N_GNEW(n, int);
+ jcn = N_GNEW(n, int);
+ val = N_GNEW(n, real);
+ nzc = 0;
+ for (i = 0; i < ncluster; i++) {
+ for (j = clusterp[i]; j < clusterp[i + 1]; j++) {
+ assert(clusterp[i + 1] > clusterp[i]);
+ irn[nzc] = cluster[j];
+ jcn[nzc] = i;
+ val[nzc++] = 1.;
+ }
+ }
+ assert(nzc == n);
+ *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn,
+ (void *) val,
+ MATRIX_TYPE_REAL);
+ *R = SparseMatrix_transpose(*P);
+ B = SparseMatrix_multiply(*R, A);
+ *cA = SparseMatrix_multiply(B, *P);
+ SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
+ *R = SparseMatrix_divide_row_by_degree(*R);
+ SparseMatrix_set_symmetric(*cA);
+ SparseMatrix_set_pattern_symmetric(*cA);
+ *cA = SparseMatrix_remove_diagonal(*cA);
+ break;
+ case COARSEN_INDEPENDENT_EDGE_SET:
+ maximal_independent_edge_set(A, ctrl->randomize, &matching,
+ &nmatch);
+ case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE:
+ if (ctrl->coarsen_scheme ==
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE)
+ maximal_independent_edge_set_heavest_edge_pernode(A,
+ ctrl->
+ randomize,
+ &matching,
+ &nmatch);
+ case COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED:
+ if (ctrl->coarsen_scheme ==
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED)
+ maximal_independent_edge_set_heavest_edge_pernode_scaled(A,
+ ctrl->
+ randomize,
+ &matching,
+ &nmatch);
+ nc = nmatch;
+ if (nc > ctrl->min_coarsen_factor * n || nc < ctrl->minsize) {
#ifdef DEBUG_PRINT
- if (Verbose)
- fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
+ if (Verbose)
+ fprintf(stderr,
+ "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",
+ nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
#endif
- goto RETURN;
- }
- irn = MALLOC(sizeof(int)*n);
- jcn = MALLOC(sizeof(int)*n);
- val = MALLOC(sizeof(real)*n);
- nzc = 0; nc = 0;
- for (i = 0; i < n; i++){
- if (matching[i] >= 0){
- if (matching[i] == i){
- irn[nzc] = i;
- jcn[nzc] = nc;
- val[nzc++] = 1.;
- } else {
- irn[nzc] = i;
- jcn[nzc] = nc;
- val[nzc++] = 1;
- irn[nzc] = matching[i];
- jcn[nzc] = nc;
- val[nzc++] = 1;
- matching[matching[i]] = -1;
- }
- nc++;
- matching[i] = -1;
- }
- }
- assert(nc == nmatch);
- assert(nzc == n);
- *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL);
- *R = SparseMatrix_transpose(*P);
- B = SparseMatrix_multiply(*R, A);
- *cA = SparseMatrix_multiply(B, *P);
- SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
- *R = SparseMatrix_divide_row_by_degree(*R);
- SparseMatrix_set_symmetric(*cA);
- SparseMatrix_set_pattern_symmetric(*cA);
- *cA = SparseMatrix_remove_diagonal(*cA);
- break;
- case COARSEN_INDEPENDENT_VERTEX_SET:
- maximal_independent_vertex_set(A, ctrl->randomize, &vset, &nvset, &nzc);
- ia = A->ia;
- ja = A->ja;
- nc = nvset;
- if (nc > ctrl->min_coarsen_factor*n || nc < ctrl->minsize) {
+ goto RETURN;
+ }
+ irn = N_GNEW(n, int);
+ jcn = N_GNEW(n, int);
+ val = N_GNEW(n, real);
+ nzc = 0;
+ nc = 0;
+ for (i = 0; i < n; i++) {
+ if (matching[i] >= 0) {
+ if (matching[i] == i) {
+ irn[nzc] = i;
+ jcn[nzc] = nc;
+ val[nzc++] = 1.;
+ } else {
+ irn[nzc] = i;
+ jcn[nzc] = nc;
+ val[nzc++] = 1;
+ irn[nzc] = matching[i];
+ jcn[nzc] = nc;
+ val[nzc++] = 1;
+ matching[matching[i]] = -1;
+ }
+ nc++;
+ matching[i] = -1;
+ }
+ }
+ assert(nc == nmatch);
+ assert(nzc == n);
+ *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn,
+ (void *) val,
+ MATRIX_TYPE_REAL);
+ *R = SparseMatrix_transpose(*P);
+ B = SparseMatrix_multiply(*R, A);
+ *cA = SparseMatrix_multiply(B, *P);
+ SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
+ *R = SparseMatrix_divide_row_by_degree(*R);
+ SparseMatrix_set_symmetric(*cA);
+ SparseMatrix_set_pattern_symmetric(*cA);
+ *cA = SparseMatrix_remove_diagonal(*cA);
+ break;
+ case COARSEN_INDEPENDENT_VERTEX_SET:
+ maximal_independent_vertex_set(A, ctrl->randomize, &vset, &nvset,
+ &nzc);
+ ia = A->ia;
+ ja = A->ja;
+ nc = nvset;
+ if (nc > ctrl->min_coarsen_factor * n || nc < ctrl->minsize) {
#ifdef DEBUG_PRINT
- if (Verbose)
- fprintf(stderr, "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
+ if (Verbose)
+ fprintf(stderr,
+ "nc = %d, nf = %d, minsz = %d, coarsen_factor = %f coarsening stops\n",
+ nc, n, ctrl->minsize, ctrl->min_coarsen_factor);
#endif
- goto RETURN;
- }
- irn = MALLOC(sizeof(int)*nzc);
- jcn = MALLOC(sizeof(int)*nzc);
- val = MALLOC(sizeof(real)*nzc);
- nzc = 0;
- for (i = 0; i < n; i++){
- if (vset[i] == MAX_IND_VTX_SET_U){
- ncov = 0;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (vset[ja[j]] >= MAX_IND_VTX_SET_C){
- ncov++;
- }
- }
- assert(ncov > 0);
- for (j = ia[i]; j < ia[i+1]; j++){
- if (vset[ja[j]] >= MAX_IND_VTX_SET_C){
- irn[nzc] = i;
- jcn[nzc] = vset[ja[j]];
- val[nzc++] = 1./(double) ncov;
- }
- }
- } else {
- assert(vset[i] >= MAX_IND_VTX_SET_C);
- irn[nzc] = i;
- jcn[nzc] = vset[i];
- val[nzc++] = 1.;
- }
- }
+ goto RETURN;
+ }
+ irn = N_GNEW(nzc, int);
+ jcn = N_GNEW(nzc, int);
+ val = N_GNEW(nzc, real);
+ nzc = 0;
+ for (i = 0; i < n; i++) {
+ if (vset[i] == MAX_IND_VTX_SET_U) {
+ ncov = 0;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (vset[ja[j]] >= MAX_IND_VTX_SET_C) {
+ ncov++;
+ }
+ }
+ assert(ncov > 0);
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (vset[ja[j]] >= MAX_IND_VTX_SET_C) {
+ irn[nzc] = i;
+ jcn[nzc] = vset[ja[j]];
+ val[nzc++] = 1. / (double) ncov;
+ }
+ }
+ } else {
+ assert(vset[i] >= MAX_IND_VTX_SET_C);
+ irn[nzc] = i;
+ jcn[nzc] = vset[i];
+ val[nzc++] = 1.;
+ }
+ }
- *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL);
- *R = SparseMatrix_transpose(*P);
- B = SparseMatrix_multiply(*R, A);
- *cA = SparseMatrix_multiply(B, *P);
- SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
- SparseMatrix_set_symmetric(*cA);
- SparseMatrix_set_pattern_symmetric(*cA);
- *cA = SparseMatrix_remove_diagonal(*cA);
- break;
- default:
- goto RETURN;
- }
- RETURN:
- if (matching) FREE(matching);
- if (vset) FREE(vset);
- if (irn) FREE(irn);
- if (jcn) FREE(jcn);
- if (val) FREE(val);
- if (B) SparseMatrix_delete(B);
+ *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn,
+ (void *) val,
+ MATRIX_TYPE_REAL);
+ *R = SparseMatrix_transpose(*P);
+ B = SparseMatrix_multiply(*R, A);
+ *cA = SparseMatrix_multiply(B, *P);
+ SparseMatrix_multiply_vector(*R, node_wgt, cnode_wgt, FALSE);
+ SparseMatrix_set_symmetric(*cA);
+ SparseMatrix_set_pattern_symmetric(*cA);
+ *cA = SparseMatrix_remove_diagonal(*cA);
+ break;
+ default:
+ goto RETURN;
+ }
+ RETURN:
+ if (matching)
+ free(matching);
+ if (vset)
+ free(vset);
+ if (irn)
+ free(irn);
+ if (jcn)
+ free(jcn);
+ if (val)
+ free(val);
+ if (B)
+ SparseMatrix_delete(B);
}
-void print_padding(int n){
- int i;
- for (i = 0; i < n; i++) fputs (" ", stderr);
+void print_padding(int n)
+{
+ int i;
+ for (i = 0; i < n; i++)
+ fputs(" ", stderr);
}
-static Multilevel Multilevel_establish(Multilevel grid, Multilevel_control ctrl){
- Multilevel cgrid;
- int coarsen_scheme_used;
- real *cnode_weights = NULL;
- SparseMatrix P, R, A, cA;
+static Multilevel *Multilevel_establish(Multilevel * grid,
+ Multilevel_control * ctrl)
+{
+ Multilevel *cgrid;
+ int coarsen_scheme_used;
+ real *cnode_weights = NULL;
+ SparseMatrix *P, *R, *A, *cA;
#ifdef DEBUG_PRINT
- if (Verbose) {
- print_padding(grid->level);
- fprintf(stderr, "level -- %d, n = %d, nz = %d nz/n = %f\n", grid->level, grid->n, grid->A->nz, grid->A->nz/(double) grid->n);
- }
+ if (Verbose) {
+ print_padding(grid->level);
+ fprintf(stderr, "level -- %d, n = %d, nz = %d nz/n = %f\n",
+ grid->level, grid->n, grid->A->nz,
+ grid->A->nz / (double) grid->n);
+ }
#endif
- A = grid->A;
- if (grid->level >= ctrl->maxlevel - 1) {
+ A = grid->A;
+ if (grid->level >= ctrl->maxlevel - 1) {
#ifdef DEBUG_PRINT
- if (Verbose) {
- print_padding(grid->level);
- fprintf(stderr, " maxlevel reached, coarsening stops\n");
- }
+ if (Verbose) {
+ print_padding(grid->level);
+ fprintf(stderr, " maxlevel reached, coarsening stops\n");
+ }
#endif
+ return grid;
+ }
+ Multilevel_coarsen(A, &cA, grid->node_weights, &cnode_weights, &P, &R,
+ ctrl, &coarsen_scheme_used);
+ if (!cA)
+ return grid;
+
+ cgrid = Multilevel_init(cA, cnode_weights);
+ grid->next = cgrid;
+ cgrid->coarsen_scheme_used = coarsen_scheme_used;
+ cgrid->level = grid->level + 1;
+ cgrid->n = cA->m;
+ cgrid->A = cA;
+ cgrid->P = P;
+ grid->R = R;
+ cgrid->prev = grid;
+ cgrid = Multilevel_establish(cgrid, ctrl);
return grid;
- }
- Multilevel_coarsen(A, &cA, grid->node_weights, &cnode_weights, &P, &R, ctrl, &coarsen_scheme_used);
- if (!cA) return grid;
-
- cgrid = Multilevel_init(cA, cnode_weights);
- grid->next = cgrid;
- cgrid->coarsen_scheme_used = coarsen_scheme_used;
- cgrid->level = grid->level + 1;
- cgrid->n = cA->m;
- cgrid->A = cA;
- cgrid->P = P;
- grid->R = R;
- cgrid->prev = grid;
- cgrid = Multilevel_establish(cgrid, ctrl);
- return grid;
-
-}
-Multilevel Multilevel_new(SparseMatrix A0, real *node_weights, Multilevel_control ctrl){
- Multilevel grid;
- SparseMatrix A = A0;
-
- if (!SparseMatrix_is_symmetric(A, FALSE) || A->type != MATRIX_TYPE_REAL){
- A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
- }
- grid = Multilevel_init(A, node_weights);
- grid = Multilevel_establish(grid, ctrl);
- if (A != A0) grid->delete_top_level_A = TRUE;/* be sure to clean up later */
- return grid;
}
+Multilevel *Multilevel_new(SparseMatrix * A0, real * node_weights,
+ Multilevel_control * ctrl)
+{
+ Multilevel *grid;
+ SparseMatrix *A = A0;
-Multilevel Multilevel_get_coarsest(Multilevel grid){
- while (grid->next){
- grid = grid->next;
- }
- return grid;
+ if (!SparseMatrix_is_symmetric(A, FALSE)
+ || A->type != MATRIX_TYPE_REAL) {
+ A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
+ }
+ grid = Multilevel_init(A, node_weights);
+ grid = Multilevel_establish(grid, ctrl);
+ if (A != A0)
+ grid->delete_top_level_A = TRUE; /* be sure to clean up later */
+ return grid;
}
+
+Multilevel *Multilevel_get_coarsest(Multilevel * grid)
+{
+ while (grid->next) {
+ grid = grid->next;
+ }
+ return grid;
+}
+/* vim:set shiftwidth=4 ts=8: */
+
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
#ifndef MULTILEVEL_H
#define MULTILEVEL_H
#include "SparseMatrix.h"
-typedef struct Multilevel_struct *Multilevel;
-
-struct Multilevel_struct {
- int level;/* 0, 1, ... */
- int n;
- SparseMatrix A;
- SparseMatrix P;
- SparseMatrix R;
- real *node_weights;
- Multilevel next;
- Multilevel prev;
- int delete_top_level_A;
- int coarsen_scheme_used;/* to get from previous level to here */
-};
+typedef struct Multilevel_s {
+ int level; /* 0, 1, ... */
+ int n;
+ SparseMatrix *A;
+ SparseMatrix *P;
+ SparseMatrix *R;
+ real *node_weights;
+ struct Multilevel_s *next;
+ struct Multilevel_s *prev;
+ int delete_top_level_A;
+ int coarsen_scheme_used; /* to get from previous level to here */
+} Multilevel;
-enum {MAX_IND_VTX_SET_0 = -100, MAX_IND_VTX_SET_U = -1, MAX_IND_VTX_SET_C = 0};
+enum { MAX_IND_VTX_SET_0 = -100, MAX_IND_VTX_SET_U =
+ -1, MAX_IND_VTX_SET_C = 0 };
-enum {MAX_CLUSTER_SIZE = 4};
+enum { MAX_CLUSTER_SIZE = 4 };
-enum {EDGE_BASED_STA, COARSEN_INDEPENDENT_EDGE_SET, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED, COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST, EDGE_BASED_STO, VERTEX_BASED_STA, COARSEN_INDEPENDENT_VERTEX_SET, VERTEX_BASED_STO, COARSEN_HYBRID};
+enum { EDGE_BASED_STA, COARSEN_INDEPENDENT_EDGE_SET,
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE,
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_LEAVES_FIRST,
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_SUPERNODES_FIRST,
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_EDGE_PERNODE_DEGREE_SCALED,
+ COARSEN_INDEPENDENT_EDGE_SET_HEAVEST_CLUSTER_PERNODE_LEAVES_FIRST,
+ EDGE_BASED_STO, VERTEX_BASED_STA, COARSEN_INDEPENDENT_VERTEX_SET,
+ VERTEX_BASED_STO, COARSEN_HYBRID };
-struct Multilevel_control_struct {
- int minsize;
- real min_coarsen_factor;
- int maxlevel;
- int randomize;
- int coarsen_scheme;
-};
+typedef struct {
+ int minsize;
+ real min_coarsen_factor;
+ int maxlevel;
+ int randomize;
+ int coarsen_scheme;
+} Multilevel_control;
-typedef struct Multilevel_control_struct *Multilevel_control;
-Multilevel_control Multilevel_control_new();
+Multilevel_control *Multilevel_control_new();
-void Multilevel_control_delete(Multilevel_control ctrl);
+void Multilevel_control_delete(Multilevel_control * ctrl);
-void Multilevel_delete(Multilevel grid);
+void Multilevel_delete(Multilevel * grid);
-Multilevel Multilevel_new(SparseMatrix A, real *node_weights, Multilevel_control ctrl);
+Multilevel *Multilevel_new(SparseMatrix * A, real * node_weights,
+ Multilevel_control * ctrl);
-Multilevel Multilevel_get_coarsest(Multilevel grid);
+Multilevel *Multilevel_get_coarsest(Multilevel * grid);
void print_padding(int n);
-#include "general.h"
-#include "LinkedList.h"
-#include "QuadTree.h"
-
-struct node_data_struct {
- real node_weight;
- real *coord;
- real id;
-};
-
-typedef struct node_data_struct *node_data;
+/* vim:set shiftwidth=4 ts=8: */
+
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
-real point_distance(real *p1, real *p2, int dim){
- int i;
- real dist;
- dist = 0;
- for (i = 0; i < dim; i++) dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
- return sqrt(dist);
+#include "QuadTree.h"
+#include "memory.h"
+#include "arith.h"
+#include <math.h>
+#include <assert.h>
+
+typedef struct {
+ real node_weight;
+ real *coord;
+ real id;
+} node_data;
+
+real point_distance(real * p1, real * p2, int dim)
+{
+ int i;
+ real dist;
+ dist = 0;
+ for (i = 0; i < dim; i++)
+ dist += (p1[i] - p2[i]) * (p1[i] - p2[i]);
+ return sqrt(dist);
}
-static node_data node_data_new(int dim, real weight, real *coord, int id){
- node_data nd;
- int i;
- nd = MALLOC(sizeof(struct node_data_struct));
- nd->node_weight = weight;
- nd->coord = MALLOC(sizeof(real)*dim);
- nd->id = id;
- for (i = 0; i < dim; i++) nd->coord[i] = coord[i];
- return nd;
+static node_data *node_data_new(int dim, real weight, real * coord, int id)
+{
+ node_data *nd;
+ int i;
+ nd = GNEW(node_data);
+ nd->node_weight = weight;
+ nd->coord = N_GNEW(dim, real);
+ nd->id = id;
+ for (i = 0; i < dim; i++)
+ nd->coord[i] = coord[i];
+ return nd;
}
-void node_data_delete(void *d){
- node_data nd = (node_data) d;
- FREE(nd->coord);
- FREE(nd);
+void node_data_delete(void *d)
+{
+ node_data *nd = (node_data *) d;
+ free(nd->coord);
+ free(nd);
}
-real node_data_get_weight(void *d){
- node_data nd = (node_data) d;
- return nd->node_weight;
+real node_data_get_weight(void *d)
+{
+ node_data *nd = (node_data *) d;
+ return nd->node_weight;
}
-real* node_data_get_coord(void *d){
- node_data nd = (node_data) d;
- return nd->coord;
+real *node_data_get_coord(void *d)
+{
+ node_data *nd = (node_data *) d;
+ return nd->coord;
}
-int node_data_get_id(void *d){
- node_data nd = (node_data) d;
- return nd->id;
+int node_data_get_id(void *d)
+{
+ node_data *nd = (node_data *) d;
+ return nd->id;
}
-void check_or_realloc_arrays(int dim, int *nsuper, int *nsupermax, real **center, real **supernode_wgts, real **distances){
-
- if (*nsuper >= *nsupermax) {
- *nsupermax = *nsuper + MAX(10, (int) 0.2*(*nsuper));
- *center = REALLOC(*center, sizeof(real)*(*nsupermax)*dim);
- *supernode_wgts = REALLOC(*supernode_wgts, sizeof(real)*(*nsupermax));
- *distances = REALLOC(*distances, sizeof(real)*(*nsupermax));
- }
+int check_or_realloc_arrays(int dim, int nsuper, int nsupermax,
+ real ** center, real ** supernode_wgts,
+ real ** distances)
+{
+
+ if (nsuper >= nsupermax) {
+ nsupermax = nsuper + MAX(10, (int) 0.2 * (nsuper));
+ *center = RALLOC(nsupermax * dim, *center, real);
+ *supernode_wgts = RALLOC(nsupermax, *supernode_wgts, real);
+ *distances = RALLOC(nsupermax, *distances, real);
+ }
+ return nsupermax;
}
-void QuadTree_get_supernodes_internal(QuadTree qt, real bh, real *point, int nodeid, int *nsuper, int *nsupermax, real **center, real **supernode_wgts, real **distances, real *counts, int *flag){
- SingleLinkedList l;
- real *coord, dist;
- int dim = qt->dim, i;
-
- (*counts)++;
- l = qt->l;
- if (l){
- while (l){
- check_or_realloc_arrays(dim, nsuper, nsupermax, center, supernode_wgts, distances);
- if (node_data_get_id(SingleLinkedList_get_data(l)) != nodeid){
- coord = node_data_get_coord(SingleLinkedList_get_data(l));
- for (i = 0; i < dim; i++){
- (*center)[dim*(*nsuper)+i] = coord[i];
+void QuadTree_get_supernodes_internal(QuadTree * qt, real bh, real * point,
+ int nodeid, int *nsuper,
+ int *nsupermax, real ** center,
+ real ** supernode_wgts,
+ real ** distances, real * counts,
+ int *flag)
+{
+ SingleLinkedList *l;
+ real *coord, dist;
+ int dim = qt->dim, i;
+
+ (*counts)++;
+ l = qt->l;
+ if (l) {
+ while (l) {
+ *nsupermax =
+ check_or_realloc_arrays(dim, *nsuper, *nsupermax, center,
+ supernode_wgts, distances);
+ if (node_data_get_id(SingleLinkedList_get_data(l)) != nodeid) {
+ coord = node_data_get_coord(SingleLinkedList_get_data(l));
+ for (i = 0; i < dim; i++) {
+ (*center)[dim * (*nsuper) + i] = coord[i];
+ }
+ (*supernode_wgts)[*nsuper] =
+ node_data_get_weight(SingleLinkedList_get_data(l));
+ (*distances)[*nsuper] = point_distance(point, coord, dim);
+ (*nsuper)++;
+ }
+ l = SingleLinkedList_get_next(l);
}
- (*supernode_wgts)[*nsuper] = node_data_get_weight(SingleLinkedList_get_data(l));
- (*distances)[*nsuper] = point_distance(point, coord, dim);
- (*nsuper)++;
- }
- l = SingleLinkedList_get_next(l);
}
- }
-
- if (qt->qts){
- dist = point_distance(qt->center, point, dim);
- if (qt->width < bh*dist){
- check_or_realloc_arrays(dim, nsuper, nsupermax, center, supernode_wgts, distances);
- for (i = 0; i < dim; i++){
- (*center)[dim*(*nsuper)+i] = qt->average[i];
- }
- (*supernode_wgts)[*nsuper] = qt->total_weight;
- (*distances)[*nsuper] = point_distance(qt->average, point, dim);
- (*nsuper)++;
- } else {
- for (i = 0; i < 1<<dim; i++){
- QuadTree_get_supernodes_internal(qt->qts[i], bh, point, nodeid, nsuper, nsupermax, center,
- supernode_wgts, distances, counts, flag);
- }
+
+ if (qt->qts) {
+ dist = point_distance(qt->center, point, dim);
+ if (qt->width < bh * dist) {
+ *nsupermax =
+ check_or_realloc_arrays(dim, *nsuper, *nsupermax, center,
+ supernode_wgts, distances);
+ for (i = 0; i < dim; i++) {
+ (*center)[dim * (*nsuper) + i] = qt->average[i];
+ }
+ (*supernode_wgts)[*nsuper] = qt->total_weight;
+ (*distances)[*nsuper] =
+ point_distance(qt->average, point, dim);
+ (*nsuper)++;
+ } else {
+ for (i = 0; i < 1 << dim; i++) {
+ QuadTree_get_supernodes_internal(qt->qts[i], bh, point,
+ nodeid, nsuper, nsupermax,
+ center, supernode_wgts,
+ distances, counts, flag);
+ }
+ }
}
- }
}
-void QuadTree_get_supernodes(QuadTree qt, real bh, real *point, int nodeid, int *nsuper,
- int *nsupermax, real **center, real **supernode_wgts, real **distances, double *counts, int *flag){
- int dim = qt->dim;
-
- (*counts) = 0;
-
- *nsuper = 0;
-
- *flag = 0;
- *nsupermax = 10;
- if (!*center) *center = MALLOC(sizeof(real)*(*nsupermax)*dim);
- if (!*supernode_wgts) *supernode_wgts = MALLOC(sizeof(real)*(*nsupermax));
- if (!*distances) *distances = MALLOC(sizeof(real)*(*nsupermax));
- QuadTree_get_supernodes_internal(qt, bh, point, nodeid, nsuper, nsupermax, center, supernode_wgts, distances, counts, flag);
+void QuadTree_get_supernodes(QuadTree * qt, real bh, real * point,
+ int nodeid, int *nsuper, int *nsupermax,
+ real ** center, real ** supernode_wgts,
+ real ** distances, double *counts, int *flag)
+{
+ int dim = qt->dim;
+
+ (*counts) = 0;
+
+ *nsuper = 0;
+
+ *flag = 0;
+ *nsupermax = 10;
+ if (!*center)
+ *center = N_GNEW((*nsupermax) * dim, real);
+ if (!*supernode_wgts)
+ *supernode_wgts = N_GNEW(*nsupermax, real);
+ if (!*distances)
+ *distances = N_GNEW(*nsupermax, real);
+ QuadTree_get_supernodes_internal(qt, bh, point, nodeid, nsuper,
+ nsupermax, center, supernode_wgts,
+ distances, counts, flag);
}
-QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord, real *weight){
- /* form a new QuadTree data structure from a list of coordinates of n points
- coord: of length n*dim, point i sits at [i*dim, i*dim+dim - 1]
- weight: node weight of lentgth n. If NULL, unit weight assumed.
- */
- real *xmin, *xmax, *center, width;
- QuadTree qt = NULL;
- int i, k;
-
- xmin = MALLOC(sizeof(real)*dim);
- xmax = MALLOC(sizeof(real)*dim);
- center = MALLOC(sizeof(real)*dim);
- if (!xmin || !xmax || !center) return NULL;
-
- for (i = 0; i < dim; i++) xmin[i] = coord[i];
- for (i = 0; i < dim; i++) xmax[i] = coord[i];
-
- for (i = 1; i < n; i++){
- for (k = 0; k < dim; k++){
- xmin[k] = MIN(xmin[k], coord[i*dim+k]);
- xmax[k] = MAX(xmax[k], coord[i*dim+k]);
+QuadTree *QuadTree_new_from_point_list(int dim, int n, int max_level,
+ real * coord, real * weight)
+{
+ /* form a new QuadTree data structure from a list of coordinates of n points
+ coord: of length n*dim, point i sits at [i*dim, i*dim+dim - 1]
+ weight: node weight of lentgth n. If NULL, unit weight assumed.
+ */
+ real *xmin, *xmax, *center, width;
+ QuadTree *qt = NULL;
+ int i, k;
+
+ xmin = N_GNEW(dim, real);
+ xmax = N_GNEW(dim, real);
+ center = N_GNEW(dim, real);
+ if (!xmin || !xmax || !center)
+ return NULL;
+
+ for (i = 0; i < dim; i++)
+ xmin[i] = coord[i];
+ for (i = 0; i < dim; i++)
+ xmax[i] = coord[i];
+
+ for (i = 1; i < n; i++) {
+ for (k = 0; k < dim; k++) {
+ xmin[k] = MIN(xmin[k], coord[i * dim + k]);
+ xmax[k] = MAX(xmax[k], coord[i * dim + k]);
+ }
}
- }
-
- width = xmax[0] - xmin[0];
- for (i = 0; i < dim; i++) {
- center[i] = (xmin[i] + xmax[i])*0.5;
- width = MAX(width, xmax[i] - xmin[i]);
- }
- width *= 0.52;
- qt = QuadTree_new(dim, center, width, max_level);
-
- if (weight){
- for (i = 0; i < n; i++){
- qt = QuadTree_add(qt, &(coord[i*dim]), weight[i], i);
+
+ width = xmax[0] - xmin[0];
+ for (i = 0; i < dim; i++) {
+ center[i] = (xmin[i] + xmax[i]) * 0.5;
+ width = MAX(width, xmax[i] - xmin[i]);
}
- } else {
- for (i = 0; i < n; i++){
- qt = QuadTree_add(qt, &(coord[i*dim]), 1, i);
+ width *= 0.52;
+ qt = QuadTree_new(dim, center, width, max_level);
+
+ if (weight) {
+ for (i = 0; i < n; i++) {
+ qt = QuadTree_add(qt, &(coord[i * dim]), weight[i], i);
+ }
+ } else {
+ for (i = 0; i < n; i++) {
+ qt = QuadTree_add(qt, &(coord[i * dim]), 1, i);
+ }
}
- }
- FREE(xmin);
- FREE(xmax);
- FREE(center);
- return qt;
+ free(xmin);
+ free(xmax);
+ free(center);
+ return qt;
}
-QuadTree QuadTree_new(int dim, real *center, real width, int max_level){
- QuadTree q;
- int i;
- q = MALLOC(sizeof(struct QuadTree_struct));
- q->dim = dim;
- q->n = 0;
- q->center = MALLOC(sizeof(real)*dim);
- for (i = 0; i < dim; i++) q->center[i] = center[i];
- assert(width > 0);
- q->width = width;
- q->total_weight = 0;
- q->average = NULL;
- q->qts = NULL;
- q->l = NULL;
- q->max_level = max_level;
- return q;
+QuadTree *QuadTree_new(int dim, real * center, real width, int max_level)
+{
+ QuadTree *q;
+ int i;
+ q = GNEW(QuadTree);
+ q->dim = dim;
+ q->n = 0;
+ q->center = N_GNEW(dim, real);
+ for (i = 0; i < dim; i++)
+ q->center[i] = center[i];
+ assert(width > 0);
+ q->width = width;
+ q->total_weight = 0;
+ q->average = NULL;
+ q->qts = NULL;
+ q->l = NULL;
+ q->max_level = max_level;
+ return q;
}
-void QuadTree_delete(QuadTree q){
- int i, dim = q->dim;
- if (!q) return;
+void QuadTree_delete(QuadTree * q)
+{
+ int i, dim = q->dim;
+ if (!q)
+ return;
+
+ free(q->center);
+ free(q->average);
+ if (q->qts) {
+ for (i = 0; i < 1 << dim; i++) {
+ QuadTree_delete(q->qts[i]);
+ }
+ free(q->qts);
+ }
+ SingleLinkedList_delete(q->l, node_data_delete);
+ free(q);
+}
- FREE(q->center);
- FREE(q->average);
- if (q->qts){
- for (i = 0; i < 1<<dim; i++){
- QuadTree_delete(q->qts[i]);
+static int QuadTree_get_quadrant(int dim, real * center, real * coord)
+{
+ /* find the quadrant that a point of coordinates coord is going into with reference to the center.
+ if coord - center == {+,-,+,+} = {1,0,1,1}, then it will sit in the i-quadrant where
+ i's binary representation is 1011 (that is, decimal 11).
+ */
+ int d = 0, i;
+
+ for (i = dim - 1; i >= 0; i--) {
+ if (coord[i] - center[i] < 0) {
+ d = 2 * d;
+ } else {
+ d = 2 * d + 1;
+ }
}
- FREE(q->qts);
- }
- SingleLinkedList_delete(q->l, node_data_delete);
- FREE(q);
+ return d;
}
-static int QuadTree_get_quadrant(int dim, real *center, real *coord){
- /* find the quadrant that a point of coordinates coord is going into with reference to the center.
- if coord - center == {+,-,+,+} = {1,0,1,1}, then it will sit in the i-quadrant where
- i's binary representation is 1011 (that is, decimal 11).
- */
- int d = 0, i;
+static QuadTree *QuadTree_add_internal(QuadTree * q, real * coord,
+ real weight, int id, int level)
+{
+ int i, dim = q->dim, ii, k;
+ node_data *nd = NULL;
+ real *center, width;
+ int max_level = q->max_level;
+
+ /* Make sure that coord is within bounding box */
+ for (i = 0; i < q->dim; i++) {
+ if (coord[i] < q->center[i] - q->width
+ || coord[i] > q->center[i] + q->width)
+ return NULL;
+ }
- for (i = dim - 1; i >= 0; i--){
- if (coord[i] - center[i] < 0){
- d = 2*d;
+ if (q->n == 0) {
+ /* if this node is empty right now */
+ q->n = 1;
+ q->total_weight = weight;
+ q->average = N_GNEW(dim, real);
+ for (i = 0; i < q->dim; i++)
+ q->average[i] = coord[i];
+ nd = node_data_new(q->dim, weight, coord, id);
+ assert(!(q->l));
+ q->l = SingleLinkedList_new(nd);
+ } else if (level < max_level) {
+ /* otherwise open up into 2^dim quadtrees unless the level is too high */
+ q->total_weight += weight;
+ for (i = 0; i < q->dim; i++)
+ q->average[i] =
+ ((q->average[i]) * q->n + coord[i]) / (q->n + 1);
+ if (!q->qts) {
+ q->qts = N_GNEW(1 << dim, QuadTree *);
+ for (i = 0; i < 1 << dim; i++) {
+ width = (q->width) / 2;
+ q->qts[i] =
+ QuadTree_new(q->dim, q->center, width, max_level);
+ center = (q->qts[i])->center; /* right now this has the center for the parent */
+ ii = i;
+ for (k = 0; k < dim; k++) { /* decompose child id into binary, if {1,0}, say, then
+ add {width/2, -width/2} to the parents' center
+ to get the child's center. */
+ if (ii % 2 == 0) {
+ center[k] -= width;
+ } else {
+ center[k] += width;
+ }
+ ii = (ii - ii % 2) / 2;
+ }
+ }
+ }
+
+ /* done adding new quadtree, now add points to them */
+ /* insert the old node (if exist) and the current node into the appropriate child quadtree */
+ ii = QuadTree_get_quadrant(dim, q->center, coord);
+ assert(ii < 1 << dim && ii >= 0);
+ q->qts[ii] =
+ QuadTree_add_internal(q->qts[ii], coord, weight, id,
+ level + 1);
+ assert(q->qts[ii]);
+
+ if (q->l) {
+ assert(q->n == 1);
+ coord = node_data_get_coord(SingleLinkedList_get_data(q->l));
+ weight = node_data_get_weight(SingleLinkedList_get_data(q->l));
+ ii = QuadTree_get_quadrant(dim, q->center, coord);
+ assert(ii < 1 << dim && ii >= 0);
+ q->qts[ii] =
+ QuadTree_add_internal(q->qts[ii], coord, weight, id,
+ level + 1);
+ assert(q->qts[ii]);
+
+ /* delete the old node data on parent */
+ SingleLinkedList_delete(q->l, node_data_delete);
+ q->l = NULL;
+ }
+
+ (q->n)++;
} else {
- d = 2*d+1;
+ assert(!(q->qts));
+ /* level is too high, append data in the linked list */
+ (q->n)++;
+ q->total_weight += weight;
+ for (i = 0; i < q->dim; i++)
+ q->average[i] =
+ ((q->average[i]) * q->n + coord[i]) / (q->n + 1);
+ nd = node_data_new(q->dim, weight, coord, id);
+ assert(q->l);
+ q->l = SingleLinkedList_prepend(q->l, nd);
}
- }
- return d;
+ return q;
}
-static QuadTree QuadTree_add_internal(QuadTree q, real *coord, real weight, int id, int level){
- int i, dim = q->dim, ii, k;
- node_data nd = NULL;
- real *center, width;
- int max_level = q->max_level;
-
- /* Make sure that coord is within bounding box */
- for (i = 0; i < q->dim; i++) {
- if (coord[i] < q->center[i] - q->width || coord[i] > q->center[i] + q->width) return NULL;
- }
-
- if (q->n == 0){
- /* if this node is empty right now */
- q->n = 1;
- q->total_weight = weight;
- q->average = MALLOC(sizeof(real)*dim);
- for (i = 0; i < q->dim; i++) q->average[i] = coord[i];
- nd = node_data_new(q->dim, weight, coord, id);
- assert(!(q->l));
- q->l = SingleLinkedList_new(nd);
- } else if (level < max_level){
- /* otherwise open up into 2^dim quadtrees unless the level is too high */
- q->total_weight += weight;
- for (i = 0; i < q->dim; i++) q->average[i] = ((q->average[i])*q->n + coord[i])/(q->n + 1);
- if (!q->qts){
- q->qts = MALLOC(sizeof(QuadTree)*(1<<dim));
- for (i = 0; i < 1<<dim; i++){
- width = (q->width)/2;
- q->qts[i] = QuadTree_new(q->dim, q->center, width, max_level);
- center = (q->qts[i])->center;/* right now this has the center for the parent */
- ii = i;
- for (k = 0; k < dim; k++){/* decompose child id into binary, if {1,0}, say, then
- add {width/2, -width/2} to the parents' center
- to get the child's center. */
- if (ii%2 == 0){
- center[k] -= width;
- } else {
- center[k] += width;
- }
- ii = (ii - ii%2)/2;
- }
- }
- }/* done adding new quadtree, now add points to them */
-
- /* insert the old node (if exist) and the current node into the appropriate child quadtree */
- ii = QuadTree_get_quadrant(dim, q->center, coord);
- assert(ii < 1<<dim && ii >= 0);
- q->qts[ii] = QuadTree_add_internal(q->qts[ii], coord, weight, id, level + 1);
- assert(q->qts[ii]);
-
- if (q->l){
- assert(q->n == 1);
- coord = node_data_get_coord(SingleLinkedList_get_data(q->l));
- weight = node_data_get_weight(SingleLinkedList_get_data(q->l));
- ii = QuadTree_get_quadrant(dim, q->center, coord);
- assert(ii < 1<<dim && ii >= 0);
- q->qts[ii] = QuadTree_add_internal(q->qts[ii], coord, weight, id, level + 1);
- assert(q->qts[ii]);
-
- /* delete the old node data on parent */
- SingleLinkedList_delete(q->l, node_data_delete);
- q->l = NULL;
- }
- (q->n)++;
- } else {
- assert(!(q->qts));
- /* level is too high, append data in the linked list */
- (q->n)++;
- q->total_weight += weight;
- for (i = 0; i < q->dim; i++) q->average[i] = ((q->average[i])*q->n + coord[i])/(q->n + 1);
- nd = node_data_new(q->dim, weight, coord, id);
- assert(q->l);
- q->l = SingleLinkedList_prepend(q->l, nd);
- }
- return q;
+QuadTree *QuadTree_add(QuadTree * q, real * coord, real weight, int id)
+{
+ if (!q)
+ return q;
+ return QuadTree_add_internal(q, coord, weight, id, 0);
+
}
+static void draw_polygon(FILE * fp, int dim, real * center, real width)
+{
+ /* pliot the enclosing square */
+ if (dim < 2 || dim > 3)
+ return;
+ fprintf(fp, "(*in c*){Line[{");
+
+ if (dim == 2) {
+ fprintf(fp, "{%f, %f}", center[0] + width, center[1] + width);
+ fprintf(fp, ",{%f, %f}", center[0] - width, center[1] + width);
+ fprintf(fp, ",{%f, %f}", center[0] - width, center[1] - width);
+ fprintf(fp, ",{%f, %f}", center[0] + width, center[1] - width);
+ fprintf(fp, ",{%f, %f}", center[0] + width, center[1] + width);
+ } else if (dim == 3) {
+ /* top */
+ fprintf(fp, "{");
+ fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width,
+ center[2] + width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width,
+ center[2] + width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width,
+ center[2] + width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width,
+ center[2] + width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width,
+ center[2] + width);
+ fprintf(fp, "},");
+ /* bot */
+ fprintf(fp, "{");
+ fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width,
+ center[2] - width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width,
+ center[2] - width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width,
+ center[2] - width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width,
+ center[2] - width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width,
+ center[2] - width);
+ fprintf(fp, "},");
+ /* for sides */
+ fprintf(fp, "{");
+ fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width,
+ center[2] - width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width,
+ center[2] + width);
+ fprintf(fp, "},");
+
+ fprintf(fp, "{");
+ fprintf(fp, "{%f, %f, %f}", center[0] - width, center[1] + width,
+ center[2] - width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width,
+ center[2] + width);
+ fprintf(fp, "},");
+
+ fprintf(fp, "{");
+ fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] - width,
+ center[2] - width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width,
+ center[2] + width);
+ fprintf(fp, "},");
+
+ fprintf(fp, "{");
+ fprintf(fp, "{%f, %f, %f}", center[0] - width, center[1] - width,
+ center[2] - width);
+ fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width,
+ center[2] + width);
+ fprintf(fp, "}");
+ }
-QuadTree QuadTree_add(QuadTree q, real *coord, real weight, int id){
- if (!q) return q;
- return QuadTree_add_internal(q, coord, weight, id, 0);
-
-}
-static void draw_polygon(FILE *fp, int dim, real *center, real width){
- /* pliot the enclosing square */
- if (dim < 2 || dim > 3) return;
- fprintf(fp,"(*in c*){Line[{");
-
- if (dim == 2){
- fprintf(fp, "{%f, %f}", center[0] + width, center[1] + width);
- fprintf(fp, ",{%f, %f}", center[0] - width, center[1] + width);
- fprintf(fp, ",{%f, %f}", center[0] - width, center[1] - width);
- fprintf(fp, ",{%f, %f}", center[0] + width, center[1] - width);
- fprintf(fp, ",{%f, %f}", center[0] + width, center[1] + width);
- } else if (dim == 3){
- /* top */
- fprintf(fp,"{");
- fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width, center[2] + width);
- fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width, center[2] + width);
- fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width, center[2] + width);
- fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width, center[2] + width);
- fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width, center[2] + width);
- fprintf(fp,"},");
- /* bot */
- fprintf(fp,"{");
- fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width, center[2] - width);
- fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width, center[2] - width);
- fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width, center[2] - width);
- fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width, center[2] - width);
- fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width, center[2] - width);
- fprintf(fp,"},");
- /* for sides */
- fprintf(fp,"{");
- fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] + width, center[2] - width);
- fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] + width, center[2] + width);
- fprintf(fp,"},");
-
- fprintf(fp,"{");
- fprintf(fp, "{%f, %f, %f}", center[0] - width, center[1] + width, center[2] - width);
- fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] + width, center[2] + width);
- fprintf(fp,"},");
-
- fprintf(fp,"{");
- fprintf(fp, "{%f, %f, %f}", center[0] + width, center[1] - width, center[2] - width);
- fprintf(fp, ",{%f, %f, %f}", center[0] + width, center[1] - width, center[2] + width);
- fprintf(fp,"},");
-
- fprintf(fp,"{");
- fprintf(fp, "{%f, %f, %f}", center[0] - width, center[1] - width, center[2] - width);
- fprintf(fp, ",{%f, %f, %f}", center[0] - width, center[1] - width, center[2] + width);
- fprintf(fp,"}");
- }
-
-
- fprintf(fp, "}]}(*end C*)");
+ fprintf(fp, "}]}(*end C*)");
}
-static void QuadTree_print_internal(FILE *fp, QuadTree q, int level){
- /* dump a quad tree in Mathematica format. */
- SingleLinkedList l, l0;
- real *coord;
- int i, dim = q->dim;
-
- draw_polygon(fp, q->dim, q->center, q->width);
-
- l0 = l = q->l;
- if (l){
- printf(",(*a*) {Red,");
- while (l){
- if (l != l0) printf(",");
- coord = node_data_get_coord(SingleLinkedList_get_data(l));
- fprintf(fp, "Point[{");
- for (i = 0; i < dim; i++){
- if (i != 0) printf(",");
- fprintf(fp, "%f",coord[i]);
- }
- fprintf(fp, "}]");
- l = SingleLinkedList_get_next(l);
+static void QuadTree_print_internal(FILE * fp, QuadTree * q, int level)
+{
+ /* dump a quad tree in Mathematica format. */
+ SingleLinkedList *l;
+ SingleLinkedList *l0;
+ real *coord;
+ int i, dim = q->dim;
+
+ draw_polygon(fp, q->dim, q->center, q->width);
+
+ l0 = l = q->l;
+ if (l) {
+ printf(",(*a*) {Red,");
+ while (l) {
+ if (l != l0)
+ printf(",");
+ coord = node_data_get_coord(SingleLinkedList_get_data(l));
+ fprintf(fp, "Point[{");
+ for (i = 0; i < dim; i++) {
+ if (i != 0)
+ printf(",");
+ fprintf(fp, "%f", coord[i]);
+ }
+ fprintf(fp, "}]");
+ l = SingleLinkedList_get_next(l);
+ }
+ fprintf(fp, "}");
}
- fprintf(fp, "}");
- }
-
- if (q->qts){
- for (i = 0; i < 1<<dim; i++){
- fprintf(fp, ",(*b*){");
- QuadTree_print_internal(fp, q->qts[i], level + 1);
- fprintf(fp, "}");
+
+ if (q->qts) {
+ for (i = 0; i < 1 << dim; i++) {
+ fprintf(fp, ",(*b*){");
+ QuadTree_print_internal(fp, q->qts[i], level + 1);
+ fprintf(fp, "}");
+ }
}
- }
}
-void QuadTree_print(FILE *fp, QuadTree q){
- if (!fp) return;
- if (q->dim == 2){
- fprintf(fp, "Graphics[{");
- } else if (q->dim == 3){
- fprintf(fp, "Graphics3D[{");
- } else {
- return;
- }
- QuadTree_print_internal(fp, q, 0);
- if (q->dim == 2){
- fprintf(fp, "}, PlotRange -> All, Frame -> True, FrameTicks -> True]\n");
- } else {
- fprintf(fp, "}, PlotRange -> All]\n");
- }
+void QuadTree_print(FILE * fp, QuadTree * q)
+{
+ if (!fp)
+ return;
+ if (q->dim == 2) {
+ fprintf(fp, "Graphics[{");
+ } else if (q->dim == 3) {
+ fprintf(fp, "Graphics3D[{");
+ } else {
+ return;
+ }
+ QuadTree_print_internal(fp, q, 0);
+ if (q->dim == 2) {
+ fprintf(fp,
+ "}, PlotRange -> All, Frame -> True, FrameTicks -> True]\n");
+ } else {
+ fprintf(fp, "}, PlotRange -> All]\n");
+ }
}
+/* vim:set shiftwidth=4 ts=8: */
+
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
+#ifndef QUAD_TREE_H
+#define QUAD_TREE_H
+
#include "LinkedList.h"
+#include "sfdpinternal.h"
+#include <stdio.h>
-typedef struct QuadTree_struct *QuadTree;
+typedef struct QuadTree_struct QuadTree;
struct QuadTree_struct {
- /* a data structure containing coordinates of n items, their average is in "average".
- The current level is a square or cube of width "width", which is subdivided into
- 2^dim QuadTrees qts. At the last level, all coordinates are stored in a single linked list l.
- total_weight is the combined weights of the nodes */
- int n;/* number of items */
- real total_weight;
- int dim;
- real *center;/* center of the bounding box, array of dimension dim. Allocated inside */
- real width;/* center +/- width gives the lower/upper bound, so really width is the
- "radius" */
- real *average;/* the average coordinates. Array of length dim. Allocated inside */
- QuadTree *qts;/* subtree . If dim = 2, there are 4, dim = 3 gives 8 */
- SingleLinkedList l;
- int max_level;
+ /* a data structure containing coordinates of n items, their average is in "average".
+ The current level is a square or cube of width "width", which is subdivided into
+ 2^dim QuadTrees qts. At the last level, all coordinates are stored in a single linked list l.
+ total_weight is the combined weights of the nodes */
+ int n; /* number of items */
+ real total_weight;
+ int dim;
+ real *center; /* center of the bounding box, array of dimension dim. Allocated inside */
+ real width; /* center +/- width gives the lower/upper bound, so really width is the
+ "radius" */
+ real *average; /* the average coordinates. Array of length dim. Allocated inside */
+ QuadTree **qts; /* subtree . If dim = 2, there are 4, dim = 3 gives 8 */
+ SingleLinkedList *l;
+ int max_level;
};
-QuadTree QuadTree_new(int dim, real *center, real width, int max_level);
+QuadTree *QuadTree_new(int dim, real * center, real width, int max_level);
-void QuadTree_delete(QuadTree q);
+void QuadTree_delete(QuadTree * q);
-QuadTree QuadTree_add(QuadTree q, real *coord, real weight, int id);/* coord is copied in */
+QuadTree *QuadTree_add(QuadTree * q, real * coord, real weight, int id); /* coord is copied in */
-void QuadTree_print(FILE *fp, QuadTree q);
+void QuadTree_print(FILE * fp, QuadTree * q);
-QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord, real *weight);
+QuadTree *QuadTree_new_from_point_list(int dim, int n, int max_level,
+ real * coord, real * weight);
-real point_distance(real *p1, real *p2, int dim);
+real point_distance(real * p1, real * p2, int dim);
-void QuadTree_get_supernodes(QuadTree qt, real bh, real *point, int nodeid, int *nsuper,
- int *nsupermax, real **center, real **supernode_wgts, real **distances, real *counts, int *flag);
+void QuadTree_get_supernodes(QuadTree * qt, real bh, real * point,
+ int nodeid, int *nsuper, int *nsupermax,
+ real ** center, real ** supernode_wgts,
+ real ** distances, real * counts, int *flag);
+#endif
-#include "general.h"
-#include "SparseMatrix.h"
-#include "spring_electrical.h"
+/* vim:set shiftwidth=4 ts=8: */
+
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
#include "post_process.h"
+#include "spring_electrical.h"
#include "sparse_solve.h"
#include "call_tri.h"
+#include "types.h"
+#include "memory.h"
+#include "globals.h"
#include <time.h>
+#include <string.h>
+#include <math.h>
#define node_degree(i) (ia[(i)+1] - ia[(i)])
-
-
-SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, real *x){
- /* find the ideal distance between edges, either 1, or |N[i] \Union N[j]| - |N[i] \Intersection N[j]|
- */
- SparseMatrix D;
- int *ia, *ja, i, j, k, l, nz;
- real *d;
- int *mask = NULL;
- real len, di, sum, sumd;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
-
- D = SparseMatrix_copy(A);
- ia = D->ia;
- ja = D->ja;
- if (D->type != MATRIX_TYPE_REAL){
- FREE(D->a);
- D->type = MATRIX_TYPE_REAL;
- D->a = MALLOC(sizeof(real)*D->nz);
- }
- d = (real*) D->a;
-
- mask = MALLOC(sizeof(int)*D->m);
- for (i = 0; i < D->m; i++) mask[i] = -1;
-
- for (i = 0; i < D->m; i++){
- di = node_degree(i);
- mask[i] = i;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- mask[ja[j]] = i;
+SparseMatrix *ideal_distance_matrix(SparseMatrix * A, int dim, real * x)
+{
+ /* find the ideal distance between edges, either 1, or |N[i] \Union N[j]| - |N[i] \Intersection N[j]|
+ */
+ SparseMatrix *D;
+ int *ia, *ja, i, j, k, l, nz;
+ real *d;
+ int *mask = NULL;
+ real len, di, sum, sumd;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+
+ D = SparseMatrix_copy(A);
+ ia = D->ia;
+ ja = D->ja;
+ if (D->type != MATRIX_TYPE_REAL) {
+ free(D->a);
+ D->type = MATRIX_TYPE_REAL;
+ D->a = N_GNEW(D->nz, real);
}
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- if (i == k) continue;
- len = di + node_degree(k);
- for (l = ia[k]; l < ia[k+1]; l++){
- if (mask[ja[l]] == i) len--;
- }
- d[j] = len;
- assert(len > 0);
+ d = (real *) D->a;
+
+ mask = N_GNEW(D->m, int);
+ for (i = 0; i < D->m; i++)
+ mask[i] = -1;
+
+ for (i = 0; i < D->m; i++) {
+ di = node_degree(i);
+ mask[i] = i;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ mask[ja[j]] = i;
+ }
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ if (i == k)
+ continue;
+ len = di + node_degree(k);
+ for (l = ia[k]; l < ia[k + 1]; l++) {
+ if (mask[ja[l]] == i)
+ len--;
+ }
+ d[j] = len;
+ assert(len > 0);
+ }
+
}
-
- }
-
- sum = 0; sumd = 0;
- nz = 0;
- for (i = 0; i < D->m; i++){
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- nz++;
- sum += distance(x, dim, i, ja[j]);
- sumd += d[j];
+
+ sum = 0;
+ sumd = 0;
+ nz = 0;
+ for (i = 0; i < D->m; i++) {
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ nz++;
+ sum += distance(x, dim, i, ja[j]);
+ sumd += d[j];
+ }
}
- }
- sum /= nz; sumd /= nz;
- sum = sum/sumd;
-
- for (i = 0; i < D->m; i++){
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- d[j] = sum*d[j];
+ sum /= nz;
+ sumd /= nz;
+ sum = sum / sumd;
+
+ for (i = 0; i < D->m; i++) {
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ d[j] = sum * d[j];
+ }
}
- }
- return D;
+ return D;
}
-StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda0, real *x,
- int ideal_dist_scheme){
- StressMajorizationSmoother sm;
- int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd;
- int *mask, nz;
- real *d, *w, *lambda;
- real *avg_dist, diag_d, diag_w, *dd, dist, s = 0, stop = 0, sbot = 0;
- SparseMatrix ID;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
-
- ID = ideal_distance_matrix(A, dim, x);
- dd = (real*) ID->a;
+StressMajorizationSmoother *StressMajorizationSmoother_new(SparseMatrix *
+ A, int dim,
+ real lambda0,
+ real * x,
+ int
+ ideal_dist_scheme)
+{
+ StressMajorizationSmoother *sm;
+ int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd;
+ int *mask, nz;
+ real *d, *w, *lambda;
+ real *avg_dist, diag_d, diag_w, *dd, dist, s = 0, stop = 0, sbot = 0;
+ SparseMatrix *ID;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+
+ ID = ideal_distance_matrix(A, dim, x);
+ dd = (real *) ID->a;
+
+ sm = GNEW(StressMajorizationSmoother);
+ lambda = sm->lambda = N_GNEW(m, real);
+ for (i = 0; i < m; i++)
+ sm->lambda[i] = lambda0;
+ mask = N_GNEW(m, int);
+
+ avg_dist = N_GNEW(m, real);
+
+ for (i = 0; i < m; i++) {
+ avg_dist[i] = 0;
+ nz = 0;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ avg_dist[i] += distance(x, dim, i, ja[j]);
+ nz++;
+ }
+ assert(nz > 0);
+ avg_dist[i] /= nz;
+ }
- sm = MALLOC(sizeof(struct StressMajorizationSmoother_struct));
- lambda = sm->lambda = MALLOC(sizeof(real)*m);
- for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
- mask = MALLOC(sizeof(int)*m);
- avg_dist = MALLOC(sizeof(real)*m);
+ for (i = 0; i < m; i++)
+ mask[i] = -1;
- for (i = 0; i < m ;i++){
- avg_dist[i] = 0;
nz = 0;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- avg_dist[i] += distance(x, dim, i, ja[j]);
- nz++;
+ for (i = 0; i < m; i++) {
+ mask[i] = i;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ if (mask[k] != i) {
+ mask[k] = i;
+ nz++;
+ }
+ }
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ for (l = ia[k]; l < ia[k + 1]; l++) {
+ if (mask[ja[l]] != i) {
+ mask[ja[l]] = i;
+ nz++;
+ }
+ }
+ }
}
- assert(nz > 0);
- avg_dist[i] /= nz;
- }
+ sm->Lw = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
+ sm->Lwd = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
+ if (!(sm->Lw) || !(sm->Lwd)) {
+ StressMajorizationSmoother_delete(sm);
+ return NULL;
+ }
- for (i = 0; i < m; i++) mask[i] = -1;
+ iw = sm->Lw->ia;
+ jw = sm->Lw->ja;
+ id = sm->Lwd->ia;
+ jd = sm->Lwd->ja;
+ w = (real *) sm->Lw->a;
+ d = (real *) sm->Lwd->a;
+ iw[0] = id[0] = 0;
- nz = 0;
- for (i = 0; i < m; i++){
- mask[i] = i;
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- if (mask[k] != i){
- mask[k] = i;
- nz++;
- }
- }
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- for (l = ia[k]; l < ia[k+1]; l++){
- if (mask[ja[l]] != i){
- mask[ja[l]] = i;
- nz++;
+ nz = 0;
+ for (i = 0; i < m; i++) {
+ mask[i] = i + m;
+ diag_d = diag_w = 0;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ if (mask[k] != i + m) {
+ mask[k] = i + m;
+
+ jw[nz] = k;
+ if (ideal_dist_scheme == IDEAL_GRAPH_DIST) {
+ dist = 1;
+ } else if (ideal_dist_scheme == IDEAL_AVG_DIST) {
+ dist = (avg_dist[i] + avg_dist[k]) * 0.5;
+ } else if (ideal_dist_scheme == IDEAL_POWER_DIST) {
+ dist = pow(distance_cropped(x, dim, i, k), .4);
+ } else {
+ fprintf(stderr, "ideal_dist_scheme value wrong");
+ assert(0);
+ exit(1);
+ }
+
+ /*
+ w[nz] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]);
+ w[nz] = -2./(avg_dist[i]+avg_dist[k]); */
+ /* w[nz] = -1.; *//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
+ w[nz] = -1 / (dist * dist);
+
+ diag_w += w[nz];
+
+ jd[nz] = k;
+ /*
+ d[nz] = w[nz]*distance(x,dim,i,k);
+ d[nz] = w[nz]*dd[j];
+ d[nz] = w[nz]*(avg_dist[i] + avg_dist[k])*0.5;
+ */
+ d[nz] = w[nz] * dist;
+ stop += d[nz] * distance(x, dim, i, k);
+ sbot += d[nz] * dist;
+ diag_d += d[nz];
+
+ nz++;
+ }
}
- }
- }
- }
- sm->Lw = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
- sm->Lwd = SparseMatrix_new(m, m, nz + m, MATRIX_TYPE_REAL, FORMAT_CSR);
- if (!(sm->Lw) || !(sm->Lwd)) {
- StressMajorizationSmoother_delete(sm);
- return NULL;
- }
-
- iw = sm->Lw->ia; jw = sm->Lw->ja;
- id = sm->Lwd->ia; jd = sm->Lwd->ja;
- w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
- iw[0] = id[0] = 0;
-
- nz = 0;
- for (i = 0; i < m; i++){
- mask[i] = i+m;
- diag_d = diag_w = 0;
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- if (mask[k] != i+m){
- mask[k] = i+m;
-
- jw[nz] = k;
- if (ideal_dist_scheme == IDEAL_GRAPH_DIST){
- dist = 1;
- } else if (ideal_dist_scheme == IDEAL_AVG_DIST){
- dist = (avg_dist[i] + avg_dist[k])*0.5;
- } else if (ideal_dist_scheme == IDEAL_POWER_DIST){
- dist = pow(distance_cropped(x,dim,i,k),.4);
- } else {
- fprintf(stderr,"ideal_dist_scheme value wrong");
- assert(0);
- exit(1);
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ for (l = ia[k]; l < ia[k + 1]; l++) {
+ if (mask[ja[l]] != i + m) {
+ mask[ja[l]] = i + m;
+
+ if (ideal_dist_scheme == IDEAL_GRAPH_DIST) {
+ dist = 2;
+ } else if (ideal_dist_scheme == IDEAL_AVG_DIST) {
+ dist =
+ (avg_dist[i] + 2 * avg_dist[k] +
+ avg_dist[ja[l]]) * 0.5;
+ } else if (ideal_dist_scheme == IDEAL_POWER_DIST) {
+ dist = pow(distance_cropped(x, dim, i, ja[l]), .4);
+ } else {
+ fprintf(stderr, "ideal_dist_scheme value wrong");
+ assert(0);
+ exit(1);
+ }
+
+ jw[nz] = ja[l];
+ /*
+ w[nz] = -1/(ia[i+1]-ia[i]+ia[ja[l]+1]-ia[ja[l]]);
+ w[nz] = -2/(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]]); */
+ /* w[nz] = -1.; *//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
+
+ w[nz] = -1 / (dist * dist);
+
+ diag_w += w[nz];
+
+ jd[nz] = ja[l];
+ /*
+ d[nz] = w[nz]*(distance(x,dim,i,k)+distance(x,dim,k,ja[l]));
+ d[nz] = w[nz]*(dd[j]+dd[l]);
+ d[nz] = w[nz]*(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
+ */
+ d[nz] = w[nz] * dist;
+ stop += d[nz] * distance(x, dim, ja[l], k);
+ sbot += d[nz] * dist;
+ diag_d += d[nz];
+
+ nz++;
+ }
+ }
}
+ jw[nz] = i;
+ lambda[i] *= (-diag_w); /* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
- /*
- w[nz] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]);
- w[nz] = -2./(avg_dist[i]+avg_dist[k]);*/
- /* w[nz] = -1.;*//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
- w[nz] = -1/(dist*dist);
-
- diag_w += w[nz];
-
- jd[nz] = k;
- /*
- d[nz] = w[nz]*distance(x,dim,i,k);
- d[nz] = w[nz]*dd[j];
- d[nz] = w[nz]*(avg_dist[i] + avg_dist[k])*0.5;
- */
- d[nz] = w[nz]*dist;
- stop += d[nz]*distance(x,dim,i,k);
- sbot += d[nz]*dist;
- diag_d += d[nz];
-
+ w[nz] = -diag_w + lambda[i];
+ jd[nz] = i;
+ d[nz] = -diag_d;
nz++;
- }
- }
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- for (l = ia[k]; l < ia[k+1]; l++){
- if (mask[ja[l]] != i+m){
- mask[ja[l]] = i+m;
-
- if (ideal_dist_scheme == IDEAL_GRAPH_DIST){
- dist = 2;
- } else if (ideal_dist_scheme == IDEAL_AVG_DIST){
- dist = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
- } else if (ideal_dist_scheme == IDEAL_POWER_DIST){
- dist = pow(distance_cropped(x,dim,i,ja[l]),.4);
- } else {
- fprintf(stderr,"ideal_dist_scheme value wrong");
- assert(0);
- exit(1);
- }
-
- jw[nz] = ja[l];
- /*
- w[nz] = -1/(ia[i+1]-ia[i]+ia[ja[l]+1]-ia[ja[l]]);
- w[nz] = -2/(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]]);*/
- /* w[nz] = -1.;*//* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
-
- w[nz] = -1/(dist*dist);
-
- diag_w += w[nz];
-
- jd[nz] = ja[l];
- /*
- d[nz] = w[nz]*(distance(x,dim,i,k)+distance(x,dim,k,ja[l]));
- d[nz] = w[nz]*(dd[j]+dd[l]);
- d[nz] = w[nz]*(avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
- */
- d[nz] = w[nz]*dist;
- stop += d[nz]*distance(x,dim,ja[l],k);
- sbot += d[nz]*dist;
- diag_d += d[nz];
-
- nz++;
- }
- }
+ iw[i + 1] = nz;
+ id[i + 1] = nz;
}
- jw[nz] = i;
- lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
-
- w[nz] = -diag_w + lambda[i];
- jd[nz] = i;
- d[nz] = -diag_d;
- nz++;
-
- iw[i+1] = nz;
- id[i+1] = nz;
- }
- s = stop/sbot;
- for (i = 0; i < nz; i++) d[i] *= s;
-
- sm->Lw->nz = nz;
- sm->Lwd->nz = nz;
-
- FREE(mask);
- FREE(avg_dist);
- SparseMatrix_delete(ID);
- return sm;
-}
-
-static real total_distance(int m, int dim, real* x, real* y){
- real total = 0, dist = 0;
- int i, j;
+ s = stop / sbot;
+ for (i = 0; i < nz; i++)
+ d[i] *= s;
- for (i = 0; i < m; i++){
- dist = 0.;
- for (j = 0; j < dim; j++){
- dist += (y[i*dim+j] - x[i*dim+j])*(y[i*dim+j] - x[i*dim+j]);
- }
- total += sqrt(dist);
- }
- return total;
+ sm->Lw->nz = nz;
+ sm->Lwd->nz = nz;
+ free(mask);
+ free(avg_dist);
+ SparseMatrix_delete(ID);
+ return sm;
}
+static real total_distance(int m, int dim, real * x, real * y)
+{
+ real total = 0, dist = 0;
+ int i, j;
+
+ for (i = 0; i < m; i++) {
+ dist = 0.;
+ for (j = 0; j < dim; j++) {
+ dist +=
+ (y[i * dim + j] - x[i * dim + j]) * (y[i * dim + j] -
+ x[i * dim + j]);
+ }
+ total += sqrt(dist);
+ }
+ return total;
+}
-void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit_sm) {
- SparseMatrix Lw = sm->Lw, Lwd = sm->Lwd, Lwdd = NULL;
- int i, j, m, *id, *jd, idiag, flag = 0, iter = 0;
- real *dd, *d, *y = NULL, *x0 = NULL, diag, diff = 1, tol = 0.001, *lambda = sm->lambda, maxit, res;
-
- Lwdd = SparseMatrix_copy(Lwd);
- m = Lw->m;
- x0 = MALLOC(sizeof(real)*dim*m);
- if (!x0) goto RETURN;
-
- x0 = MEMCPY(x0, x, sizeof(real)*dim*m);
- y = MALLOC(sizeof(real)*dim*m);
- if (!y) goto RETURN;
- id = Lwd->ia; jd = Lwd->ja;
- d = (real*) Lwd->a;
- dd = (real*) Lwdd->a;
- while (iter++ < maxit_sm && diff > tol){
- for (i = 0; i < m; i++){
- idiag = -1;
- diag = 0.;
- for (j = id[i]; j < id[i+1]; j++){
- if (i == jd[j]) {
- idiag = j;
- continue;
+void StressMajorizationSmoother_smooth(StressMajorizationSmoother * sm,
+ int dim, real * x, int maxit_sm)
+{
+ SparseMatrix *Lw = sm->Lw;
+ SparseMatrix *Lwd = sm->Lwd;
+ SparseMatrix *Lwdd = NULL;
+ int i, j, m, *id, *jd, idiag, flag = 0, iter = 0;
+ real *dd, *d, *y = NULL, *x0 = NULL, diag, diff = 1, tol =
+ 0.001, *lambda = sm->lambda, maxit, res;
+
+ Lwdd = SparseMatrix_copy(Lwd);
+ m = Lw->m;
+ x0 = N_GNEW(dim * m, real);
+ if (!x0)
+ goto RETURN;
+
+ x0 = memcpy(x0, x, sizeof(real) * dim * m);
+ y = N_GNEW(dim * m, real);
+ if (!y)
+ goto RETURN;
+
+ id = Lwd->ia;
+ jd = Lwd->ja;
+ d = (real *) Lwd->a;
+ dd = (real *) Lwdd->a;
+
+ while (iter++ < maxit_sm && diff > tol) {
+ for (i = 0; i < m; i++) {
+ idiag = -1;
+ diag = 0.;
+ for (j = id[i]; j < id[i + 1]; j++) {
+ if (i == jd[j]) {
+ idiag = j;
+ continue;
+ }
+ dd[j] = d[j] / distance_cropped(x, dim, i, jd[j]);
+ diag += dd[j];
+ }
+ assert(idiag >= 0);
+ dd[idiag] = -diag;
}
- dd[j] = d[j]/distance_cropped(x, dim, i, jd[j]);
- diag += dd[j];
- }
- assert(idiag >= 0);
- dd[idiag] = -diag;
- }
- /* solve (Lw+lambda*I) x = Lwdd y + lambda x0 */
+ /* solve (Lw+lambda*I) x = Lwdd y + lambda x0 */
- SparseMatrix_multiply_dense(Lwdd, FALSE, x, FALSE, &y, FALSE, dim);
- for (i = 0; i < m; i++){
- for (j = 0; j < dim; j++){
- y[i*dim+j] += lambda[i]*x0[i*dim+j];
- }
- }
+ SparseMatrix_multiply_dense(Lwdd, FALSE, x, FALSE, &y, FALSE, dim);
+ for (i = 0; i < m; i++) {
+ for (j = 0; j < dim; j++) {
+ y[i * dim + j] += lambda[i] * x0[i * dim + j];
+ }
+ }
- maxit = sqrt((double) m);
- res = SparseMatrix_solve(Lw, dim, x, y, 0.01, maxit, SOLVE_METHOD_CG, &flag);
+ maxit = sqrt((double) m);
+ res =
+ SparseMatrix_solve(Lw, dim, x, y, 0.01, maxit, SOLVE_METHOD_CG,
+ &flag);
- if (flag) goto RETURN;
+ if (flag)
+ goto RETURN;
- diff = total_distance(m, dim, x, y)/m;
+ diff = total_distance(m, dim, x, y) / m;
#ifdef DEBUG_PRINT
- if (Verbose){
- fprintf(stderr, "Outer iter = %d, res = %g Stress Majorization diff = %g\n",iter, res, diff);
- }
+ if (Verbose) {
+ fprintf(stderr,
+ "Outer iter = %d, res = %g Stress Majorization diff = %g\n",
+ iter, res, diff);
+ }
#endif
- MEMCPY(x, y, sizeof(real)*m*dim);
- }
+ memcpy(x, y, sizeof(real) * m * dim);
+ }
-#ifdef OVERLAP
- _statistics[1] += iter-1;
+#ifdef DEBUG
+ _statistics[1] += iter - 1;
#endif
- RETURN:
- SparseMatrix_delete(Lwdd);
- if (!x0) FREE(x0);
- if (!y) FREE(y);
+ RETURN:
+ SparseMatrix_delete(Lwdd);
+ if (!x0)
+ free(x0);
+ if (!y)
+ free(y);
}
-void StressMajorizationSmoother_delete(StressMajorizationSmoother sm){
- if (!sm) return;
- if (sm->Lw) SparseMatrix_delete(sm->Lw);
- if (sm->Lwd) SparseMatrix_delete(sm->Lwd);
- if (sm->lambda) FREE(sm->lambda);
+void StressMajorizationSmoother_delete(StressMajorizationSmoother * sm)
+{
+ if (!sm)
+ return;
+ if (sm->Lw)
+ SparseMatrix_delete(sm->Lw);
+ if (sm->Lwd)
+ SparseMatrix_delete(sm->Lwd);
+ if (sm->lambda)
+ free(sm->lambda);
}
-TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, int use_triangularization){
- TriangleSmoother sm;
- int i, j, k, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd, jdiag, nz;
- SparseMatrix B;
- real *avg_dist, *lambda, *d, *w, diag_d, diag_w, dist;
- real s = 0, stop = 0, sbot = 0;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
+TriangleSmoother *TriangleSmoother_new(SparseMatrix * A, int dim,
+ real lambda0, real * x,
+ int use_triangularization)
+{
+ TriangleSmoother *sm;
+ int i, j, k, m = A->m, *ia = A->ia, *ja =
+ A->ja, *iw, *jw, *id, *jd, jdiag, nz;
+ SparseMatrix *B;
+ real *avg_dist, *lambda, *d, *w, diag_d, diag_w, dist;
+ real s = 0, stop = 0, sbot = 0;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+
+ avg_dist = N_GNEW(m, real);
+
+ for (i = 0; i < m; i++) {
+ avg_dist[i] = 0;
+ nz = 0;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ avg_dist[i] += distance(x, dim, i, ja[j]);
+ nz++;
+ }
+ assert(nz > 0);
+ avg_dist[i] /= nz;
+ }
- avg_dist = MALLOC(sizeof(real)*m);
+ sm = GNEW(TriangleSmoother);
+ lambda = sm->lambda = N_GNEW(m, real);
+ for (i = 0; i < m; i++)
+ sm->lambda[i] = lambda0;
- for (i = 0; i < m ;i++){
- avg_dist[i] = 0;
- nz = 0;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- avg_dist[i] += distance(x, dim, i, ja[j]);
- nz++;
- }
- assert(nz > 0);
- avg_dist[i] /= nz;
- }
-
- sm = MALLOC(sizeof(struct TriangleSmoother_struct));
- lambda = sm->lambda = MALLOC(sizeof(real)*m);
- for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
-
- if (m > 2){
- if (use_triangularization){
- B= call_tri(m, dim, x);
+ if (m > 2) {
+ if (use_triangularization) {
+ B = call_tri(m, dim, x);
+ } else {
+ B = call_tri2(m, dim, x);
+ }
} else {
- B= call_tri2(m, dim, x);
+ B = SparseMatrix_copy(A);
}
- } else {
- B = SparseMatrix_copy(A);
- }
-
- /*
- {FILE *fp;
- fp = fopen("/tmp/111","w");
- export_embedding(fp, dim, B, x, NULL);
- fclose(fp);}
- */
-
-
- sm->Lw = SparseMatrix_add(A, B);
-
- SparseMatrix_delete(B);
- sm->Lwd = SparseMatrix_copy(sm->Lw);
- if (!(sm->Lw) || !(sm->Lwd)) {
- TriangleSmoother_delete(sm);
- return NULL;
- }
-
- iw = sm->Lw->ia; jw = sm->Lw->ja;
- id = sm->Lwd->ia; jd = sm->Lwd->ja;
- w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
-
- for (i = 0; i < m; i++){
- diag_d = diag_w = 0;
- jdiag = -1;
- for (j = iw[i]; j < iw[i+1]; j++){
- k = jw[j];
- if (k == i){
- jdiag = j;
- continue;
- }
- /* w[j] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]);
- w[j] = -2./(avg_dist[i]+avg_dist[k]);
- w[j] = -1.*/;/* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
- dist = pow(distance_cropped(x,dim,i,k),0.6);
- w[j] = 1/(dist*dist);
- diag_w += w[j];
-
- /* d[j] = w[j]*distance(x,dim,i,k);
- d[j] = w[j]*(avg_dist[i] + avg_dist[k])*0.5;*/
- d[j] = w[j]*dist;
- stop += d[j]*distance(x,dim,i,k);
- sbot += d[j]*dist;
- diag_d += d[j];
+ /*
+ {FILE *fp;
+ fp = fopen("/tmp/111","w");
+ export_embedding(fp, dim, B, x, NULL);
+ fclose(fp);}
+ */
+
+
+ sm->Lw = SparseMatrix_add(A, B);
+
+ SparseMatrix_delete(B);
+ sm->Lwd = SparseMatrix_copy(sm->Lw);
+ if (!(sm->Lw) || !(sm->Lwd)) {
+ TriangleSmoother_delete(sm);
+ return NULL;
}
- lambda[i] *= (-diag_w);/* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
+ iw = sm->Lw->ia;
+ jw = sm->Lw->ja;
+ id = sm->Lwd->ia;
+ jd = sm->Lwd->ja;
+ w = (real *) sm->Lw->a;
+ d = (real *) sm->Lwd->a;
+
+ for (i = 0; i < m; i++) {
+ diag_d = diag_w = 0;
+ jdiag = -1;
+ for (j = iw[i]; j < iw[i + 1]; j++) {
+ k = jw[j];
+ if (k == i) {
+ jdiag = j;
+ continue;
+ }
+ /* w[j] = -1./(ia[i+1]-ia[i]+ia[ja[j]+1]-ia[ja[j]]);
+ w[j] = -2./(avg_dist[i]+avg_dist[k]);
+ w[j] = -1. */ ;
+ /* use unit weight for now, later can try 1/(deg(i)+deg(k)) */
+ dist = pow(distance_cropped(x, dim, i, k), 0.6);
+ w[j] = 1 / (dist * dist);
+ diag_w += w[j];
+
+ /* d[j] = w[j]*distance(x,dim,i,k);
+ d[j] = w[j]*(avg_dist[i] + avg_dist[k])*0.5; */
+ d[j] = w[j] * dist;
+ stop += d[j] * distance(x, dim, i, k);
+ sbot += d[j] * dist;
+ diag_d += d[j];
+
+ }
+
+ lambda[i] *= (-diag_w); /* alternatively don't do that then we have a constant penalty term scaled by lambda0 */
- assert(jdiag >= 0);
- w[jdiag] = -diag_w + lambda[i];
- d[jdiag] = -diag_d;
- }
+ assert(jdiag >= 0);
+ w[jdiag] = -diag_w + lambda[i];
+ d[jdiag] = -diag_d;
+ }
- s = stop/sbot;
- for (i = 0; i < iw[m]; i++) d[i] *= s;
+ s = stop / sbot;
+ for (i = 0; i < iw[m]; i++)
+ d[i] *= s;
- FREE(avg_dist);
+ free(avg_dist);
- return sm;
+ return sm;
}
-void TriangleSmoother_delete(TriangleSmoother sm){
+void TriangleSmoother_delete(TriangleSmoother * sm)
+{
- StressMajorizationSmoother_delete(sm);
+ StressMajorizationSmoother_delete(sm);
}
-void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x){
+void TriangleSmoother_smooth(TriangleSmoother * sm, int dim, real * x)
+{
- StressMajorizationSmoother_smooth(sm, dim, x, 50);
+ StressMajorizationSmoother_smooth(sm, dim, x, 50);
}
/* ================================ spring and spring-electrical based smoother ================ */
-SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, real *x){
- SpringSmoother sm;
- int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *id, *jd;
- int *mask, nz;
- real *d, *dd;
- real *avg_dist;
- SparseMatrix ID = NULL;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
-
- ID = ideal_distance_matrix(A, dim, x);
- dd = (real*) ID->a;
+SpringSmoother *SpringSmoother_new(SparseMatrix * A, int dim,
+ spring_electrical_control * ctrl,
+ real * x)
+{
+ SpringSmoother *sm;
+ int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *id, *jd;
+ int *mask, nz;
+ real *d, *dd;
+ real *avg_dist;
+ SparseMatrix *ID = NULL;
+
+ assert(SparseMatrix_is_symmetric(A, FALSE));
+
+ ID = ideal_distance_matrix(A, dim, x);
+ dd = (real *) ID->a;
+
+ sm = GNEW(SpringSmoother);
+ mask = N_GNEW(m, int);
+
+ avg_dist = N_GNEW(m, real);
+
+ for (i = 0; i < m; i++) {
+ avg_dist[i] = 0;
+ nz = 0;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ if (i == ja[j])
+ continue;
+ avg_dist[i] += distance(x, dim, i, ja[j]);
+ nz++;
+ }
+ assert(nz > 0);
+ avg_dist[i] /= nz;
+ }
- sm = MALLOC(sizeof(struct SpringSmoother_struct));
- mask = MALLOC(sizeof(int)*m);
- avg_dist = MALLOC(sizeof(real)*m);
+ for (i = 0; i < m; i++)
+ mask[i] = -1;
- for (i = 0; i < m ;i++){
- avg_dist[i] = 0;
nz = 0;
- for (j = ia[i]; j < ia[i+1]; j++){
- if (i == ja[j]) continue;
- avg_dist[i] += distance(x, dim, i, ja[j]);
- nz++;
+ for (i = 0; i < m; i++) {
+ mask[i] = i;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ if (mask[k] != i) {
+ mask[k] = i;
+ nz++;
+ }
+ }
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ for (l = ia[k]; l < ia[k + 1]; l++) {
+ if (mask[ja[l]] != i) {
+ mask[ja[l]] = i;
+ nz++;
+ }
+ }
+ }
}
- assert(nz > 0);
- avg_dist[i] /= nz;
- }
+ sm->D = SparseMatrix_new(m, m, nz, MATRIX_TYPE_REAL, FORMAT_CSR);
+ if (!(sm->D)) {
+ SpringSmoother_delete(sm);
+ return NULL;
+ }
- for (i = 0; i < m; i++) mask[i] = -1;
+ id = sm->D->ia;
+ jd = sm->D->ja;
+ d = (real *) sm->D->a;
+ id[0] = 0;
- nz = 0;
- for (i = 0; i < m; i++){
- mask[i] = i;
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- if (mask[k] != i){
- mask[k] = i;
- nz++;
- }
- }
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- for (l = ia[k]; l < ia[k+1]; l++){
- if (mask[ja[l]] != i){
- mask[ja[l]] = i;
- nz++;
+ nz = 0;
+ for (i = 0; i < m; i++) {
+ mask[i] = i + m;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ if (mask[k] != i + m) {
+ mask[k] = i + m;
+ jd[nz] = k;
+ d[nz] = (avg_dist[i] + avg_dist[k]) * 0.5;
+ d[nz] = dd[j];
+ nz++;
+ }
}
- }
- }
- }
-
- sm->D = SparseMatrix_new(m, m, nz, MATRIX_TYPE_REAL, FORMAT_CSR);
- if (!(sm->D)){
- SpringSmoother_delete(sm);
- return NULL;
- }
-
- id = sm->D->ia; jd = sm->D->ja;
- d = (real*) sm->D->a;
- id[0] = 0;
-
- nz = 0;
- for (i = 0; i < m; i++){
- mask[i] = i+m;
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- if (mask[k] != i+m){
- mask[k] = i+m;
- jd[nz] = k;
- d[nz] = (avg_dist[i] + avg_dist[k])*0.5;
- d[nz] = dd[j];
- nz++;
- }
- }
- for (j = ia[i]; j < ia[i+1]; j++){
- k = ja[j];
- for (l = ia[k]; l < ia[k+1]; l++){
- if (mask[ja[l]] != i+m){
- mask[ja[l]] = i+m;
- jd[nz] = ja[l];
- d[nz] = (avg_dist[i] + 2*avg_dist[k] + avg_dist[ja[l]])*0.5;
- d[nz] = dd[j]+dd[l];
- nz++;
+ for (j = ia[i]; j < ia[i + 1]; j++) {
+ k = ja[j];
+ for (l = ia[k]; l < ia[k + 1]; l++) {
+ if (mask[ja[l]] != i + m) {
+ mask[ja[l]] = i + m;
+ jd[nz] = ja[l];
+ d[nz] =
+ (avg_dist[i] + 2 * avg_dist[k] +
+ avg_dist[ja[l]]) * 0.5;
+ d[nz] = dd[j] + dd[l];
+ nz++;
+ }
+ }
}
- }
+ id[i + 1] = nz;
}
- id[i+1] = nz;
- }
- sm->D->nz = nz;
- sm->ctrl = spring_electrical_control_new();
- *(sm->ctrl) = *ctrl;
- sm->ctrl->random_start = FALSE;
- sm->ctrl->multilevels = 1;
- sm->ctrl->step /= 2;
- sm->ctrl->maxiter = 20;
-
- FREE(mask);
- FREE(avg_dist);
- SparseMatrix_delete(ID);
-
- return sm;
+ sm->D->nz = nz;
+ sm->ctrl = spring_electrical_control_new();
+ *(sm->ctrl) = *ctrl;
+ sm->ctrl->random_start = FALSE;
+ sm->ctrl->multilevels = 1;
+ sm->ctrl->step /= 2;
+ sm->ctrl->maxiter = 20;
+
+ free(mask);
+ free(avg_dist);
+ SparseMatrix_delete(ID);
+
+ return sm;
}
-void SpringSmoother_delete(SpringSmoother sm){
- if (!sm) return;
- if (sm->D) SparseMatrix_delete(sm->D);
- if (sm->ctrl) spring_electrical_control_delete(sm->ctrl);
+void SpringSmoother_delete(SpringSmoother * sm)
+{
+ if (!sm)
+ return;
+ if (sm->D)
+ SparseMatrix_delete(sm->D);
+ if (sm->ctrl)
+ spring_electrical_control_delete(sm->ctrl);
}
-void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights, int dim, real *x){
- int flag = 0;
+void SpringSmoother_smooth(SpringSmoother * sm, SparseMatrix * A,
+ real * node_weights, int dim, real * x)
+{
+ int flag = 0;
- spring_electrical_spring_embedding(dim, A, sm->D, sm->ctrl, node_weights, x, &flag);
- assert(!flag);
+ spring_electrical_spring_embedding(dim, A, sm->D, sm->ctrl,
+ node_weights, x, &flag);
+ assert(!flag);
}
/*=============================== end of spring and spring-electrical based smoother =========== */
-void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
-#ifdef TIME
- clock_t cpu;
+int post_process_smoothing(int dim, SparseMatrix * A,
+ spring_electrical_control * ctrl,
+ real * node_weights, real * x)
+{
+#ifdef DEBUG
+ clock_t cpu = clock();
#endif
-#ifdef TIME
- cpu = clock();
-#endif
- *flag = 0;
-
- switch (ctrl->smoothing){
- case SMOOTHING_RNG:
- case SMOOTHING_TRIANGLE:{
- TriangleSmoother sm;
-
- if (ctrl->smoothing == SMOOTHING_RNG){
- sm = TriangleSmoother_new(A, dim, 0, x, FALSE);
- } else {
- sm = TriangleSmoother_new(A, dim, 0, x, TRUE);
- }
- TriangleSmoother_smooth(sm, dim, x);
- TriangleSmoother_delete(sm);
- break;
- }
- case SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST:
- case SMOOTHING_STRESS_MAJORIZATION_POWER_DIST:
- case SMOOTHING_STRESS_MAJORIZATION_AVG_DIST:
- {
- StressMajorizationSmoother sm;
- int k, dist_scheme = IDEAL_AVG_DIST;
-
- if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST){
- dist_scheme = IDEAL_GRAPH_DIST;
- } else if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_AVG_DIST){
- dist_scheme = IDEAL_AVG_DIST;
- } else if (ctrl->smoothing == SMOOTHING_STRESS_MAJORIZATION_POWER_DIST){
- dist_scheme = IDEAL_POWER_DIST;
- }
-
- for (k = 0; k < 1; k++){
- sm = StressMajorizationSmoother_new(A, dim, 0.05, x, dist_scheme);
- StressMajorizationSmoother_smooth(sm, dim, x, 50);
- StressMajorizationSmoother_delete(sm);
- }
- break;
- }
- case SMOOTHING_SPRING:{
- SpringSmoother sm;
- int k;
-
- for (k = 0; k < 1; k++){
- sm = SpringSmoother_new(A, dim, ctrl, x);
- SpringSmoother_smooth(sm, A, node_weights, dim, x);
- SpringSmoother_delete(sm);
- }
+ int flag = 0;
+
+ switch (ctrl->smoothing) {
+ case SMOOTHING_RNG:
+ case SMOOTHING_TRIANGLE:{
+ TriangleSmoother *sm;
+
+ if (ctrl->smoothing == SMOOTHING_RNG) {
+ sm = TriangleSmoother_new(A, dim, 0, x, FALSE);
+ } else {
+ sm = TriangleSmoother_new(A, dim, 0, x, TRUE);
+ }
+ TriangleSmoother_smooth(sm, dim, x);
+ TriangleSmoother_delete(sm);
+ break;
+ }
+ case SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST:
+ case SMOOTHING_STRESS_MAJORIZATION_POWER_DIST:
+ case SMOOTHING_STRESS_MAJORIZATION_AVG_DIST:
+ {
+ StressMajorizationSmoother *sm;
+ int k, dist_scheme = IDEAL_AVG_DIST;
+
+ if (ctrl->smoothing ==
+ SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST) {
+ dist_scheme = IDEAL_GRAPH_DIST;
+ } else if (ctrl->smoothing ==
+ SMOOTHING_STRESS_MAJORIZATION_AVG_DIST) {
+ dist_scheme = IDEAL_AVG_DIST;
+ } else if (ctrl->smoothing ==
+ SMOOTHING_STRESS_MAJORIZATION_POWER_DIST) {
+ dist_scheme = IDEAL_POWER_DIST;
+ }
+
+ for (k = 0; k < 1; k++) {
+ sm = StressMajorizationSmoother_new(A, dim, 0.05, x,
+ dist_scheme);
+ StressMajorizationSmoother_smooth(sm, dim, x, 50);
+ StressMajorizationSmoother_delete(sm);
+ }
+ break;
+ }
+ case SMOOTHING_SPRING:{
+ SpringSmoother *sm;
+ int k;
- break;
- }
+ for (k = 0; k < 1; k++) {
+ sm = SpringSmoother_new(A, dim, ctrl, x);
+ SpringSmoother_smooth(sm, A, node_weights, dim, x);
+ SpringSmoother_delete(sm);
+ }
+
+ break;
+ }
- }/* end switch between smoothing methods */
+ } /* end switch between smoothing methods */
-#ifdef TIME
- if (Verbose) fprintf(stderr, "post processing %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC);
+#ifdef DEBUG
+ if (Verbose)
+ fprintf(stderr, "post processing %f\n",
+ ((real) (clock() - cpu)) / CLOCKS_PER_SEC);
#endif
+ return flag;
}
-struct StressMajorizationSmoother_struct {
- SparseMatrix Lw;
- SparseMatrix Lwd;
- real* lambda;
-};
-
-typedef struct StressMajorizationSmoother_struct *StressMajorizationSmoother;
-
-void StressMajorizationSmoother_delete(StressMajorizationSmoother sm);
-
-enum {IDEAL_GRAPH_DIST, IDEAL_AVG_DIST, IDEAL_POWER_DIST};
-StressMajorizationSmoother StressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda, real *x, int ideal_dist_scheme);
-
-void StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit);
+/* vim:set shiftwidth=4 ts=8: */
+
+/**********************************************************
+* This software is part of the graphviz package *
+* http://www.graphviz.org/ *
+* *
+* Copyright (c) 1994-2004 AT&T Corp. *
+* and is licensed under the *
+* Common Public License, Version 1.0 *
+* by AT&T Corp. *
+* *
+* Information and Software Systems Research *
+* AT&T Research, Florham Park NJ *
+**********************************************************/
+
+#ifndef POST_PROCESS_H
+#define POST_PROCESS_H
+
+#include <spring_electrical.h>
+
+typedef struct {
+ SparseMatrix *Lw;
+ SparseMatrix *Lwd;
+ real *lambda;
+} StressMajorizationSmoother;
+
+void StressMajorizationSmoother_delete(StressMajorizationSmoother * sm);
+
+enum { IDEAL_GRAPH_DIST, IDEAL_AVG_DIST, IDEAL_POWER_DIST };
+StressMajorizationSmoother *StressMajorizationSmoother_new(SparseMatrix *
+ A, int dim,
+ real lambda,
+ real * x,
+ int
+ ideal_dist_scheme);
+
+void StressMajorizationSmoother_smooth(StressMajorizationSmoother * sm,
+ int dim, real * x, int maxit);
/*-------------------- triangle/neirhborhood graph based smoother ------------------- */
-typedef StressMajorizationSmoother TriangleSmoother;
+typedef StressMajorizationSmoother TriangleSmoother;
-#define TriangleSmoother_struct StressMajorizationSmoother_struct
+/* #define TriangleSmoother_struct StressMajorizationSmoother_struct */
-void TriangleSmoother_delete(TriangleSmoother sm);
+void TriangleSmoother_delete(TriangleSmoother * sm);
-TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda, real *x, int use_triangularization);
+TriangleSmoother *TriangleSmoother_new(SparseMatrix * A, int dim,
+ real lambda, real * x,
+ int use_triangularization);
-void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x);
+void TriangleSmoother_smooth(TriangleSmoother * sm, int dim, real * x);
/*------------------ spring and spring-electrical based smoother */
-struct SpringSmoother_struct {
- SparseMatrix D;
- spring_electrical_control ctrl;
-};
-
-typedef struct SpringSmoother_struct *SpringSmoother;
+typedef struct {
+ SparseMatrix *D;
+ spring_electrical_control *ctrl;
+} SpringSmoother;
-SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, real *x);
+SpringSmoother *SpringSmoother_new(SparseMatrix * A, int dim,
+ spring_electrical_control * ctrl,
+ real * x);
-void SpringSmoother_delete(SpringSmoother sm);
+void SpringSmoother_delete(SpringSmoother * sm);
-void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights, int dim, real *x);
+void SpringSmoother_smooth(SpringSmoother * sm, SparseMatrix * A,
+ real * node_weights, int dim, real * x);
/*------------------------------------------------------------------*/
-void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag);
+int post_process_smoothing(int dim, SparseMatrix * A,
+ spring_electrical_control * ctrl,
+ real * node_weights, real * x);
+
+#endif
projects / graphviz / commitdiff