From 5008fc2f749895e8e4cdce80b7cb07d3a8874e58 Mon Sep 17 00:00:00 2001 From: yifanhu Date: Wed, 14 May 2008 13:04:13 +0000 Subject: [PATCH] delete from head and move to branch att_07932 --- lib/sfdpgen/LinkedList.c | 63 - lib/sfdpgen/LinkedList.h | 36 - lib/sfdpgen/Makefile.am | 26 - lib/sfdpgen/Multilevel.c | 1193 ---------------- lib/sfdpgen/Multilevel.h | 77 - lib/sfdpgen/QuadTree.c | 519 ------- lib/sfdpgen/QuadTree.h | 61 - lib/sfdpgen/SparseMatrix.c | 2336 ------------------------------- lib/sfdpgen/SparseMatrix.h | 131 -- lib/sfdpgen/call_tri.c | 111 -- lib/sfdpgen/call_tri.h | 22 - lib/sfdpgen/layout_similarity.c | 197 --- lib/sfdpgen/overlap.c | 545 ------- lib/sfdpgen/overlap.h | 39 - lib/sfdpgen/post_process.c | 756 ---------- lib/sfdpgen/post_process.h | 76 - lib/sfdpgen/red_black_tree.c | 922 ------------ lib/sfdpgen/red_black_tree.h | 76 - lib/sfdpgen/sfdp.h | 29 - lib/sfdpgen/sfdpinit.c | 382 ----- lib/sfdpgen/sfdpinternal.h | 27 - lib/sfdpgen/sparse_solve.c | 230 --- lib/sfdpgen/sparse_solve.h | 36 - lib/sfdpgen/spring_electrical.c | 1340 ------------------ lib/sfdpgen/spring_electrical.h | 86 -- 25 files changed, 9316 deletions(-) delete mode 100644 lib/sfdpgen/LinkedList.c delete mode 100644 lib/sfdpgen/LinkedList.h delete mode 100644 lib/sfdpgen/Makefile.am delete mode 100644 lib/sfdpgen/Multilevel.c delete mode 100644 lib/sfdpgen/Multilevel.h delete mode 100644 lib/sfdpgen/QuadTree.c delete mode 100644 lib/sfdpgen/QuadTree.h delete mode 100644 lib/sfdpgen/SparseMatrix.c delete mode 100644 lib/sfdpgen/SparseMatrix.h delete mode 100644 lib/sfdpgen/call_tri.c delete mode 100644 lib/sfdpgen/call_tri.h delete mode 100644 lib/sfdpgen/layout_similarity.c delete mode 100644 lib/sfdpgen/overlap.c delete mode 100644 lib/sfdpgen/overlap.h delete mode 100644 lib/sfdpgen/post_process.c delete mode 100644 lib/sfdpgen/post_process.h delete mode 100644 lib/sfdpgen/red_black_tree.c delete mode 100644 lib/sfdpgen/red_black_tree.h delete mode 100644 lib/sfdpgen/sfdp.h delete mode 100644 lib/sfdpgen/sfdpinit.c delete mode 100644 lib/sfdpgen/sfdpinternal.h delete mode 100644 lib/sfdpgen/sparse_solve.c delete mode 100644 lib/sfdpgen/sparse_solve.h delete mode 100644 lib/sfdpgen/spring_electrical.c delete mode 100644 lib/sfdpgen/spring_electrical.h diff --git a/lib/sfdpgen/LinkedList.c b/lib/sfdpgen/LinkedList.c deleted file mode 100644 index ecd6c7c1c..000000000 --- a/lib/sfdpgen/LinkedList.c +++ /dev/null @@ -1,63 +0,0 @@ -/* 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 "LinkedList.h" -#include "memory.h" - -SingleLinkedList SingleLinkedList_new(void *data) -{ - SingleLinkedList head; - head = GNEW(struct SingleLinkedList_s); - head->data = data; - head->next = NULL; - return head; -} - -void SingleLinkedList_delete(SingleLinkedList head, - void (*linklist_deallocator) (void *)) -{ - SingleLinkedList next; - - if (!head) - return; - do { - next = head->next; - if (head->data) - linklist_deallocator(head->data); - if (head) - free(head); - head = next; - } while (head); - -} - - -SingleLinkedList SingleLinkedList_prepend(SingleLinkedList l, - void *data) -{ - SingleLinkedList head = SingleLinkedList_new(data); - head->next = l; - return head; -} - -void *SingleLinkedList_get_data(SingleLinkedList l) -{ - return l->data; -} - -SingleLinkedList SingleLinkedList_get_next(SingleLinkedList l) -{ - return l->next; -} diff --git a/lib/sfdpgen/LinkedList.h b/lib/sfdpgen/LinkedList.h deleted file mode 100644 index f0568c2c7..000000000 --- a/lib/sfdpgen/LinkedList.h +++ /dev/null @@ -1,36 +0,0 @@ -/* 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 LINKED_LIST_H -#define LINKED_LIST_H - -struct SingleLinkedList_s { - void *data; - struct SingleLinkedList_s *next; -}; - -typedef struct SingleLinkedList_s* SingleLinkedList; - -SingleLinkedList SingleLinkedList_new(void *data); -void SingleLinkedList_delete(SingleLinkedList head, - void (*linklist_deallocator) (void *)); -SingleLinkedList SingleLinkedList_prepend(SingleLinkedList l, - void *data); - -void *SingleLinkedList_get_data(SingleLinkedList l); - -SingleLinkedList SingleLinkedList_get_next(SingleLinkedList l); - -#endif diff --git a/lib/sfdpgen/Makefile.am b/lib/sfdpgen/Makefile.am deleted file mode 100644 index 57cf69d7b..000000000 --- a/lib/sfdpgen/Makefile.am +++ /dev/null @@ -1,26 +0,0 @@ -# $Id$ $Revision$ -## Process this file with automake to produce Makefile.in - -AM_CPPFLAGS = \ - -I$(top_srcdir) \ - -I$(top_srcdir)/lib/common \ - -I$(top_srcdir)/lib/gvc \ - -I$(top_srcdir)/lib/neatogen \ - -I$(top_srcdir)/lib/pack \ - -I$(top_srcdir)/lib/pathplan \ - -I$(top_srcdir)/lib/graph \ - -I$(top_srcdir)/lib/cdt - -noinst_HEADERS = SparseMatrix.h sfdpinternal.h spring_electrical.h \ - LinkedList.h overlap.h red_black_tree.h sparse_solve.h post_process.h \ - QuadTree.h call_tri.h Multilevel.h sfdp.h - -if WITH_SFDP -noinst_LTLIBRARIES = libsfdpgen_C.la -endif - -libsfdpgen_C_la_SOURCES = SparseMatrix.c sfdpinit.c spring_electrical.c \ - LinkedList.c overlap.c red_black_tree.c sparse_solve.c post_process.c \ - QuadTree.c call_tri.c Multilevel.c - -EXTRA_DIST = Makefile.old diff --git a/lib/sfdpgen/Multilevel.c b/lib/sfdpgen/Multilevel.c deleted file mode 100644 index ee6e2fb5b..000000000 --- a/lib/sfdpgen/Multilevel.c +++ /dev/null @@ -1,1193 +0,0 @@ -/* 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(struct Multilevel_control_s); - 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; -} - -void Multilevel_control_delete(Multilevel_control ctrl) -{ - free(ctrl); -} - -static Multilevel Multilevel_init(SparseMatrix A, real * node_weights) -{ - Multilevel grid; - if (!A) - return NULL; - assert(A->m == A->n); - grid = GNEW(struct Multilevel_s); - 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_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 int irand(int n){ - /* 0, 1, ..., n-1 */ - assert(n > 1); - /*return (int) MIN(floor(drand()*n),n-1);*/ - return rand()%n; -} - -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--; - } - 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 = 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); - } - (*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 = 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); - } -} - -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); - } -} - - - - - -#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 = 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); -#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; - } - } - - /* 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); -#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; - } - } - - /* dan yi dian, wu ban */ - for (i = 0; i < m; i++) { - if (matched[i] == i) { - (*cluster)[nz++] = i; - (*clusterp)[++(*ncluster)] = nz; - } - } - - free(p); - } - - 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 = 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; - } - - 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; - } - } - - /* dan yi dian, wu ban */ - for (i = 0; i < m; i++) { - if (matched[i] == i) { - (*cluster)[nz++] = i; - (*clusterp)[++(*ncluster)] = nz; - } - } - free(p); - - } - - free(super); - - free(superp); - - free(matched); -} - -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; -} - -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; - } - } - } - - 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; - } - } - free(p); - - - 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 = 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); - } -} - -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"); -#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); - - if (!(*cA)) { -#ifdef DEBUG_PRINT - 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); - } - - if (!(*cA)) { -#ifdef DEBUG_PRINT - 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); - } - - if (!(*cA)) { -#ifdef DEBUG_PRINT - 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); - } - - - if (!(*cA)) { -#ifdef DEBUG_PRINT - 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) { -#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); -#endif - 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); -#endif - 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); -#endif - 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); -} - -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; - -#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); - } -#endif - 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"); - } -#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 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_get_coarsest(Multilevel grid) -{ - while (grid->next) { - grid = grid->next; - } - return grid; -} diff --git a/lib/sfdpgen/Multilevel.h b/lib/sfdpgen/Multilevel.h deleted file mode 100644 index cc8f56b88..000000000 --- a/lib/sfdpgen/Multilevel.h +++ /dev/null @@ -1,77 +0,0 @@ -/* 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" - -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 */ -}; - -typedef struct Multilevel_s* Multilevel; - -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 { 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_s { - int minsize; - real min_coarsen_factor; - int maxlevel; - int randomize; - int coarsen_scheme; -}; - -typedef struct Multilevel_control_s* Multilevel_control; - -Multilevel_control Multilevel_control_new(); - -void Multilevel_control_delete(Multilevel_control ctrl); - -void Multilevel_delete(Multilevel grid); - -Multilevel Multilevel_new(SparseMatrix A, real * node_weights, - Multilevel_control ctrl); - -Multilevel Multilevel_get_coarsest(Multilevel grid); - -void print_padding(int n); - -#define Multilevel_is_finest(grid) (!((grid)->prev)) -#define Multilevel_is_coarsest(grid) (!((grid)->next)) - -#endif diff --git a/lib/sfdpgen/QuadTree.c b/lib/sfdpgen/QuadTree.c deleted file mode 100644 index 01810f2bc..000000000 --- a/lib/sfdpgen/QuadTree.c +++ /dev/null @@ -1,519 +0,0 @@ -/* 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 "QuadTree.h" -#include "memory.h" -#include "arith.h" -#include -#include - -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 = 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); -} - -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; -} - -int node_data_get_id(void *d) -{ - node_data *nd = (node_data *) d; - return nd->id; -} - -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) { - *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); - } - } - - 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 = 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 = 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); - } - } else { - for (i = 0; i < n; i++) { - qt = QuadTree_add(qt, &(coord[i * dim]), 1, i); - } - } - - - 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 = GNEW(struct QuadTree_s); - 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; - - 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); -} - -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; - } - } - return d; -} - -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 = 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 { - 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, "}"); - } - - - fprintf(fp, "}]}(*end C*)"); - - -} -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, "}"); - } - - 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"); - } -} diff --git a/lib/sfdpgen/QuadTree.h b/lib/sfdpgen/QuadTree.h deleted file mode 100644 index 07387892b..000000000 --- a/lib/sfdpgen/QuadTree.h +++ /dev/null @@ -1,61 +0,0 @@ -/* 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 - -typedef struct QuadTree_s* QuadTree; - -struct QuadTree_s { - /* 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); - -void QuadTree_delete(QuadTree q); - -QuadTree QuadTree_add(QuadTree q, real * coord, real weight, int id); /* coord is copied in */ - -void QuadTree_print(FILE * fp, QuadTree q); - -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); - -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 diff --git a/lib/sfdpgen/SparseMatrix.c b/lib/sfdpgen/SparseMatrix.c deleted file mode 100644 index bb016d5a3..000000000 --- a/lib/sfdpgen/SparseMatrix.c +++ /dev/null @@ -1,2336 +0,0 @@ -/* 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 "SparseMatrix.h" -#include "logic.h" -#include "math.h" -#include "string.h" -#include "memory.h" -#include "arith.h" -#include "assert.h" - -static size_t size_of_matrix_type(int type) -{ - int size = 0; - switch (type) { - case MATRIX_TYPE_REAL: - size = sizeof(real); - break; - case MATRIX_TYPE_COMPLEX: - size = 2 * sizeof(real); - break; - case MATRIX_TYPE_INTEGER: - size = sizeof(int); - break; - case MATRIX_TYPE_PATTERN: - size = 0; - break; - case MATRIX_TYPE_UNKNOWN: - size = 0; - break; - default: - size = 0; - break; - } - - return size; -} - -SparseMatrix SparseMatrix_make_undirected(SparseMatrix A) -{ - /* make it strictly low diag only, and set flag to undirected */ - SparseMatrix B; - B = SparseMatrix_symmetrize(A, FALSE); - SparseMatrix_set_undirected(B); - return SparseMatrix_remove_upper(B); -} - -SparseMatrix SparseMatrix_transpose(SparseMatrix A) -{ - int *ia = A->ia, *ja = A->ja, *ib, *jb, nz = A->nz, m = A->m, n = - A->n, type = A->type, format = A->format; - SparseMatrix B; - int i, j; - - if (!A) - return NULL; - assert(A->format == FORMAT_CSR); /* only implemented for CSR right now */ - - B = SparseMatrix_new(n, m, nz, type, format); - B->nz = nz; - ib = B->ia; - jb = B->ja; - - for (i = 0; i <= n; i++) - ib[i] = 0; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - ib[ja[j] + 1]++; - } - } - - for (i = 0; i < n; i++) - ib[i + 1] += ib[i]; - - switch (A->type) { - case MATRIX_TYPE_REAL:{ - real *a = (real *) A->a; - real *b = (real *) B->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jb[ib[ja[j]]] = i; - b[ib[ja[j]]++] = a[j]; - } - } - break; - } - case MATRIX_TYPE_COMPLEX:{ - real *a = (real *) A->a; - real *b = (real *) B->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jb[ib[ja[j]]] = i; - b[2 * ib[ja[j]]] = a[2 * j]; - b[2 * ib[ja[j]] + 1] = a[2 * j + 1]; - ib[ja[j]]++; - } - } - break; - } - case MATRIX_TYPE_INTEGER:{ - int *ai = (int *) A->a; - int *bi = (int *) B->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jb[ib[ja[j]]] = i; - bi[ib[ja[j]]++] = ai[j]; - } - } - break; - } - case MATRIX_TYPE_PATTERN: - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jb[ib[ja[j]]++] = i; - } - } - break; - case MATRIX_TYPE_UNKNOWN: - SparseMatrix_delete(B); - return NULL; - default: - SparseMatrix_delete(B); - return NULL; - } - - - for (i = n - 1; i >= 0; i--) - ib[i + 1] = ib[i]; - ib[0] = 0; - - - return B; -} - -SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, - int pattern_symmetric_only) -{ - SparseMatrix B; - if (SparseMatrix_is_symmetric(A, pattern_symmetric_only)) - return SparseMatrix_copy(A); - B = SparseMatrix_transpose(A); - if (!B) - return NULL; - A = SparseMatrix_add(A, B); - SparseMatrix_delete(B); - SparseMatrix_set_symmetric(A); - SparseMatrix_set_pattern_symmetric(A); - return A; -} - -int SparseMatrix_is_symmetric(SparseMatrix A, - int test_pattern_symmetry_only) -{ - /* assume no repeated entries! */ - SparseMatrix B; - int *ia, *ja, *ib, *jb, type, m; - int *mask; - int res = FALSE; - int i, j; - assert(A->format == FORMAT_CSR); /* only implemented for CSR right now */ - - if (!A) - return FALSE; - - if (SparseMatrix_known_symmetric(A)) - return TRUE; - if (test_pattern_symmetry_only - && SparseMatrix_known_strucural_symmetric(A)) - return TRUE; - - if (A->m != A->n) - return FALSE; - - B = SparseMatrix_transpose(A); - if (!B) - return FALSE; - - ia = A->ia; - ja = A->ja; - ib = B->ia; - jb = B->ja; - m = A->m; - - mask = N_NEW(m, int); - for (i = 0; i < m; i++) - mask[i] = -1; - - type = A->type; - if (test_pattern_symmetry_only) - type = MATRIX_TYPE_PATTERN; - - switch (type) { - case MATRIX_TYPE_REAL:{ - real *a = (real *) A->a; - real *b = (real *) B->a; - for (i = 0; i <= m; i++) - if (ia[i] != ib[i]) - goto RETURN; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - mask[ja[j]] = j; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (mask[jb[j]] < ia[i]) - goto RETURN; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (ABS(b[j] - a[mask[jb[j]]]) > SYMMETRY_EPSILON) - goto RETURN; - } - } - res = TRUE; - break; - } - case MATRIX_TYPE_COMPLEX:{ - real *a = (real *) A->a; - real *b = (real *) B->a; - for (i = 0; i <= m; i++) - if (ia[i] != ib[i]) - goto RETURN; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - mask[ja[j]] = j; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (mask[jb[j]] < ia[i]) - goto RETURN; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (ABS(b[2 * j] - a[2 * mask[jb[j]]]) > - SYMMETRY_EPSILON) - goto RETURN; - if (ABS(b[2 * j + 1] - a[2 * mask[jb[j]] + 1]) > - SYMMETRY_EPSILON) - goto RETURN; - } - } - res = TRUE; - break; - } - case MATRIX_TYPE_INTEGER:{ - int *ai = (int *) A->a; - int *bi = (int *) B->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - mask[ja[j]] = j; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (mask[jb[j]] < ia[i]) - goto RETURN; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (bi[j] != ai[mask[jb[j]]]) - goto RETURN; - } - } - break; - } - case MATRIX_TYPE_PATTERN: - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - mask[ja[j]] = j; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (mask[jb[j]] < ia[i]) - goto RETURN; - } - } - break; - case MATRIX_TYPE_UNKNOWN: - goto RETURN; - break; - default: - goto RETURN; - break; - } - - if (test_pattern_symmetry_only) { - SparseMatrix_set_pattern_symmetric(A); - } else { - SparseMatrix_set_symmetric(A); - SparseMatrix_set_pattern_symmetric(A); - } - RETURN: - free(mask); - - SparseMatrix_delete(B); - return res; -} - -static SparseMatrix SparseMatrix_init(int m, int n, int type, int format) -{ - SparseMatrix A; - - A = GNEW(struct SparseMatrix_s); - A->m = m; - A->n = n; - A->nz = 0; - A->nzmax = 0; - A->type = type; - switch (format) { - case FORMAT_COORD: - A->ia = NULL; - break; - case FORMAT_CSC: - case FORMAT_CSR: - default: - A->ia = N_GNEW(m + 1, int); - } - A->ja = NULL; - A->a = NULL; - A->format = format; - A->type = type; - A->property = 0; - clear_flag(A->property, MATRIX_PATTERN_SYMMETRIC); - clear_flag(A->property, MATRIX_SYMMETRIC); - clear_flag(A->property, MATRIX_SKEW); - clear_flag(A->property, MATRIX_HERMITIAN); - return A; -} - -static SparseMatrix SparseMatrix_alloc(SparseMatrix A, int nz) -{ - int type = A->type, format = A->format; - - A->a = NULL; - switch (format) { - case FORMAT_COORD: - A->ia = N_NEW(nz, int); - A->ja = N_NEW(nz, int); - A->a = gmalloc(size_of_matrix_type(type) * nz); - break; - case FORMAT_CSR: - case FORMAT_CSC: - default: - A->ja = N_NEW(nz, int); - if (size_of_matrix_type(type) && nz > 0) - A->a = gmalloc(size_of_matrix_type(type) * nz); - break; - } - A->nzmax = nz; - return A; -} - -static SparseMatrix SparseMatrix_realloc(SparseMatrix A, int nz) -{ - int type = A->type, format = A->format; - switch (format) { - case FORMAT_COORD: - A->ia = RALLOC(nz, A->ia, int); - A->ja = RALLOC(nz, A->ja, int); - if (size_of_matrix_type(type) > 0) { - if (A->a) { - A->a = grealloc(A->a, size_of_matrix_type(type) * nz); - } else { - A->a = gmalloc(size_of_matrix_type(type) * nz); - } - } - break; - case FORMAT_CSR: - case FORMAT_CSC: - default: - A->ja = RALLOC(nz, A->ja, int); - if (size_of_matrix_type(type) > 0) { - if (A->a) { - A->a = grealloc(A->a, size_of_matrix_type(type) * nz); - } else { - A->a = gmalloc(size_of_matrix_type(type) * nz); - } - } - break; - } - A->nzmax = nz; - return A; -} - -SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format) -{ - /* return a sparse matrix skeleton with row dimension m and storage nz. If nz == 0, - only row pointers are allocated */ - SparseMatrix A; - - A = SparseMatrix_init(m, n, type, format); - if (nz > 0) - A = SparseMatrix_alloc(A, nz); - return A; - -} - -void SparseMatrix_delete(SparseMatrix A) -{ - /* return a sparse matrix skeleton with row dimension m and storage nz. If nz == 0, - only row pointers are allocated */ - if (!A) - return; - if (A->ia) - free(A->ia); - if (A->ja) - free(A->ja); - if (A->a) - free(A->a); - free(A); -} -void SparseMatrix_print_csr(char *c, SparseMatrix A) -{ - int *ia, *ja; - real *a; - int *ai; - int i, j, m = A->m; - - printf("%s\n SparseArray[{", c); - ia = A->ia; - ja = A->ja; - a = A->a; - switch (A->type) { - case MATRIX_TYPE_REAL: - a = (real *) A->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - printf("{%d, %d}->%f", i + 1, ja[j] + 1, a[j]); - if (j != ia[m] - 1) - printf(","); - } - } - break; - case MATRIX_TYPE_COMPLEX: - a = (real *) A->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - printf("{%d, %d}->%f + %f I", i + 1, ja[j] + 1, a[2 * j], - a[2 * j + 1]); - if (j != ia[m] - 1) - printf(","); - } - } - printf("\n"); - break; - case MATRIX_TYPE_INTEGER: - ai = (int *) A->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - printf("{%d, %d}->%d", i + 1, ja[j] + 1, ai[j]); - if (j != ia[m] - 1) - printf(","); - } - } - printf("\n"); - break; - case MATRIX_TYPE_PATTERN: - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - printf("{%d, %d}->_", i + 1, ja[j] + 1); - if (j != ia[m] - 1) - printf(","); - } - } - printf("\n"); - break; - case MATRIX_TYPE_UNKNOWN: - return; - default: - return; - } - printf("},{%d, %d}]\n", m, A->n); - -} - - - -void SparseMatrix_print_coord(char *c, SparseMatrix A) -{ - int *ia, *ja; - real *a; - int *ai; - int i, m = A->m; - - printf("%s\n SparseArray[{", c); - ia = A->ia; - ja = A->ja; - a = A->a; - switch (A->type) { - case MATRIX_TYPE_REAL: - a = (real *) A->a; - for (i = 0; i < A->nz; i++) { - printf("{%d, %d}->%f", ia[i] + 1, ja[i] + 1, a[i]); - if (i != A->nz - 1) - printf(","); - } - printf("\n"); - break; - case MATRIX_TYPE_COMPLEX: - a = (real *) A->a; - for (i = 0; i < A->nz; i++) { - printf("{%d, %d}->%f + %f I", ia[i] + 1, ja[i] + 1, a[2 * i], - a[2 * i + 1]); - if (i != A->nz - 1) - printf(","); - } - printf("\n"); - break; - case MATRIX_TYPE_INTEGER: - ai = (int *) A->a; - for (i = 0; i < A->nz; i++) { - printf("{%d, %d}->%d", ia[i] + 1, ja[i] + 1, ai[i]); - if (i != A->nz) - printf(","); - } - printf("\n"); - break; - case MATRIX_TYPE_PATTERN: - for (i = 0; i < A->nz; i++) { - printf("{%d, %d}->_", ia[i] + 1, ja[i] + 1); - if (i != A->nz - 1) - printf(","); - } - printf("\n"); - break; - case MATRIX_TYPE_UNKNOWN: - return; - default: - return; - } - printf("},{%d, %d}]\n", m, A->n); - -} - - - - -void SparseMatrix_print(char *c, SparseMatrix A) -{ - switch (A->format) { - case FORMAT_CSR: - SparseMatrix_print_csr(c, A); - break; - case FORMAT_CSC: - assert(0); /* not implemented yet... - SparseMatrix_print_csc(c, A); */ - break; - case FORMAT_COORD: - SparseMatrix_print_coord(c, A); - break; - default: - assert(0); - } -} - - - - - -static void SparseMatrix_export_csr(FILE * f, SparseMatrix A) -{ - int *ia, *ja; - real *a; - int *ai; - int i, j, m = A->m; - - fprintf(f, "%d %d %d\n", A->m, A->n, A->nz); - ia = A->ia; - ja = A->ja; - a = A->a; - switch (A->type) { - case MATRIX_TYPE_REAL: - a = (real *) A->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - fprintf(f, "%d %d %16.8g\n", i + 1, ja[j] + 1, a[j]); - } - } - break; - case MATRIX_TYPE_COMPLEX: - a = (real *) A->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - fprintf(f, "%d %d %16.8g %16.8g\n", i + 1, ja[j] + 1, - a[2 * j], a[2 * j + 1]); - } - } - break; - case MATRIX_TYPE_INTEGER: - ai = (int *) A->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - fprintf(f, "%d %d %d\n", i + 1, ja[j] + 1, ai[j]); - } - } - break; - case MATRIX_TYPE_PATTERN: - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - fprintf(f, "%d %d\n", i + 1, ja[j] + 1); - } - } - break; - case MATRIX_TYPE_UNKNOWN: - return; - default: - return; - } - -} - -void SparseMatrix_export_binary(char *name, SparseMatrix A, int *flag) -{ - FILE *f; - - *flag = 0; - f = fopen(name, "wb"); - if (!f) { - *flag = 1; - return; - } - fwrite(&(A->m), sizeof(int), 1, f); - fwrite(&(A->n), sizeof(int), 1, f); - fwrite(&(A->nz), sizeof(int), 1, f); - fwrite(&(A->nzmax), sizeof(int), 1, f); - fwrite(&(A->type), sizeof(int), 1, f); - fwrite(&(A->format), sizeof(int), 1, f); - fwrite(&(A->property), sizeof(int), 1, f); - if (A->format == FORMAT_COORD) { - fwrite(A->ia, sizeof(int), A->nz, f); - } else { - fwrite(A->ia, sizeof(int), A->m + 1, f); - } - fwrite(A->ja, sizeof(int), A->nz, f); - if (size_of_matrix_type(A->type) > 0) - fwrite(A->a, size_of_matrix_type(A->type), A->nz, f); - fclose(f); - -} - -SparseMatrix SparseMatrix_import_binary(char *name) -{ - SparseMatrix A = NULL; - int m, n, nz, nzmax, type, format, property; - FILE *f; - - f = fopen(name, "rb"); - - if (!f) - return NULL; - fread(&m, sizeof(int), 1, f); - fread(&n, sizeof(int), 1, f); - fread(&nz, sizeof(int), 1, f); - fread(&nzmax, sizeof(int), 1, f); - fread(&type, sizeof(int), 1, f); - fread(&format, sizeof(int), 1, f); - fread(&property, sizeof(int), 1, f); - A = SparseMatrix_new(m, n, nz, type, format); - A->nz = nz; - A->property = property; - - if (format == FORMAT_COORD) { - fread(A->ia, sizeof(int), A->nz, f); - } else { - fread(A->ia, sizeof(int), A->m + 1, f); - } - fread(A->ja, sizeof(int), A->nz, f); - if (size_of_matrix_type(A->type) > 0) { - fread(A->a, size_of_matrix_type(A->type), A->nz, f); - } - fclose(f); - return A; -} - -static void SparseMatrix_export_coord(FILE * f, SparseMatrix A) -{ - int *ia, *ja; - real *a; - int *ai; - int i; - - fprintf(f, "%d %d %d\n", A->m, A->n, A->nz); - ia = A->ia; - ja = A->ja; - a = A->a; - switch (A->type) { - case MATRIX_TYPE_REAL: - a = (real *) A->a; - for (i = 0; i < A->nz; i++) { - fprintf(f, "%d %d %16.8g\n", ia[i] + 1, ja[i] + 1, a[i]); - } - break; - case MATRIX_TYPE_COMPLEX: - a = (real *) A->a; - for (i = 0; i < A->nz; i++) { - fprintf(f, "%d %d %16.8g %16.8g\n", ia[i] + 1, ja[i] + 1, - a[2 * i], a[2 * i + 1]); - } - break; - case MATRIX_TYPE_INTEGER: - ai = (int *) A->a; - for (i = 0; i < A->nz; i++) { - fprintf(f, "%d %d %d\n", ia[i] + 1, ja[i] + 1, ai[i]); - } - break; - case MATRIX_TYPE_PATTERN: - for (i = 0; i < A->nz; i++) { - fprintf(f, "%d %d\n", ia[i] + 1, ja[i] + 1); - } - break; - case MATRIX_TYPE_UNKNOWN: - return; - default: - return; - } -} - - - -void SparseMatrix_export(FILE * f, SparseMatrix A) -{ - - switch (A->format) { - case FORMAT_CSR: - SparseMatrix_export_csr(f, A); - break; - case FORMAT_CSC: - assert(0); /* not implemented yet... - SparseMatrix_print_csc(c, A); */ - break; - case FORMAT_COORD: - SparseMatrix_export_coord(f, A); - break; - default: - assert(0); - } -} - - -SparseMatrix SparseMatrix_from_coordinate_format(SparseMatrix A) -{ - /* convert a sparse matrix in coordinate form to one in compressed row form. */ - int *irn, *jcn; - - void *a = A->a; - - assert(A->format == FORMAT_COORD); - if (A->format != FORMAT_COORD) { - return NULL; - } - irn = A->ia; - jcn = A->ja; - return SparseMatrix_from_coordinate_arrays(A->nz, A->m, A->n, irn, jcn, - a, A->type); - -} - -SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, - int *irn, int *jcn, - void *val0, int type) -{ - /* convert a sparse matrix in coordinate form to one in compressed row form. - nz: number of entries - irn: row indices 0-based - jcn: column indices 0-based - val values if not NULL - type: matrix type - */ - - SparseMatrix A = NULL; - int *ia, *ja; - real *a, *val; - int *ai, *vali; - int i; - - assert(m > 0 && n > 0 && nz >= 0); - - if (m <= 0 || n <= 0 || nz < 0) - return NULL; - A = SparseMatrix_new(m, n, nz, type, FORMAT_CSR); - assert(A); - if (!A) - return NULL; - ia = A->ia; - ja = A->ja; - - for (i = 0; i <= m; i++) { - ia[i] = 0; - } - - switch (type) { - case MATRIX_TYPE_REAL: - val = (real *) val0; - a = (real *) A->a; - for (i = 0; i < nz; i++) { - if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) { - assert(0); - return NULL; - } - ia[irn[i] + 1]++; - } - for (i = 0; i < m; i++) - ia[i + 1] += ia[i]; - for (i = 0; i < nz; i++) { - a[ia[irn[i]]] = val[i]; - ja[ia[irn[i]]++] = jcn[i]; - } - for (i = m; i > 0; i--) - ia[i] = ia[i - 1]; - ia[0] = 0; - break; - case MATRIX_TYPE_COMPLEX: - val = (real *) val0; - a = (real *) A->a; - for (i = 0; i < nz; i++) { - if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) { - assert(0); - return NULL; - } - ia[irn[i] + 1]++; - } - for (i = 0; i < m; i++) - ia[i + 1] += ia[i]; - for (i = 0; i < nz; i++) { - a[2 * ia[irn[i]]] = *(val++); - a[2 * ia[irn[i]] + 1] = *(val++); - ja[ia[irn[i]]++] = jcn[i]; - } - for (i = m; i > 0; i--) - ia[i] = ia[i - 1]; - ia[0] = 0; - break; - case MATRIX_TYPE_INTEGER: - vali = (int *) val0; - ai = (int *) A->a; - for (i = 0; i < nz; i++) { - if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) { - assert(0); - return NULL; - } - ia[irn[i] + 1]++; - } - for (i = 0; i < m; i++) - ia[i + 1] += ia[i]; - for (i = 0; i < nz; i++) { - ai[ia[irn[i]]] = vali[i]; - ja[ia[irn[i]]++] = jcn[i]; - } - for (i = m; i > 0; i--) - ia[i] = ia[i - 1]; - ia[0] = 0; - break; - case MATRIX_TYPE_PATTERN: - for (i = 0; i < nz; i++) { - if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) { - assert(0); - return NULL; - } - ia[irn[i] + 1]++; - } - for (i = 0; i < m; i++) - ia[i + 1] += ia[i]; - for (i = 0; i < nz; i++) { - ja[ia[irn[i]]++] = jcn[i]; - } - for (i = m; i > 0; i--) - ia[i] = ia[i - 1]; - ia[0] = 0; - break; - case MATRIX_TYPE_UNKNOWN: - assert(0); - return NULL; - default: - assert(0); - return NULL; - } - A->nz = nz; - A = SparseMatrix_sum_repeat_entries(A); - return A; -} - - -SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B) -{ - int m, n; - SparseMatrix C = NULL; - int *mask = NULL; - int *ia = A->ia, *ja = A->ja, *ib = B->ia, *jb = B->ja, *ic, *jc; - int i, j, nz, nzmax; - - assert(A && B); - assert(A->format == B->format && A->format == FORMAT_CSR); /* other format not yet supported */ - assert(A->type == B->type); - m = A->m; - n = A->n; - if (m != B->m || n != B->n) - return NULL; - - nzmax = A->nz + B->nz; /* just assume that no entries overlaps for speed */ - - C = SparseMatrix_new(m, n, nzmax, A->type, FORMAT_CSR); - if (!C) - goto RETURN; - ic = C->ia; - jc = C->ja; - - mask = N_GNEW(n, int); - - for (i = 0; i < n; i++) - mask[i] = -1; - - nz = 0; - ic[0] = 0; - switch (A->type) { - case MATRIX_TYPE_REAL:{ - real *a = (real *) A->a; - real *b = (real *) B->a; - real *c = (real *) C->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - mask[ja[j]] = nz; - jc[nz] = ja[j]; - c[nz] = a[j]; - nz++; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (mask[jb[j]] < ic[i]) { - jc[nz] = jb[j]; - c[nz++] = b[j]; - } else { - c[mask[jb[j]]] += b[j]; - } - } - ic[i + 1] = nz; - } - break; - } - case MATRIX_TYPE_COMPLEX:{ - real *a = (real *) A->a; - real *b = (real *) B->a; - real *c = (real *) C->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - mask[ja[j]] = nz; - jc[nz] = ja[j]; - c[2 * nz] = a[2 * j]; - c[2 * nz + 1] = a[2 * j + 1]; - nz++; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (mask[jb[j]] < ic[i]) { - jc[nz] = jb[j]; - c[2 * nz] = b[2 * j]; - c[2 * nz + 1] = b[2 * j + 1]; - nz++; - } else { - c[2 * mask[jb[j]]] += b[2 * j]; - c[2 * mask[jb[j]] + 1] += b[2 * j + 1]; - } - } - ic[i + 1] = nz; - } - break; - } - case MATRIX_TYPE_INTEGER:{ - int *a = (int *) A->a; - int *b = (int *) B->a; - int *c = (int *) C->a; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - mask[ja[j]] = nz; - jc[nz] = ja[j]; - c[nz] = a[j]; - nz++; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (mask[jb[j]] < ic[i]) { - jc[nz] = jb[j]; - c[nz] = b[j]; - nz++; - } else { - c[mask[jb[j]]] += b[j]; - } - } - ic[i + 1] = nz; - } - break; - } - case MATRIX_TYPE_PATTERN:{ - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - mask[ja[j]] = nz; - jc[nz] = ja[j]; - nz++; - } - for (j = ib[i]; j < ib[i + 1]; j++) { - if (mask[jb[j]] < ic[i]) { - jc[nz] = jb[j]; - nz++; - } - } - ic[i + 1] = nz; - } - break; - } - case MATRIX_TYPE_UNKNOWN: - break; - default: - break; - } - C->nz = nz; - - RETURN: - if (mask) - free(mask); - - return C; -} - - - -static void dense_transpose(real * v, int m, int n) -{ - /* transpose an m X n matrix in place. Well, we do no really do it without xtra memory. This is possibe, but too complicated for ow */ - int i, j; - real *u; - u = N_GNEW(m * n, real); - memcpy(u, v, sizeof(real) * m * n); - - for (i = 0; i < m; i++) { - for (j = 0; j < n; j++) { - v[j * m + i] = u[i * n + j]; - } - } - free(u); -} - - -static void SparseMatrix_multiply_dense1(SparseMatrix A, real * v, - real ** res, int dim, - int transposed, - int res_transposed) -{ - /* A v or A^T v where v a dense matrix of second dimension dim. Real only for now. */ - int i, j, k, *ia, *ja, n, m; - real *a, *u; - - assert(A->format == FORMAT_CSR); - assert(A->type == MATRIX_TYPE_REAL); - - a = (real *) A->a; - ia = A->ia; - ja = A->ja; - m = A->m; - n = A->n; - u = *res; - - if (!transposed) { - if (!u) - u = N_GNEW(m * dim, real); - for (i = 0; i < m; i++) { - for (k = 0; k < dim; k++) - u[i * dim + k] = 0.; - for (j = ia[i]; j < ia[i + 1]; j++) { - for (k = 0; k < dim; k++) - u[i * dim + k] += a[j] * v[ja[j] * dim + k]; - } - } - if (res_transposed) - dense_transpose(u, m, dim); - } else { - if (!u) - u = N_GNEW(n * dim, real); - for (i = 0; i < n * dim; i++) - u[i] = 0.; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - for (k = 0; k < dim; k++) - u[ja[j] * dim + k] += a[j] * v[i * dim + k]; - } - } - if (res_transposed) - dense_transpose(u, n, dim); - } - - *res = u; - - -} - -static void SparseMatrix_multiply_dense2(SparseMatrix A, real * v, - real ** res, int dim, - int transposed, - int res_transposed) -{ - /* A v^T or A^T v^T where v a dense matrix of second dimension n or m. Real only for now. - transposed = FALSE: A*V^T, with A dimension m x n, V dimension dim x n, v[i*n+j] gives V[i,j]. Result of dimension m x dim - transposed = TRUE: A^T*V^T, with A dimension m x n, V dimension dim x m. v[i*m+j] gives V[i,j]. Result of dimension n x dim - */ - real *u, *rr; - int i, m, n; - assert(A->format == FORMAT_CSR); - assert(A->type == MATRIX_TYPE_REAL); - u = *res; - m = A->m; - n = A->n; - - if (!transposed) { - if (!u) - u = N_GNEW(m * dim, real); - for (i = 0; i < dim; i++) { - rr = &(u[m * i]); - SparseMatrix_multiply_vector(A, &(v[n * i]), &rr, transposed); - } - if (!res_transposed) - dense_transpose(u, dim, m); - } else { - if (!u) - u = N_GNEW(n * dim, real); - for (i = 0; i < dim; i++) { - rr = &(u[n * i]); - SparseMatrix_multiply_vector(A, &(v[m * i]), &rr, transposed); - } - if (!res_transposed) - dense_transpose(u, dim, n); - } - - *res = u; - - -} - - - -void SparseMatrix_multiply_dense(SparseMatrix A, int ATransposed, - real * v, int vTransposed, real ** res, - int res_transposed, int dim) -{ - /* depend on value of {ATranspose, vTransposed}, assume res_transposed == FALSE - {FALSE, FALSE}: A * V, with A dimension m x n, with V of dimension n x dim. v[i*dim+j] gives V[i,j]. Result of dimension m x dim - {TRUE, FALSE}: A^T * V, with A dimension m x n, with V of dimension m x dim. v[i*dim+j] gives V[i,j]. Result of dimension n x dim - {FALSE, TRUE}: A*V^T, with A dimension m x n, V dimension dim x n, v[i*n+j] gives V[i,j]. Result of dimension m x dim - {TRUE, TRUE}: A^T*V^T, with A dimension m x n, V dimension dim x m. v[i*m+j] gives V[i,j]. Result of dimension n x dim - - furthermore, if res_transpose d== TRUE, then the result is transposed. Hence if res_transposed == TRUE - - {FALSE, FALSE}: V^T A^T, with A dimension m x n, with V of dimension n x dim. v[i*dim+j] gives V[i,j]. Result of dimension dim x dim - {TRUE, FALSE}: V^T A, with A dimension m x n, with V of dimension m x dim. v[i*dim+j] gives V[i,j]. Result of dimension dim x n - {FALSE, TRUE}: V*A^T, with A dimension m x n, V dimension dim x n, v[i*n+j] gives V[i,j]. Result of dimension dim x m - {TRUE, TRUE}: V A, with A dimension m x n, V dimension dim x m. v[i*m+j] gives V[i,j]. Result of dimension dim x n - */ - - if (!vTransposed) { - SparseMatrix_multiply_dense1(A, v, res, dim, ATransposed, - res_transposed); - } else { - SparseMatrix_multiply_dense2(A, v, res, dim, ATransposed, - res_transposed); - } - -} - - - -void SparseMatrix_multiply_vector(SparseMatrix A, real * v, real ** res, - int transposed) -{ - /* A v or A^T v. Real only for now. */ - int i, j, *ia, *ja, n, m; - real *a, *u = NULL; - assert(A->format == FORMAT_CSR); - assert(A->type == MATRIX_TYPE_REAL); - - a = (real *) A->a; - ia = A->ia; - ja = A->ja; - m = A->m; - n = A->n; - u = *res; - - if (v) { - if (!transposed) { - if (!u) - u = N_GNEW(m, real); - for (i = 0; i < m; i++) { - u[i] = 0.; - for (j = ia[i]; j < ia[i + 1]; j++) { - u[i] += a[j] * v[ja[j]]; - } - } - } else { - if (!u) - u = N_GNEW(n, real); - for (i = 0; i < n; i++) - u[i] = 0.; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - u[ja[j]] += a[j] * v[i]; - } - } - } - } else { - /* v is assumed to be all 1's */ - if (!transposed) { - if (!u) - u = N_GNEW(m, real); - for (i = 0; i < m; i++) { - u[i] = 0.; - for (j = ia[i]; j < ia[i + 1]; j++) { - u[i] += a[j]; - } - } - } else { - if (!u) - u = N_GNEW(n, real); - for (i = 0; i < n; i++) - u[i] = 0.; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - u[ja[j]] += a[j]; - } - } - } - } - *res = u; - -} - - - -SparseMatrix SparseMatrix_scaled_by_vector(SparseMatrix A, real * v, - int apply_to_row) -{ - /* A SCALED BY VECOTR V IN ROW/COLUMN. Real only for now. */ - int i, j, *ia, *ja, m; - real *a; - assert(A->format == FORMAT_CSR); - assert(A->type == MATRIX_TYPE_REAL); - - a = (real *) A->a; - ia = A->ia; - ja = A->ja; - m = A->m; - - - if (!apply_to_row) { - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - a[j] *= v[ja[j]]; - } - } - } else { - for (i = 0; i < m; i++) { - if (v[i] != 0) { - for (j = ia[i]; j < ia[i + 1]; j++) { - a[j] *= v[i]; - } - } - } - } - return A; - -} - -SparseMatrix SparseMatrix_multiply_by_scaler(SparseMatrix A, real s) -{ - /* A scaled by a number */ - int i, j, *ia, m; - real *a; - assert(A->format == FORMAT_CSR); - assert(A->type == MATRIX_TYPE_REAL); - - a = (real *) A->a; - ia = A->ia; - - m = A->m; - - - - - - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - a[j] *= s; - } - } - - return A; - -} - - -SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B) -{ - int m, n; - SparseMatrix C = NULL; - int *mask = NULL; - int *ia = A->ia, *ja = A->ja, *ib = B->ia, *jb = B->ja, *ic, *jc; - int i, j, k, jj, type, nz; - - assert(A->format == B->format && A->format == FORMAT_CSR); /* other format not yet supported */ - - m = A->m; - n = A->n; - if (n != B->m) - return NULL; - if (A->type != B->type) { -#ifdef DEBUG - printf - ("in SparseMatrix_multiply, the matrix types do not match, right now only multiplication of matrices of the same type is supported\n"); -#endif - return NULL; - } - type = A->type; - - mask = N_GNEW(n, int); - if (!mask) - return NULL; - - for (i = 0; i < n; i++) - mask[i] = -1; - - nz = 0; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jj = ja[j]; - for (k = ib[jj]; k < ib[jj + 1]; k++) { - if (mask[jb[k]] != -i - 2) { - nz++; - mask[jb[k]] = -i - 2; - } - } - } - } - - C = SparseMatrix_new(m, B->n, nz, type, FORMAT_CSR); - if (!C) - goto RETURN; - ic = C->ia; - jc = C->ja; - - nz = 0; - - switch (type) { - case MATRIX_TYPE_REAL: - { - real *a = (real *) A->a; - real *b = (real *) B->a; - real *c = (real *) C->a; - ic[0] = 0; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jj = ja[j]; - for (k = ib[jj]; k < ib[jj + 1]; k++) { - if (mask[jb[k]] < ic[i]) { - mask[jb[k]] = nz; - jc[nz] = jb[k]; - c[nz] = a[j] * b[k]; - nz++; - } else { - assert(jc[mask[jb[k]]] == jb[k]); - c[mask[jb[k]]] += a[j] * b[k]; - } - } - } - ic[i + 1] = nz; - } - } - break; - case MATRIX_TYPE_COMPLEX: - { - real *a = (real *) A->a; - real *b = (real *) B->a; - real *c = (real *) C->a; - a = (real *) A->a; - b = (real *) B->a; - c = (real *) C->a; - ic[0] = 0; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jj = ja[j]; - for (k = ib[jj]; k < ib[jj + 1]; k++) { - if (mask[jb[k]] < ic[i]) { - mask[jb[k]] = nz; - jc[nz] = jb[k]; - c[2 * nz] = a[2 * j] * b[2 * k] - a[2 * j + 1] * b[2 * k + 1]; /*real part */ - c[2 * nz + 1] = a[2 * j] * b[2 * k + 1] + a[2 * j + 1] * b[2 * k]; /*img part */ - nz++; - } else { - assert(jc[mask[jb[k]]] == jb[k]); - c[2 * mask[jb[k]]] += a[2 * j] * b[2 * k] - a[2 * j + 1] * b[2 * k + 1]; /*real part */ - c[2 * mask[jb[k]] + 1] += a[2 * j] * b[2 * k + 1] + a[2 * j + 1] * b[2 * k]; /*img part */ - } - } - } - ic[i + 1] = nz; - } - } - break; - case MATRIX_TYPE_INTEGER: - { - int *a = (int *) A->a; - int *b = (int *) B->a; - int *c = (int *) C->a; - ic[0] = 0; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jj = ja[j]; - for (k = ib[jj]; k < ib[jj + 1]; k++) { - if (mask[jb[k]] < ic[i]) { - mask[jb[k]] = nz; - jc[nz] = jb[k]; - c[nz] += a[j] * b[k]; - nz++; - } else { - assert(jc[mask[jb[k]]] == jb[k]); - c[mask[jb[k]]] += a[j] * b[k]; - } - } - } - ic[i + 1] = nz; - } - } - break; - case MATRIX_TYPE_PATTERN: - ic[0] = 0; - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jj = ja[j]; - for (k = ib[jj]; k < ib[jj + 1]; k++) { - if (mask[jb[k]] < ic[i]) { - mask[jb[k]] = nz; - jc[nz] = jb[k]; - nz++; - } else { - assert(jc[mask[jb[k]]] == jb[k]); - } - } - } - ic[i + 1] = nz; - } - break; - case MATRIX_TYPE_UNKNOWN: - default: - SparseMatrix_delete(C); - C = NULL; - goto RETURN; - break; - } - - C->nz = nz; - - RETURN: - free(mask); - return C; - -} - -SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A) -{ - /* sum repeated entries in the same row, i.e., {1,1}->1, {1,1}->2 becomes {1,1}->3 */ - int *ia = A->ia, *ja = A->ja, type = A->type, n = A->n; - int *mask = NULL, nz, i, j, sta; - - mask = N_GNEW(n, int); - for (i = 0; i < n; i++) - mask[i] = -1; - - switch (type) { - case MATRIX_TYPE_REAL: - { - real *a = (real *) A->a; - nz = 0; - sta = ia[0]; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (mask[ja[j]] < ia[i]) { - ja[nz] = ja[j]; - a[nz] = a[j]; - mask[ja[j]] = nz++; - } else { - assert(ja[mask[ja[j]]] == ja[j]); - a[mask[ja[j]]] += a[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - } - break; - case MATRIX_TYPE_COMPLEX: - { - real *a = (real *) A->a; - nz = 0; - sta = ia[0]; - for (i = 0; i < A->m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - if (mask[ja[j]] < ia[i]) { - ja[nz] = ja[j]; - a[2 * nz] = a[2 * j]; - a[2 * nz + 1] = a[2 * j + 1]; - mask[ja[j]] = nz++; - } else { - assert(ja[mask[ja[j]]] == ja[j]); - a[2 * mask[ja[j]]] += a[2 * j]; - a[2 * mask[ja[j]] + 1] += a[2 * j + 1]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - } - break; - case MATRIX_TYPE_INTEGER: - { - int *a = (int *) A->a; - nz = 0; - sta = ia[0]; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (mask[ja[j]] < ia[i]) { - ja[nz] = ja[j]; - a[nz] = a[j]; - mask[ja[j]] = nz++; - } else { - assert(ja[mask[ja[j]]] == ja[j]); - a[mask[ja[j]]] += a[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - } - break; - case MATRIX_TYPE_PATTERN: - { - nz = 0; - sta = ia[0]; - for (i = 0; i < A->m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - if (mask[ja[j]] < ia[i]) { - ja[nz] = ja[j]; - mask[ja[j]] = nz++; - } else { - assert(ja[mask[ja[j]]] == ja[j]); - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - } - break; - case MATRIX_TYPE_UNKNOWN: - return NULL; - break; - default: - return NULL; - break; - } - A->nz = nz; - free(mask); - return A; -} - -SparseMatrix SparseMatrix_coordinate_form_add_entries(SparseMatrix A, - int nentries, - int *irn, int *jcn, - void *val) -{ - int nz, type = A->type, nzmax; - - assert(A->format == FORMAT_COORD); - if (nentries <= 0) - return A; - nz = A->nz; - nzmax = A->nzmax; - - if (nz + nentries >= A->nzmax) { - nzmax = nz + nentries; - A->nzmax = nzmax = MAX(10, (int) 0.2 * nzmax) + nzmax; - A = SparseMatrix_realloc(A, nzmax); - } - memcpy((char *) A->ia + nz * sizeof(int) / sizeof(char), irn, - sizeof(int) * nentries); - memcpy((char *) A->ja + nz * sizeof(int) / sizeof(char), jcn, - sizeof(int) * nentries); - if (size_of_matrix_type(type)) - memcpy((char *) A->a + - nz * size_of_matrix_type(type) / sizeof(char), val, - size_of_matrix_type(type) * nentries); - A->nz += nentries; - return A; -} - - -SparseMatrix SparseMatrix_remove_diagonal(SparseMatrix A) -{ - int i, j, *ia, *ja, nz, sta; - - if (!A) - return A; - - nz = 0; - ia = A->ia; - ja = A->ja; - sta = ia[0]; - switch (A->type) { - case MATRIX_TYPE_REAL:{ - real *a = (real *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ja[j] != i) { - ja[nz] = ja[j]; - a[nz++] = a[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_COMPLEX:{ - real *a = (real *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ja[j] != i) { - ja[nz] = ja[j]; - a[2 * nz] = a[2 * j]; - a[2 * nz + 1] = a[2 * j + 1]; - nz++; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_INTEGER:{ - int *a = (int *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ja[j] != i) { - ja[nz] = ja[j]; - a[nz++] = a[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_PATTERN:{ - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ja[j] != i) { - ja[nz++] = ja[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_UNKNOWN: - return NULL; - default: - return NULL; - } - - return A; -} - - -SparseMatrix SparseMatrix_remove_upper(SparseMatrix A) -{ /* remove diag and upper diag */ - int i, j, *ia, *ja, nz, sta; - - if (!A) - return A; - - nz = 0; - ia = A->ia; - ja = A->ja; - sta = ia[0]; - switch (A->type) { - case MATRIX_TYPE_REAL:{ - real *a = (real *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ja[j] < i) { - ja[nz] = ja[j]; - a[nz++] = a[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_COMPLEX:{ - real *a = (real *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ja[j] < i) { - ja[nz] = ja[j]; - a[2 * nz] = a[2 * j]; - a[2 * nz + 1] = a[2 * j + 1]; - nz++; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_INTEGER:{ - int *a = (int *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ja[j] < i) { - ja[nz] = ja[j]; - a[nz++] = a[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_PATTERN:{ - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ja[j] < i) { - ja[nz++] = ja[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_UNKNOWN: - return NULL; - default: - return NULL; - } - - clear_flag(A->property, MATRIX_PATTERN_SYMMETRIC); - clear_flag(A->property, MATRIX_SYMMETRIC); - clear_flag(A->property, MATRIX_SKEW); - clear_flag(A->property, MATRIX_HERMITIAN); - return A; -} - - - - -SparseMatrix SparseMatrix_divide_row_by_degree(SparseMatrix A) -{ - int i, j, *ia, *ja; - real deg; - - if (!A) - return A; - - ia = A->ia; - ja = A->ja; - switch (A->type) { - case MATRIX_TYPE_REAL:{ - real *a = (real *) A->a; - for (i = 0; i < A->m; i++) { - deg = ia[i + 1] - ia[i]; - for (j = ia[i]; j < ia[i + 1]; j++) { - a[j] = a[j] / deg; - } - } - break; - } - case MATRIX_TYPE_COMPLEX:{ - real *a = (real *) A->a; - for (i = 0; i < A->m; i++) { - deg = ia[i + 1] - ia[i]; - for (j = ia[i]; j < ia[i + 1]; j++) { - if (ja[j] != i) { - a[2 * j] = a[2 * j] / deg; - a[2 * j + 1] = a[2 * j + 1] / deg; - } - } - } - break; - } - case MATRIX_TYPE_INTEGER:{ - assert(0); /* this operation would not make sense for int matrix */ - break; - } - case MATRIX_TYPE_PATTERN:{ - break; - } - case MATRIX_TYPE_UNKNOWN: - return NULL; - default: - return NULL; - } - - return A; -} - - -SparseMatrix -SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A) -{ - /* symmetric, all entries to 1, diaginal removed */ - int i, *ia, *ja, nz, m, n; - real *a; - SparseMatrix B; - - if (!A) - return A; - - nz = A->nz; - ia = A->ia; - ja = A->ja; - n = A->n; - m = A->m; - - if (n != m) - return NULL; - - B = SparseMatrix_new(m, n, nz, MATRIX_TYPE_PATTERN, FORMAT_CSR); - - memcpy(B->ia, ia, sizeof(int) * (m + 1)); - memcpy(B->ja, ja, sizeof(int) * nz); - B->nz = A->nz; - - A = SparseMatrix_symmetrize(B, TRUE); - SparseMatrix_delete(B); - A = SparseMatrix_remove_diagonal(A); - A->a = N_GNEW(A->nz, real); - a = (real *) A->a; - for (i = 0; i < A->nz; i++) - a[i] = 1.; - A->type = MATRIX_TYPE_REAL; - return A; -} - - - -SparseMatrix SparseMatrix_normalize_to_rowsum1(SparseMatrix A) -{ - int i, j; - real sum, *a; - - if (!A) - return A; - if (A->format != FORMAT_CSR && A->type != MATRIX_TYPE_REAL) { -#ifdef DEBUG - printf("only CSR and real matrix supported.\n"); -#endif - return A; - } - - a = (real *) A->a; - for (i = 0; i < A->m; i++) { - sum = 0; - for (j = A->ia[i]; j < A->ia[i + 1]; j++) { - sum += a[j]; - } - if (sum != 0) { - for (j = A->ia[i]; j < A->ia[i + 1]; j++) { - a[j] /= sum; - } - } - } - return A; -} - - - -SparseMatrix SparseMatrix_normalize_by_row(SparseMatrix A) -{ - int i, j; - real max, *a; - - if (!A) - return A; - if (A->format != FORMAT_CSR && A->type != MATRIX_TYPE_REAL) { -#ifdef DEBUG - printf("only CSR and real matrix supported.\n"); -#endif - return A; - } - - a = (real *) A->a; - for (i = 0; i < A->m; i++) { - max = 0; - for (j = A->ia[i]; j < A->ia[i + 1]; j++) { - max = MAX(ABS(a[j]), max); - } - if (max != 0) { - for (j = A->ia[i]; j < A->ia[i + 1]; j++) { - a[j] /= max; - } - } - } - return A; -} - - -SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, - double (*fun) (double x)) -{ - int i, j; - real *a; - - if (!A) - return A; - if (A->format != FORMAT_CSR && A->type != MATRIX_TYPE_REAL) { -#ifdef DEBUG - printf("only CSR and real matrix supported.\n"); -#endif - return A; - } - - a = (real *) A->a; - for (i = 0; i < A->m; i++) { - for (j = A->ia[i]; j < A->ia[i + 1]; j++) { - a[j] = fun(a[j]); - } - } - return A; -} - - -SparseMatrix SparseMatrix_crop(SparseMatrix A, real epsilon) -{ - int i, j, *ia, *ja, nz, sta; - - if (!A) - return A; - - nz = 0; - ia = A->ia; - ja = A->ja; - sta = ia[0]; - switch (A->type) { - case MATRIX_TYPE_REAL:{ - real *a = (real *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ABS(a[j]) > epsilon) { - ja[nz] = ja[j]; - a[nz++] = a[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_COMPLEX:{ - real *a = (real *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (sqrt - (a[2 * j] * a[2 * j] + - a[2 * j + 1] * a[2 * j + 1]) > epsilon) { - ja[nz] = ja[j]; - a[2 * nz] = a[2 * j]; - a[2 * nz + 1] = a[2 * j + 1]; - nz++; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_INTEGER:{ - int *a = (int *) A->a; - for (i = 0; i < A->m; i++) { - for (j = sta; j < ia[i + 1]; j++) { - if (ABS(a[j]) > epsilon) { - ja[nz] = ja[j]; - a[nz++] = a[j]; - } - } - sta = ia[i + 1]; - ia[i + 1] = nz; - } - A->nz = nz; - break; - } - case MATRIX_TYPE_PATTERN:{ - break; - } - case MATRIX_TYPE_UNKNOWN: - return NULL; - default: - return NULL; - } - - return A; -} - -SparseMatrix SparseMatrix_copy(SparseMatrix A) -{ - SparseMatrix B; - if (!A) - return A; - B = SparseMatrix_new(A->m, A->n, A->nz, A->type, A->format); - memcpy(B->ia, A->ia, sizeof(int) * (A->m + 1)); - memcpy(B->ja, A->ja, sizeof(int) * (A->ia[A->m])); - if (A->a) - memcpy(B->a, A->a, size_of_matrix_type(A->type) * A->nz); - B->property = A->property; - B->nz = A->nz; - return B; -} - -int SparseMatrix_has_diagonal(SparseMatrix A) -{ - - int i, j, m = A->m, *ia = A->ia, *ja = A->ja; - - for (i = 0; i < m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - if (i == ja[j]) - return TRUE; - } - } - return FALSE; -} - -void SparseMatrix_level_sets(SparseMatrix A, int root, int *nlevel, - int **levelset_ptr, int **levelset, - int **mask, int reinitialize_mask) -{ - /* mask is assumed to be initialized to negative if provided. */ - int i, j, sta = 0, sto = 1, nz, ii; - int m = A->m, *ia = A->ia, *ja = A->ja; - - if (!(*levelset_ptr)) - *levelset_ptr = N_GNEW(m + 1, int); - if (!(*levelset)) - *levelset = N_GNEW(m, int); - if (!(*mask)) { - *mask = malloc(sizeof(int) * m); - for (i = 0; i < m; i++) - (*mask)[i] = UNMASKED; - } -#ifdef DEBUG - for (i = 0; i < m; i++) - assert((*mask)[i] < 0); -#endif - *nlevel = 0; - - assert(root >= 0 && root < m); - (*levelset_ptr)[0] = 0; - (*levelset_ptr)[1] = 1; - (*levelset)[0] = root; - (*mask)[root] = 1; - *nlevel = 1; - nz = 1; - sta = 0; - sto = 1; - while (sto > sta) { - for (i = sta; i < sto; i++) { - ii = (*levelset)[i]; - for (j = ia[ii]; j < ia[ii + 1]; j++) { - if (ii == ja[j]) - continue; - if ((*mask)[ja[j]] < 0) { - (*levelset)[nz++] = ja[j]; - (*mask)[ja[j]] = 1; - } - } - } - (*levelset_ptr)[++(*nlevel)] = nz; - sta = sto; - sto = nz; - } - (*nlevel)--; - if (reinitialize_mask) - for (i = 0; i < (*levelset_ptr)[*nlevel]; i++) - (*mask)[(*levelset)[i]] = UNMASKED; -} - -void SparseMatrix_weakly_connected_components(SparseMatrix A0, - int *ncomp, int **comps, - int **comps_ptr) -{ - SparseMatrix A = A0; - int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL, nlevel; - int m = A->m, i, nn; - - if (!SparseMatrix_is_symmetric(A, TRUE)) { - A = SparseMatrix_symmetrize(A, TRUE); - } - if (!(*comps_ptr)) - *comps_ptr = N_GNEW(m + 1, int); - - *ncomp = 0; - (*comps_ptr)[0] = 0; - for (i = 0; i < m; i++) { - if (i == 0 || mask[i] < 0) { - SparseMatrix_level_sets(A, i, &nlevel, &levelset_ptr, - &levelset, &mask, FALSE); - if (i == 0) - *comps = levelset; - nn = levelset_ptr[nlevel]; - levelset += nn; - (*comps_ptr)[(*ncomp) + 1] = (*comps_ptr)[(*ncomp)] + nn; - (*ncomp)++; - } - - } - if (A != A0) - SparseMatrix_delete(A); - free(mask); -} - -int SparseMatrix_pseudo_diameter(SparseMatrix A0, int root, - int aggressive, int *end1, int *end2, - int *connectedQ) -{ - /* assume unit edge length, unsymmetric matrix ill be symmetrized */ - SparseMatrix A = A0; - int m = A->m, i; - int nlevel; - int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL; - int nlevel0 = 0; - int roots[5], iroots, enda, endb; - - if (!SparseMatrix_is_symmetric(A, TRUE)) { - A = SparseMatrix_symmetrize(A, TRUE); - } - - assert(SparseMatrix_is_symmetric(A, TRUE)); - - SparseMatrix_level_sets(A, root, &nlevel, &levelset_ptr, &levelset, - &mask, TRUE); - *connectedQ = (levelset_ptr[nlevel] == m); - while (nlevel0 < nlevel) { - nlevel0 = nlevel; - root = levelset[levelset_ptr[nlevel] - 1]; - SparseMatrix_level_sets(A, root, &nlevel, &levelset_ptr, &levelset, - &mask, TRUE); - } - *end1 = levelset[0]; - *end2 = levelset[levelset_ptr[nlevel] - 1]; - - if (aggressive) { - nlevel0 = nlevel; - iroots = 0; - for (i = levelset_ptr[nlevel]; - i < MIN(levelset_ptr[nlevel + 1], levelset_ptr[nlevel] + 5); - i++) { - iroots++; - roots[i - levelset_ptr[nlevel]] = levelset[levelset_ptr[i]]; - } - for (i = 0; i < iroots; i++) { - root = roots[i]; - nlevel = - SparseMatrix_pseudo_diameter(A, root, FALSE, &enda, &endb, - connectedQ); - if (nlevel > nlevel0) { - nlevel0 = nlevel; - *end1 = enda; - *end2 = endb; - } - } - } - - free(levelset_ptr); - free(levelset); - free(mask); - if (A != A0) - SparseMatrix_delete(A); - return nlevel0; -} - -int SparseMatrix_pseudo_diameter_only(SparseMatrix A) -{ - int end1, end2, connectedQ; - return SparseMatrix_pseudo_diameter(A, 0, FALSE, &end1, &end2, - &connectedQ); -} - -int SparseMatrix_connectedQ(SparseMatrix A0) -{ - int root = 0, nlevel, *levelset_ptr = NULL, *levelset = NULL, *mask = - NULL, connected; - SparseMatrix A = A0; - - if (!SparseMatrix_is_symmetric(A, TRUE)) { - A = SparseMatrix_symmetrize(A, TRUE); - } - - SparseMatrix_level_sets(A, root, &nlevel, &levelset_ptr, &levelset, - &mask, TRUE); - connected = (levelset_ptr[nlevel] == A->m); - - free(levelset_ptr); - free(levelset); - free(mask); - if (A != A0) - SparseMatrix_delete(A); - - return connected; -} - - -void SparseMatrix_decompose_to_supervariables(SparseMatrix A, - int *ncluster, int **cluster, - int **clusterp) -{ - /* nodes for a super variable if they share exactly the same neighbors. This is know as modules in graph theory. - We work on columns only and columns with the same pattern are grouped as a super variable - */ - int *ia = A->ia, *ja = A->ja, n = A->n, m = A->m; - int *super = NULL, *nsuper = NULL, i, j, *mask = - NULL, isup, *newmap, isuper; - - super = N_GNEW(n, int); - nsuper = N_GNEW(n + 1, int); - mask = N_GNEW(n, int); - newmap = N_GNEW(n, int); - nsuper++; - - isup = 0; - for (i = 0; i < n; i++) - super[i] = isup; /* every node belongs to super variable 0 by default */ - nsuper[0] = n; - for (i = 0; i < n; i++) - mask[i] = -1; - isup++; - - for (i = 0; i < m; i++) { -#ifdef DEBUG_PRINT1 - printf("\n"); - printf("doing row %d-----\n", i + 1); -#endif - for (j = ia[i]; j < ia[i + 1]; j++) { - isuper = super[ja[j]]; - nsuper[isuper]--; /* those entries will move to a different super vars */ - } - for (j = ia[i]; j < ia[i + 1]; j++) { - isuper = super[ja[j]]; - if (mask[isuper] < i) { - mask[isuper] = i; - if (nsuper[isuper] == 0) { /* all nodes in the isuper group exist in this row */ -#ifdef DEBUG_PRINT1 - printf("node %d keep super node id %d\n", ja[j] + 1, - isuper + 1); -#endif - nsuper[isuper] = 1; - newmap[isuper] = isuper; - } else { - newmap[isuper] = isup; - nsuper[isup] = 1; -#ifdef DEBUG_PRINT1 - printf("make node %d into supernode %d\n", ja[j] + 1, - isup + 1); -#endif - super[ja[j]] = isup++; - } - } else { -#ifdef DEBUG_PRINT1 - printf("node %d join super node %d\n", ja[j] + 1, - newmap[isuper] + 1); -#endif - super[ja[j]] = newmap[isuper]; - nsuper[newmap[isuper]]++; - } - } -#ifdef DEBUG_PRINT1 - printf("nsuper="); - for (j = 0; j < isup; j++) - printf("(%d,%d),", j + 1, nsuper[j]); - printf("\n"); -#endif - } -#ifdef DEBUG_PRINT1 - for (i = 0; i < n; i++) { - printf("node %d is in supernode %d\n", i, super[i]); - } -#endif -#ifdef PRINT - fprintf(stderr, "n = %d, nsup = %d\n", n, isup); -#endif - /* now accumulate super nodes */ - nsuper--; - nsuper[0] = 0; - for (i = 0; i < isup; i++) - nsuper[i + 1] += nsuper[i]; - - *cluster = newmap; - for (i = 0; i < n; i++) { - isuper = super[i]; - (*cluster)[nsuper[isuper]++] = i; - } - for (i = isup; i > 0; i--) - nsuper[i] = nsuper[i - 1]; - nsuper[0] = 0; - *clusterp = nsuper; - *ncluster = isup; - -#ifdef PRINT - for (i = 0; i < *ncluster; i++) { - printf("{"); - for (j = (*clusterp)[i]; j < (*clusterp)[i + 1]; j++) { - printf("%d, ", (*cluster)[j]); - } - printf("},"); - } - printf("\n"); -#endif - - free(mask); - free(super); - -} diff --git a/lib/sfdpgen/SparseMatrix.h b/lib/sfdpgen/SparseMatrix.h deleted file mode 100644 index 28cf7d3e6..000000000 --- a/lib/sfdpgen/SparseMatrix.h +++ /dev/null @@ -1,131 +0,0 @@ -/* 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 SPARSEMATRIX_H -#define SPARSEMATRIX_H - -#include -#include - -#define SYMMETRY_EPSILON 0.0000001 -enum { FORMAT_CSC, FORMAT_CSR, FORMAT_COORD }; -enum { UNMASKED = -10 }; -enum { MATRIX_PATTERN_SYMMETRIC = 1 << 0, MATRIX_SYMMETRIC = - 1 << 1, MATRIX_SKEW = 1 << 2, MATRIX_HERMITIAN = - 1 << 3, MATRIX_UNDIRECTED = 1 << 4 -}; - -struct SparseMatrix_s { - int m; /* row dimension */ - int n; /* column dimension */ - int nz; /* The actual length used is nz, for CSR/CSC matrix this is the same as ia[n] */ - int nzmax; /* the current length of ja and a (if exists) allocated. */ - int type; /* whether it is real/complex matrix, or pattern only */ - int *ia; /* row pointer for CSR format, or row indices for coordinate format. 0-based */ - int *ja; /* column indices. 0-based */ - void *a; /* entry values. If NULL, pattern matrix */ - int format; /* whether it is CSR, CSC, COORD. By default it is in CSR format */ - int property; /* pattern_symmetric/symmetric/skew/hermitian */ -}; - -typedef struct SparseMatrix_s* SparseMatrix; - -enum { MATRIX_TYPE_REAL = 1 << 0, MATRIX_TYPE_COMPLEX = - 1 << 1, MATRIX_TYPE_INTEGER = 1 << 2, MATRIX_TYPE_PATTERN = - 1 << 3, MATRIX_TYPE_UNKNOWN = 1 << 4 -}; - -SparseMatrix SparseMatrix_new(int m, int n, int nz, int type, int format); - -SparseMatrix SparseMatrix_from_coordinate_format(SparseMatrix A); - -SparseMatrix SparseMatrix_from_coordinate_arrays(int nz, int m, int n, - int *irn, int *jcn, - void *val, int type); - -void SparseMatrix_print(char *, SparseMatrix A); /*print to stdout in Mathematica format */ - -void SparseMatrix_export(FILE * f, SparseMatrix A); /* export into MM format except the header */ - -SparseMatrix SparseMatrix_import_binary(char *name); - -void SparseMatrix_export_binary(char *name, SparseMatrix A, int *flag); - -void SparseMatrix_delete(SparseMatrix A); - -SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B); -SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B); -SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A); -SparseMatrix SparseMatrix_coordinate_form_add_entries(SparseMatrix A, - int nentries, - int *irn, int *jcn, - void *val); -int SparseMatrix_is_symmetric(SparseMatrix A, - int test_pattern_symmetry_only); -SparseMatrix SparseMatrix_transpose(SparseMatrix A); -SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, - int pattern_symmetric_only); -void SparseMatrix_multiply_vector(SparseMatrix A, real * v, real ** res, int transposed); /* if v = NULL, v is assumed to be {1,1,...,1} */ -SparseMatrix SparseMatrix_remove_diagonal(SparseMatrix A); -SparseMatrix SparseMatrix_remove_upper(SparseMatrix A); /* remove diag and upper diag */ -SparseMatrix SparseMatrix_divide_row_by_degree(SparseMatrix A); -SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A); /* symmetric, all entries to 1, diaginal removed */ -SparseMatrix SparseMatrix_normalize_to_rowsum1(SparseMatrix A); /* for real only! */ -void SparseMatrix_multiply_dense(SparseMatrix A, int ATranspose, - real * v, int vTransposed, real ** res, - int res_transpose, int dim); -SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, double (*fun) (double x)); /* for real only! */ -SparseMatrix SparseMatrix_copy(SparseMatrix A); -int SparseMatrix_has_diagonal(SparseMatrix A); -SparseMatrix SparseMatrix_normalize_by_row(SparseMatrix A); /* divide by max of each row */ -SparseMatrix SparseMatrix_crop(SparseMatrix A, real epsilon); /*remove any entry <= epsilon */ -SparseMatrix SparseMatrix_scaled_by_vector(SparseMatrix A, real * v, - int apply_to_row); -SparseMatrix SparseMatrix_multiply_by_scaler(SparseMatrix A, real s); -SparseMatrix SparseMatrix_make_undirected(SparseMatrix A); /* make it strictly low diag only, and set flag to undirected */ -int SparseMatrix_connectedQ(SparseMatrix A); -int SparseMatrix_pseudo_diameter_only(SparseMatrix A); -int SparseMatrix_pseudo_diameter(SparseMatrix A0, int root, int aggressive, int *end1, int *end2, int *connectedQ); /* assume unit edge length, unsymmetric matrix ill be symmetrized */ -void SparseMatrix_level_sets(SparseMatrix A, int root, int *nlevel, - int **levelset_ptr, int **levelset, - int **mask, int reintialize_mask); -void SparseMatrix_weakly_connected_components(SparseMatrix A0, - int *ncomp, int **comps, - int **comps_ptr); -void SparseMatrix_decompose_to_supervariables(SparseMatrix A, - int *ncluster, int **cluster, - int **clusterp); - - -#define set_flag(a, flag) ((a)=((a)|(flag))) -#define test_flag(a, flag) ((a)&(flag)) -#define clear_flag(a, flag) ((a) &=(~(flag))) - -#define SparseMatrix_set_undirected(A) set_flag((A)->property, MATRIX_UNDIRECTED) -#define SparseMatrix_set_symmetric(A) set_flag((A)->property, MATRIX_SYMMETRIC) -#define SparseMatrix_set_pattern_symmetric(A) set_flag((A)->property, MATRIX_PATTERN_SYMMETRIC) -#define SparseMatrix_set_skew(A) set_flag((A)->property, MATRIX_SKEW) -#define SparseMatrix_set_hemitian(A) set_flag((A)->property, MATRIX_HERMITIAN) - -#define SparseMatrix_known_undirected(A) test_flag((A)->property, MATRIX_UNDIRECTED) -#define SparseMatrix_known_symmetric(A) test_flag((A)->property, MATRIX_SYMMETRIC) -#define SparseMatrix_known_strucural_symmetric(A) test_flag((A)->property, MATRIX_PATTERN_SYMMETRIC) -#define SparseMatrix_known_skew(A) test_flag((A)->property, MATRIX_SKEW) -#define SparseMatrix_known_hemitian(A) test_flag((A)->property, MATRIX_HERMITIAN) - - - - -#endif diff --git a/lib/sfdpgen/call_tri.c b/lib/sfdpgen/call_tri.c deleted file mode 100644 index 904f545ee..000000000 --- a/lib/sfdpgen/call_tri.c +++ /dev/null @@ -1,111 +0,0 @@ -/* 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 * -**********************************************************/ - -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif - -#include "SparseMatrix.h" -#include "logic.h" -#include "memory.h" -#include "delaunay.h" - -SparseMatrix call_tri(int n, int dim, real * x) -{ - real one = 1; - int i, ii, jj; - SparseMatrix A; - SparseMatrix B; - int* edgelist = NULL; - real* xv = N_GNEW(n, real); - real* yv = N_GNEW(n, real); - int numberofedges; - - for (i = 0; i < n; i++) { - xv[i] = x[i * 2]; - yv[i] = x[i * 2 + 1]; - } - - if (n > 2) { - edgelist = delaunay_tri (xv, yv, n, &numberofedges); - } else { - numberofedges = 0; - } - - A = SparseMatrix_new(n, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD); - for (i = 0; i < numberofedges; i++) { - ii = edgelist[i * 2]; - jj = edgelist[i * 2 + 1]; - SparseMatrix_coordinate_form_add_entries(A, 1, &(ii), &(jj), &one); - } - if (n == 2) { /* if two points, add edge i->j */ - ii = 0; - jj = 1; - SparseMatrix_coordinate_form_add_entries(A, 1, &(ii), &(jj), &one); - } - for (i = 0; i < n; i++) { - SparseMatrix_coordinate_form_add_entries(A, 1, &i, &i, &one); - } - B = SparseMatrix_from_coordinate_format(A); - B = SparseMatrix_symmetrize(B, FALSE); - SparseMatrix_delete(A); - - free (edgelist); - free (xv); - free (yv); - return B; -} - -SparseMatrix call_tri2(int n, int dim, real * xx) -{ - real *x, *y; - v_data *delaunay; - int i, j; - SparseMatrix A; - SparseMatrix B; - real one = 1; - x = N_GNEW(n, real); - y = N_GNEW(n, real); - - for (i = 0; i < n; i++) { - x[i] = xx[dim * i]; - y[i] = xx[dim * i + 1]; - } - - delaunay = UG_graph(x, y, n, 0); - - A = SparseMatrix_new(n, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD); - - for (i = 0; i < n; i++) { - for (j = 1; j < delaunay[i].nedges; j++) { - SparseMatrix_coordinate_form_add_entries(A, 1, &i, - &(delaunay[i]. - edges[j]), &one); - } - } - for (i = 0; i < n; i++) { - SparseMatrix_coordinate_form_add_entries(A, 1, &i, &i, &one); - } - B = SparseMatrix_from_coordinate_format(A); - B = SparseMatrix_symmetrize(B, FALSE); - SparseMatrix_delete(A); - - free (x); - free (y); - freeGraph (delaunay); - - return B; - -} diff --git a/lib/sfdpgen/call_tri.h b/lib/sfdpgen/call_tri.h deleted file mode 100644 index 639dca16b..000000000 --- a/lib/sfdpgen/call_tri.h +++ /dev/null @@ -1,22 +0,0 @@ -/* 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 CALL_TRI_H -#define CALL_TRI_H - -SparseMatrix call_tri(int n, int dim, real * x); -SparseMatrix call_tri2(int n, int dim, real * x); - -#endif diff --git a/lib/sfdpgen/layout_similarity.c b/lib/sfdpgen/layout_similarity.c deleted file mode 100644 index 477377fb2..000000000 --- a/lib/sfdpgen/layout_similarity.c +++ /dev/null @@ -1,197 +0,0 @@ -#include -#include -#include -#include "general.h" -#include "SparseMatrix.h" -#include "call_tri.h" - -static double boundingbox_area(int n, double *x){ - double xmax, xmin, ymax, ymin; - int i; - xmax = xmin = x[0]; - ymax = ymin = x[1]; - for (i = 0; i < n; i++){ - xmax = MAX(xmax, x[2*i]); - xmin = MIN(xmin, x[2*i]); - ymax = MAX(ymax, x[2*i+1]); - ymin = MIN(ymin, x[2*i+1]); - } - return (xmax-xmin)*(ymax-ymin); - -} - - -#define MINDIST 0.000000001 - -real distance_cropped(real *x, int dim, int i, int j){ - int k; - real dist = 0.; - for (k = 0; k < dim; k++) dist += (x[i*dim+k] - x[j*dim + k])*(x[i*dim+k] - x[j*dim + k]); - dist = sqrt(dist); - return MAX(dist, MINDIST); -} - -real distance(real *x, int dim, int i, int j){ - int k; - real dist = 0.; - for (k = 0; k < dim; k++) dist += (x[i*dim+k] - x[j*dim + k])*(x[i*dim+k] - x[j*dim + k]); - dist = sqrt(dist); - return dist; -} - -static real get_r(int n, real *x, real *y, real theta, real *r, real *p, real *q){ - /* find the best scaling factor r* such that - f[theta, r] = Sum[(x[i,1] + p - r(Sin[theta] y[i,1] + Cos[theta] y[i,2]))^2+(x[i,2] + q - r(Cos[theta] y[i,1] -Sin[theta] y[i,2]))^2, i = 1, n] is minimized. Return the value f[theta, r*], - also return r. - */ - real top = 0, bot = 0, c, s, f = 0, yy, yx, xxavg = 0, xyavg = 0, yxavg = 0, yyavg = 0; - int i; - - for (i = 0; i < n; i++){ - xxavg += x[2*i]; - xyavg += x[2*i+1]; - } - xxavg /= n; - xyavg /= n; - - for (i = 0; i < n; i++){ - s = sin(theta); - c = cos(theta); - yx = c*y[2*i] - s*y[2*i+1]; - yy = s*y[2*i] + c*y[2*i+1]; - top += (x[2*i] - xxavg)*yx + (x[2*i+1] - xyavg)*yy; - bot += yx*yx+yy*yy; - yxavg += yx; - yyavg += yy; - } - *r = top/(bot - yxavg*yxavg/n - yyavg*yyavg/n); - - *p = yxavg/n*(*r)-xxavg; - *q = yyavg/n*(*r)-xyavg; - - for (i = 0; i < n; i++){ - s = sin(theta); - c = cos(theta); - yx = c*y[2*i] - s*y[2*i+1]; - yy = s*y[2*i] + c*y[2*i+1]; - f += (x[2*i] + *p - (*r)*yx)*(x[2*i] + *p - (*r)*yx)+(x[2*i+1] + *q - (*r)*yy)*(x[2*i+1] + *q - (*r)*yy); - } - return f; -} - -static real dispacement(int n, real *x, real *y, real *rmin, real *thetamin, real *pmin, real *qmin){ - /* find the distance between two sets of cordinates by finding a suitable - rotation and scaling to minimize the distance */ - real dist = 0, theta, normx = 0, r, d, p, q; - int i; - - *thetamin = 0.; - dist = get_r(n, x, y, 0., &r, &p, &q); - *rmin = r; *pmin = p; *qmin = q; - - for (i = 0; i < 180; i++){ - theta = 3.1415626*i/180; - d = get_r(n, x, y, theta, &r, &p, &q); - if (d < dist){ - dist = d; - *rmin = r; *pmin = p; *qmin = q; - *thetamin = theta; - } - } - - - for (i = 0; i < n; i++) normx += x[2*i]*x[2*i]+x[2*i+1]*x[2*i+1]; - return sqrt(dist/normx); -} - -int main(int argc, char *argv[]) -{ - - - char *infile; - - SparseMatrix B = NULL; - int dim = 2, n = 0, i, *ia, *ja, j, k, nz = 0; - real *x, *y, mean, dev, ratio, *z, r, theta, p, q, xx; - - FILE *fp; - - - if (argc < 2){ - fprintf(stderr,"Usage: similarity graph layout1 layout2\n"); - } - - infile = argv[1]; - fp = fopen(infile,"r"); - while (fscanf(fp,"%lf %lf\n",&xx, &xx) == 2) n++; - fclose(fp); - - infile = argv[1]; - fp = fopen(infile,"r"); - x = MALLOC(sizeof(real)*n*2); - for (i = 0; i < n; i++){ - fscanf(fp,"%lf %lf\n",&(x[2*i]), &(x[2*i+1])); - } - fclose(fp); - - - infile = argv[2]; - fp = fopen(infile,"r"); - y = MALLOC(sizeof(real)*n*2); - nz = 0; - for (i = 0; i < n; i++){ - if (fscanf(fp,"%lf %lf\n",&(y[2*i]), &(y[2*i+1])) != 2) goto ERROR; - } - fclose(fp); - - B = call_tri(n, 2, x); - - z = MALLOC(sizeof(real)*B->nz); - ia = B->ia; ja = B->ja; - nz = 0; - for (i = 0; i < n; i++){ - for (j = ia[i]; j < ia[i+1]; j++){ - k = ja[j]; - if (k == i) continue; - z[nz++] = distance(y,dim,i,k)/distance_cropped(x,dim,i,k); - } - } - - /* mean */ - mean = 0; - for (i = 0; i < nz; i++) mean += z[i]; - mean /= nz; - - /* deviation*/ - dev = 0; - for (i = 0; i < nz; i++) { - dev += (z[i]-mean)*(z[i]-mean); - } - dev /= nz; - dev = sqrt(dev); - - /* bounding box area */ - ratio = boundingbox_area(n, y); - /*/MAX(boundingbox_area(n, x), 0.001);*/ - - - fprintf(stderr, "mean = %f std = %f disimilarity = %f area = %f badness = %f displacement = %f\n", - mean, dev, dev/mean, ratio, (dev/mean+1)*ratio, dispacement(n, x, y, &r, &theta, &p, &q)); - fprintf(stderr, "theta = %f scaling = %f, shift = {%f, %f}\n",theta, 1/r, p, q); - - - printf("%.2f %.2f %.2f\n",dev/mean, dispacement(n, x, y, &r, &theta, &p, &q),ratio/1000000.); - - SparseMatrix_delete(B); - FREE(x); - FREE(y); - FREE(z); - - - return 0; - - ERROR: - printf("- - -\n"); - return 0; -} - diff --git a/lib/sfdpgen/overlap.c b/lib/sfdpgen/overlap.c deleted file mode 100644 index ac6067354..000000000 --- a/lib/sfdpgen/overlap.c +++ /dev/null @@ -1,545 +0,0 @@ -/* 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 "SparseMatrix.h" -#include "overlap.h" -#include "call_tri.h" -#include "red_black_tree.h" -#include "types.h" -#include "memory.h" -#include "globals.h" -#include - -#define MACHINEACC 1.0e-16 - -static void ideal_distance_avoid_overlap(int dim, SparseMatrix A, - real * x, real * width, - real * ideal_distance, - real * tmax, real * tmin) -{ - /* if (x1>x2 && y1 > y2) we want either x1 + t (x1-x2) - x2 > (width1+width2), or y1 + t (y1-y2) - y2 > (height1+height2), - hence t = MAX(expandmin, MIN(expandmax, (width1+width2)/(x1-x2) - 1, (height1+height2)/(y1-y2) - 1)), and - new ideal distance = (1+t) old_distance. t can be negative sometimes. - The result ideal distance is set to negative if the edge needs shrinking - */ - int i, j, jj; - int *ia = A->ia, *ja = A->ja; - real dist, dx, dy, wx, wy, t; - real expandmax = 1.5, expandmin = 1; - - *tmax = 0; - *tmin = 1.e10; - assert(SparseMatrix_is_symmetric(A, FALSE)); - for (i = 0; i < A->m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - jj = ja[j]; - if (jj == i) - continue; - dist = distance(x, dim, i, jj); - dx = ABS(x[i * dim] - x[jj * dim]); - dy = ABS(x[i * dim + 1] - x[jj * dim + 1]); - wx = width[i * dim] + width[jj * dim]; - wy = width[i * dim + 1] + width[jj * dim + 1]; - if (dx < MACHINEACC * wx && dy < MACHINEACC * wy) { - ideal_distance[j] = sqrt(wx * wx + wy * wy); - *tmax = 2; - } else { - if (dx < MACHINEACC * wx) { - t = wy / dy; - } else if (dy < MACHINEACC * wy) { - t = wx / dx; - } else { - t = MIN(wx / dx, wy / dy); - } - if (t > 1) - t = MAX(t, 1.001); /* no point in things like t = 1.00000001 as this slow down convergence */ - *tmax = MAX(*tmax, t); - *tmin = MIN(*tmin, t); - t = MIN(expandmax, t); - t = MAX(expandmin, t); - if (t > 1) { - ideal_distance[j] = t * dist; - } else { - ideal_distance[j] = -t * dist; - } - } - - } - } - return; -} - -#define collide(i,j) ((ABS(x[(i)*dim] - x[(j)*dim]) < width[(i)*dim]+width[(j)*dim]) || (ABS(x[(i)*dim+1] - x[(j)*dim+1]) < width[(i)*dim+1]+width[(j)*dim+1])) - -enum { INTV_OPEN, INTV_CLOSE }; - -struct scan_point_struct { - int node; - real x; - int status; -}; - -typedef struct scan_point_struct scan_point; - - -static int comp_scan_points(const void *p, const void *q) -{ - scan_point *pp = (scan_point *) p; - scan_point *qq = (scan_point *) q; - if (pp->x > qq->x) { - return 1; - } else if (pp->x < qq->x) { - return -1; - } else { - if (pp->node > qq->node) { - return 1; - } else if (pp->node < qq->node) { - return -1; - } - return 0; - } - return 0; -} - - -void NodeDest(void *a) -{ - /* free((int*)a); */ -} - - - -int NodeComp(const void *a, const void *b) -{ - return comp_scan_points(a, b); - -} - -void NodePrint(const void *a) -{ - scan_point *aa; - - aa = (scan_point *) a; - fprintf(stderr, "node {%d, %f, %d}\n", aa->node, aa->x, aa->status); - -} - -void InfoPrint(void *a) -{ - ; -} - -void InfoDest(void *a) -{ - ; -} - -static SparseMatrix get_overlap_graph(int dim, int n, real * x, - real * width) -{ - scan_point *scanpointsx, *scanpointsy; - int i, k, neighbor; - SparseMatrix A = NULL; - SparseMatrix B = NULL; - rb_red_blk_node *newNode, *newNode0; - rb_red_blk_tree *treey; - real one = 1; - - A = SparseMatrix_new(n, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD); - - scanpointsx = N_GNEW(2 * n, scan_point); - for (i = 0; i < n; i++) { - scanpointsx[2 * i].node = i; - scanpointsx[2 * i].x = x[i * dim] - width[i * dim]; - scanpointsx[2 * i].status = INTV_OPEN; - scanpointsx[2 * i + 1].node = i + n; - scanpointsx[2 * i + 1].x = x[i * dim] + width[i * dim]; - scanpointsx[2 * i + 1].status = INTV_CLOSE; - } - qsort(scanpointsx, 2 * n, sizeof(scan_point), comp_scan_points); - - scanpointsy = N_GNEW(2 * n, scan_point); - for (i = 0; i < n; i++) { - scanpointsy[i].node = i; - scanpointsy[i].x = x[i * dim + 1] - width[i * dim + 1]; - scanpointsy[i].status = INTV_OPEN; - scanpointsy[i + n].node = i; - scanpointsy[i + n].x = x[i * dim + 1] + width[i * dim + 1]; - scanpointsy[i + n].status = INTV_CLOSE; - } - - - treey = - RBTreeCreate(NodeComp, NodeDest, InfoDest, NodePrint, InfoPrint); - - for (i = 0; i < 2 * n; i++) { -#ifdef DEBUG_RBTREE - fprintf(stderr, " k = %d node = %d x====%f\n", - (scanpointsx[i].node) % n, (scanpointsx[i].node), - (scanpointsx[i].x)); -#endif - - k = (scanpointsx[i].node) % n; - - - if (scanpointsx[i].status == INTV_OPEN) { -#ifdef DEBUG_RBTREE - fprintf(stderr, "inserting..."); - treey->PrintKey(&(scanpointsy[k])); -#endif - - RBTreeInsert(treey, &(scanpointsy[k]), NULL); /* add both open and close int for y */ - -#ifdef DEBUG_RBTREE - fprintf(stderr, "inserting2..."); - treey->PrintKey(&(scanpointsy[k + n])); -#endif - - RBTreeInsert(treey, &(scanpointsy[k + n]), NULL); - } else { - assert(scanpointsx[i].node >= n); - - newNode = newNode0 = - RBExactQuery(treey, &(scanpointsy[k + n])); - -#ifdef DEBUG_RBTREE - fprintf(stderr, "poping..%d....", scanpointsy[k + n].node); - treey->PrintKey(newNode->key); -#endif - - assert(treey->nil != newNode); - while (((newNode = - TreePredecessor(treey, newNode)) != treey->nil) - && ((scan_point *) newNode->key)->node != k) { - neighbor = (((scan_point *) newNode->key)->node) % n; - A = SparseMatrix_coordinate_form_add_entries(A, 1, - &neighbor, &k, - &one); -#ifdef DEBUG_RBTREE - fprintf(stderr, "%d %d\n", k, neighbor); -#endif - - } - -#ifdef DEBUG_RBTREE - fprintf(stderr, "deleting..."); - treey->PrintKey(newNode0->key); -#endif - - RBDelete(treey, newNode0); - if (newNode != treey->nil && newNode != newNode0) { - -#ifdef DEBUG_RBTREE - fprintf(stderr, "deleting2..."); - treey->PrintKey(newNode->key) -#endif - RBDelete(treey, newNode); - } - } - } - - free(scanpointsx); - free(scanpointsy); - RBTreeDestroy(treey); - - B = SparseMatrix_from_coordinate_format(A); - SparseMatrix_delete(A); - A = SparseMatrix_symmetrize(B, FALSE); - SparseMatrix_delete(B); - if (Verbose) - fprintf(stderr, "found %d clashes\n", A->nz); - return A; -} - - - -/* ============================== label overlap smoother ==================*/ - -#define OverlapSmoother_s StressMajorizationSmoother_s - -OverlapSmoother OverlapSmoother_new(SparseMatrix A, int dim, - real lambda0, real * x, real * width, - int include_original_graph, - int neighborhood_only, - real * max_overlap, - real * min_overlap) -{ - OverlapSmoother sm; - int i, j, k, m = A->m, *iw, *jw, *id, *jd, jdiag; - SparseMatrix B; - real *lambda, *d, *w, diag_d, diag_w, dist; - - assert((!A) || SparseMatrix_is_symmetric(A, FALSE)); - - sm = GNEW(struct OverlapSmoother_s); - lambda = sm->lambda = N_GNEW(m, real); - for (i = 0; i < m; i++) - sm->lambda[i] = lambda0; - - B = call_tri(m, dim, x); - - if (!neighborhood_only) { - SparseMatrix C; - SparseMatrix D; - C = get_overlap_graph(dim, A->m, x, width); - D = SparseMatrix_add(B, C); - SparseMatrix_delete(B); - SparseMatrix_delete(C); - B = D; - } - if (include_original_graph) { - sm->Lw = SparseMatrix_add(A, B); - SparseMatrix_delete(B); - } else { - sm->Lw = B; - } - sm->Lwd = SparseMatrix_copy(sm->Lw); - -#ifdef DEBUG - { - FILE *fp; - fp = fopen("/tmp/111", "w"); - export_embedding(fp, dim, sm->Lwd, x, NULL); - fclose(fp); - } -#endif - - if (!(sm->Lw) || !(sm->Lwd)) { - OverlapSmoother_delete(sm); - return NULL; - } - - assert((sm->Lwd)->type == MATRIX_TYPE_REAL); - - ideal_distance_avoid_overlap(dim, sm->Lwd, x, width, - (real *) (sm->Lwd->a), max_overlap, - min_overlap); - - /* no overlap at all! */ - if (*max_overlap < 1) { - if (Verbose) - fprintf(stderr, - " no overlap (overlap = %f), rescale to shrink\n", - *max_overlap - 1); - for (i = 0; i < dim * m; i++) { - x[i] *= (*max_overlap); - } - *max_overlap = 1; - goto RETURN; - } - - 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; - } - if (d[j] > 0) { /* those edges that needs expansion */ - w[j] = 100 / d[j] / d[j]; - } else { /* those that needs shrinking is set to negative in ideal_distance_avoid_overlap */ - w[j] = 1 / d[j] / d[j]; - d[j] = -d[j]; - } - dist = d[j]; - diag_w += w[j]; - d[j] = w[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; - } - RETURN: - return sm; -} - -void OverlapSmoother_delete(OverlapSmoother sm) -{ - - StressMajorizationSmoother_delete(sm); - -} - -void OverlapSmoother_smooth(OverlapSmoother sm, int dim, real * x) -{ - int maxit_sm = 1; /* only using 1 iteration of stress majorization - is found to give better results and save time! */ - StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm); -#ifdef DEBUG - { - FILE *fp; - fp = fopen("/tmp/222", "w"); - export_embedding(fp, dim, sm->Lwd, x, NULL); - fclose(fp); - } -#endif -} - -/*================================= end OverlapSmoother =============*/ - -static void scale_to_edge_length(int dim, SparseMatrix A, real * x, - real avg_label_size) -{ - real dist; - int i; - - if (!A) - return; - dist = average_edge_length(A, dim, x); - if (Verbose) - fprintf(stderr, "avg edge len=%f avg_label-size= %f\n", dist, - avg_label_size); - - - dist = avg_label_size / MAX(dist, MACHINEACC); - - for (i = 0; i < dim * A->m; i++) - x[i] *= dist; -} - -static void print_bounding_box(int n, int dim, real * x) -{ - real *xmin, *xmax; - int i, k; - - xmin = N_GNEW(dim, real); - xmax = N_GNEW(dim, real); - - for (i = 0; i < dim; i++) - xmin[i] = xmax[i] = x[i]; - - for (i = 0; i < n; i++) { - for (k = 0; k < dim; k++) { - xmin[k] = MIN(xmin[k], x[i * dim + k]); - xmax[k] = MAX(xmax[k], x[i * dim + k]); - } - } - fprintf(stderr, "bounding box = \n"); - for (i = 0; i < dim; i++) - fprintf(stderr, "{%f,%f}, ", xmin[i], xmax[i]); - fprintf(stderr, "\n"); - - free(xmin); - free(xmax); -} - -int remove_overlap(int dim, SparseMatrix A, real * x, - real * label_sizes, int ntry) -{ - real lambda = 0.00; - OverlapSmoother sm; - int include_original_graph = 0, i; - real avg_label_size; - real max_overlap = 0, min_overlap = 999; - int neighborhood_only = TRUE; - int flag = 0; - -#ifdef DEBUG - clock_t cpu = clock(); -#endif - - if (!label_sizes) - return flag; - - avg_label_size = 0; - for (i = 0; i < A->m; i++) - avg_label_size += label_sizes[i * dim] + label_sizes[i * dim + 1]; - /* for (i = 0; i < A->m; i++) avg_label_size += 2*MAX(label_sizes[i*dim],label_sizes[i*dim+1]); */ - avg_label_size /= A->m; - - scale_to_edge_length(dim, A, x, 4 * avg_label_size); - - if (!ntry) - return flag; - - -#ifdef DEBUG - _statistics[0] = _statistics[1] = 0.; - { - FILE *fp; - fp = fopen("x1", "w"); - for (i = 0; i < A->m; i++) { - fprintf(fp, "%f %f\n", x[i * 2], x[i * 2 + 1]); - } - fclose(fp); - } -#endif - - for (i = 0; i < ntry; i++) { - if (Verbose) - print_bounding_box(A->m, dim, x); - sm = OverlapSmoother_new(A, dim, lambda, x, label_sizes, - include_original_graph, neighborhood_only, - &max_overlap, &min_overlap); - if (Verbose) - fprintf(stderr, - "overlap removal neighbors only?= %d iter -- %d, overlap factor = %g underlap factor = %g\n", - neighborhood_only, i, max_overlap - 1, min_overlap); - if (max_overlap <= 1) { - OverlapSmoother_delete(sm); - if (neighborhood_only == FALSE) { - break; - } else { - neighborhood_only = FALSE; - continue; - } - } - - OverlapSmoother_smooth(sm, dim, x); - OverlapSmoother_delete(sm); - } - -#ifdef DEBUG - fprintf(stderr, - " number of cg iter = %f, number of stress majorization iter = %f number of overlap removal try = %d\n", - _statistics[0], _statistics[1], i - 1); - - { - FILE *fp; - fp = fopen("x2", "w"); - for (i = 0; i < A->m; i++) { - fprintf(fp, "%f %f\n", x[i * 2], x[i * 2 + 1]); - } - fclose(fp); - } -#endif - -#ifdef DEBUG - { - FILE *fp; - fp = fopen("/tmp/m", "w"); - export_embedding(fp, dim, A, x, label_sizes); - } -#endif -#ifdef DEBUG - fprintf(stderr, "post processing %f\n", - ((real) (clock() - cpu)) / CLOCKS_PER_SEC); -#endif - return flag; -} diff --git a/lib/sfdpgen/overlap.h b/lib/sfdpgen/overlap.h deleted file mode 100644 index b10ea3dd2..000000000 --- a/lib/sfdpgen/overlap.h +++ /dev/null @@ -1,39 +0,0 @@ - -/* 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 OVERLAP_H -#define OVERLAP_H - -#include "post_process.h" - -typedef StressMajorizationSmoother OverlapSmoother; - -void OverlapSmoother_delete(OverlapSmoother sm); - -OverlapSmoother OverlapSmoother_new(SparseMatrix A, int dim, - real lambda0, real * x, real * width, - int include_original_graph, - int neighborhood_only, - real * max_overlap, - real * min_overlap); - -void OverlapSmoother_smooth(OverlapSmoother sm, int dim, real * x); - - -int remove_overlap(int dim, SparseMatrix A, real * x, - real * label_sizes, int ntry); - -#endif diff --git a/lib/sfdpgen/post_process.c b/lib/sfdpgen/post_process.c deleted file mode 100644 index 7462747dd..000000000 --- a/lib/sfdpgen/post_process.c +++ /dev/null @@ -1,756 +0,0 @@ -/* 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 -#include -#include - -#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 = N_GNEW(D->nz, real); - } - 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 /= 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; -} - -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(struct StressMajorizationSmoother_s); - 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; - } - - - for (i = 0; i < m; i++) - mask[i] = -1; - - 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++; - } - } - } - } - - 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); - } - - /* - 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++; - } - } - - 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] = -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; - - 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; - 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; - } - /* 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]; - } - } - - maxit = sqrt((double) m); - res = - SparseMatrix_solve(Lw, dim, x, y, 0.01, maxit, SOLVE_METHOD_CG, - &flag); - - if (flag) - goto RETURN; - - 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); - } -#endif - memcpy(x, y, sizeof(real) * m * dim); - } - -#ifdef DEBUG - _statistics[1] += iter - 1; -#endif - - 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); -} - - -#define TriangleSmoother_s StressMajorizationSmoother_s - -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; - } - - sm = GNEW(struct TriangleSmoother_s); - lambda = sm->lambda = N_GNEW(m, real); - for (i = 0; i < m; i++) - sm->lambda[i] = lambda0; - - if (m > 2) { - if (use_triangularization) { - B = call_tri(m, dim, x); - } else { - B = call_tri2(m, dim, x); - } - } 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]; - - } - - 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; - } - - s = stop / sbot; - for (i = 0; i < iw[m]; i++) - d[i] *= s; - - free(avg_dist); - - return sm; -} - -void TriangleSmoother_delete(TriangleSmoother sm) -{ - - StressMajorizationSmoother_delete(sm); - -} - -void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real * x) -{ - - 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; - - sm = GNEW(struct SpringSmoother_s); - 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; - } - - - for (i = 0; i < m; i++) - mask[i] = -1; - - 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++; - } - } - } - } - - 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++; - } - } - } - 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; -} - - -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; - - spring_electrical_spring_embedding(dim, A, sm->D, sm->ctrl, - node_weights, x, &flag); - assert(!flag); - -} - -/*=============================== end of spring and spring-electrical based smoother =========== */ - -int post_process_smoothing(int dim, SparseMatrix A, - spring_electrical_control ctrl, - real * node_weights, real * x) -{ -#ifdef DEBUG - clock_t cpu = clock(); -#endif - - 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; - - 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 */ - -#ifdef DEBUG - if (Verbose) - fprintf(stderr, "post processing %f\n", - ((real) (clock() - cpu)) / CLOCKS_PER_SEC); -#endif - return flag; -} diff --git a/lib/sfdpgen/post_process.h b/lib/sfdpgen/post_process.h deleted file mode 100644 index 080f9ab46..000000000 --- a/lib/sfdpgen/post_process.h +++ /dev/null @@ -1,76 +0,0 @@ -/* 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 - -struct StressMajorizationSmoother_s { - SparseMatrix Lw; - SparseMatrix Lwd; - real *lambda; -}; - -typedef struct StressMajorizationSmoother_s* 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; - -void TriangleSmoother_delete(TriangleSmoother sm); - -TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, - real lambda, real * x, - int use_triangularization); - -void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real * x); - - - -/*------------------ spring and spring-electrical based smoother */ - -struct SpringSmoother_s { - SparseMatrix D; - spring_electrical_control ctrl; -}; -typedef struct SpringSmoother_s* SpringSmoother; - -SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, - spring_electrical_control ctrl, - real * x); - -void SpringSmoother_delete(SpringSmoother sm); - -void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, - real * node_weights, int dim, real * x); -/*------------------------------------------------------------------*/ - -int post_process_smoothing(int dim, SparseMatrix A, - spring_electrical_control ctrl, - real * node_weights, real * x); - -#endif diff --git a/lib/sfdpgen/red_black_tree.c b/lib/sfdpgen/red_black_tree.c deleted file mode 100644 index 776228b45..000000000 --- a/lib/sfdpgen/red_black_tree.c +++ /dev/null @@ -1,922 +0,0 @@ -/* 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 "red_black_tree.h" -#include -#include - -/* CONVENTIONS: All data structures for red-black trees have the prefix */ -/* "rb_" to prevent name conflicts. */ -/* */ -/* Function names: Each word in a function name begins with */ -/* a capital letter. An example funcntion name is */ -/* CreateRedTree(a,b,c). Furthermore, each function name */ -/* should begin with a capital letter to easily distinguish */ -/* them from variables. */ -/* */ -/* Variable names: Each word in a variable name begins with */ -/* a capital letter EXCEPT the first letter of the variable */ -/* name. For example, int newLongInt. Global variables have */ -/* names beginning with "g". An example of a global */ -/* variable name is gNewtonsConstant. */ - -void Assert(int assertion, char *error); -void *SafeMalloc(size_t size); - -/*--------------------- stack.h --------*/ - -/* CONVENTIONS: All data structures for stacks have the prefix */ -/* "stk_" to prevent name conflicts. */ -/* */ -/* Function names: Each word in a function name begins with */ -/* a capital letter. An example funcntion name is */ -/* CreateRedTree(a,b,c). Furthermore, each function name */ -/* should begin with a capital letter to easily distinguish */ -/* them from variables. */ -/* */ -/* Variable names: Each word in a variable name begins with */ -/* a capital letter EXCEPT the first letter of the variable */ -/* name. For example, int newLongInt. Global variables have */ -/* names beginning with "g". An example of a global */ -/* variable name is gNewtonsConstant. */ - -/* if DATA_TYPE is undefined then stack.h and stack.c will be code for */ -/* stacks of void *, if they are defined then they will be stacks of the */ -/* appropriate data_type */ - -#ifndef DATA_TYPE -#define DATA_TYPE void * -#endif - -typedef struct stk_stack_node { - DATA_TYPE info; - struct stk_stack_node *next; -} stk_stack_node; - -typedef struct stk_stack { - stk_stack_node *top; - stk_stack_node *tail; -} stk_stack; - -/* These functions are all very straightforward and self-commenting so */ -/* I didn't think additional comments would be useful */ -stk_stack *StackJoin(stk_stack * stack1, stk_stack * stack2); -stk_stack *StackCreate(); -void StackPush(stk_stack * theStack, DATA_TYPE newInfoPointer); -void *StackPop(stk_stack * theStack); -int StackNotEmpty(stk_stack *); - -/* comment out the line below to remove all the debugging assertion */ -/* checks from the compiled code. */ -#define DEBUG_ASSERT 1 - -/***********************************************************************/ -/* FUNCTION: RBTreeCreate */ - /**/ -/* INPUTS: All the inputs are names of functions. CompFunc takes to */ -/* void pointers to keys and returns 1 if the first arguement is */ -/* "greater than" the second. DestFunc takes a pointer to a key and */ -/* destroys it in the appropriate manner when the node containing that */ -/* key is deleted. InfoDestFunc is similiar to DestFunc except it */ -/* recieves a pointer to the info of a node and destroys it. */ -/* PrintFunc recieves a pointer to the key of a node and prints it. */ -/* PrintInfo recieves a pointer to the info of a node and prints it. */ -/* If RBTreePrint is never called the print functions don't have to be */ -/* defined and NullFunction can be used. */ - /**/ -/* OUTPUT: This function returns a pointer to the newly created */ -/* red-black tree. */ - /**/ -/* Modifies Input: none */ -/***********************************************************************/ -rb_red_blk_tree * -RBTreeCreate(int (*CompFunc) (const void *, const void *), - void (*DestFunc) (void *), void (*InfoDestFunc) (void *), - void (*PrintFunc) (const void *), void (*PrintInfo) (void *)) -{ - rb_red_blk_tree *newTree; - rb_red_blk_node *temp; - - newTree = (rb_red_blk_tree *) SafeMalloc(sizeof(rb_red_blk_tree)); - newTree->Compare = CompFunc; - newTree->DestroyKey = DestFunc; - newTree->PrintKey = PrintFunc; - newTree->PrintInfo = PrintInfo; - newTree->DestroyInfo = InfoDestFunc; - - /* see the comment in the rb_red_blk_tree structure in red_black_tree.h */ - /* for information on nil and root */ - temp = newTree->nil = - (rb_red_blk_node *) SafeMalloc(sizeof(rb_red_blk_node)); - temp->parent = temp->left = temp->right = temp; - temp->red = 0; - temp->key = 0; - temp = newTree->root = - (rb_red_blk_node *) SafeMalloc(sizeof(rb_red_blk_node)); - temp->parent = temp->left = temp->right = newTree->nil; - temp->key = 0; - temp->red = 0; - return (newTree); -} - -/***********************************************************************/ -/* FUNCTION: LeftRotate */ - /**/ -/* INPUTS: This takes a tree so that it can access the appropriate */ -/* root and nil pointers, and the node to rotate on. */ - /**/ -/* OUTPUT: None */ - /**/ -/* Modifies Input: tree, x */ - /**/ -/* EFFECTS: Rotates as described in _Introduction_To_Algorithms by */ -/* Cormen, Leiserson, Rivest (Chapter 14). Basically this */ -/* makes the parent of x be to the left of x, x the parent of */ -/* its parent before the rotation and fixes other pointers */ -/* accordingly. */ -/***********************************************************************/ -void LeftRotate(rb_red_blk_tree * tree, rb_red_blk_node * x) -{ - rb_red_blk_node *y; - rb_red_blk_node *nil = tree->nil; - - /* I originally wrote this function to use the sentinel for */ - /* nil to avoid checking for nil. However this introduces a */ - /* very subtle bug because sometimes this function modifies */ - /* the parent pointer of nil. This can be a problem if a */ - /* function which calls LeftRotate also uses the nil sentinel */ - /* and expects the nil sentinel's parent pointer to be unchanged */ - /* after calling this function. For example, when RBDeleteFixUP */ - /* calls LeftRotate it expects the parent pointer of nil to be */ - /* unchanged. */ - - y = x->right; - x->right = y->left; - - if (y->left != nil) - y->left->parent = x; /* used to use sentinel here */ - /* and do an unconditional assignment instead of testing for nil */ - - y->parent = x->parent; - - /* instead of checking if x->parent is the root as in the book, we */ - /* count on the root sentinel to implicitly take care of this case */ - if (x == x->parent->left) { - x->parent->left = y; - } else { - x->parent->right = y; - } - y->left = x; - x->parent = y; - -#ifdef DEBUG_ASSERT - Assert(!tree->nil->red, "nil not red in LeftRotate"); -#endif -} - - -/***********************************************************************/ -/* FUNCTION: RighttRotate */ - /**/ -/* INPUTS: This takes a tree so that it can access the appropriate */ -/* root and nil pointers, and the node to rotate on. */ - /**/ -/* OUTPUT: None */ - /**/ -/* Modifies Input?: tree, y */ - /**/ -/* EFFECTS: Rotates as described in _Introduction_To_Algorithms by */ -/* Cormen, Leiserson, Rivest (Chapter 14). Basically this */ -/* makes the parent of x be to the left of x, x the parent of */ -/* its parent before the rotation and fixes other pointers */ -/* accordingly. */ -/***********************************************************************/ -void RightRotate(rb_red_blk_tree * tree, rb_red_blk_node * y) -{ - rb_red_blk_node *x; - rb_red_blk_node *nil = tree->nil; - - /* I originally wrote this function to use the sentinel for */ - /* nil to avoid checking for nil. However this introduces a */ - /* very subtle bug because sometimes this function modifies */ - /* the parent pointer of nil. This can be a problem if a */ - /* function which calls LeftRotate also uses the nil sentinel */ - /* and expects the nil sentinel's parent pointer to be unchanged */ - /* after calling this function. For example, when RBDeleteFixUP */ - /* calls LeftRotate it expects the parent pointer of nil to be */ - /* unchanged. */ - - x = y->left; - y->left = x->right; - - if (nil != x->right) - x->right->parent = y; /*used to use sentinel here */ - /* and do an unconditional assignment instead of testing for nil */ - - /* instead of checking if x->parent is the root as in the book, we */ - /* count on the root sentinel to implicitly take care of this case */ - x->parent = y->parent; - if (y == y->parent->left) { - y->parent->left = x; - } else { - y->parent->right = x; - } - x->right = y; - y->parent = x; - -#ifdef DEBUG_ASSERT - Assert(!tree->nil->red, "nil not red in RightRotate"); -#endif -} - -/***********************************************************************/ -/* FUNCTION: TreeInsertHelp */ - /**/ -/* INPUTS: tree is the tree to insert into and z is the node to insert */ - /**/ -/* OUTPUT: none */ - /**/ -/* Modifies Input: tree, z */ - /**/ -/* EFFECTS: Inserts z into the tree as if it were a regular binary tree */ -/* using the algorithm described in _Introduction_To_Algorithms_ */ -/* by Cormen et al. This funciton is only intended to be called */ -/* by the RBTreeInsert function and not by the user */ -/***********************************************************************/ -void TreeInsertHelp(rb_red_blk_tree * tree, rb_red_blk_node * z) -{ - /* This function should only be called by InsertRBTree (see above) */ - rb_red_blk_node *x; - rb_red_blk_node *y; - rb_red_blk_node *nil = tree->nil; - - z->left = z->right = nil; - y = tree->root; - x = tree->root->left; - while (x != nil) { - y = x; - if (1 == tree->Compare(x->key, z->key)) { /* x.key > z.key */ - x = x->left; - } else { /* x,key <= z.key */ - x = x->right; - } - } - z->parent = y; - if ((y == tree->root) || (1 == tree->Compare(y->key, z->key))) { /* y.key > z.key */ - y->left = z; - } else { - y->right = z; - } - -#ifdef DEBUG_ASSERT - Assert(!tree->nil->red, "nil not red in TreeInsertHelp"); -#endif -} - -/* Before calling Insert RBTree the node x should have its key set */ - -/***********************************************************************/ -/* FUNCTION: RBTreeInsert */ - /**/ -/* INPUTS: tree is the red-black tree to insert a node which has a key */ -/* pointed to by key and info pointed to by info. */ - /**/ -/* OUTPUT: This function returns a pointer to the newly inserted node */ -/* which is guarunteed to be valid until this node is deleted. */ -/* What this means is if another data structure stores this */ -/* pointer then the tree does not need to be searched when this */ -/* is to be deleted. */ - /**/ -/* Modifies Input: tree */ - /**/ -/* EFFECTS: Creates a node node which contains the appropriate key and */ -/* info pointers and inserts it into the tree. */ -/***********************************************************************/ - rb_red_blk_node * RBTreeInsert(rb_red_blk_tree * tree, void *key, - void *info) -{ - rb_red_blk_node *y; - rb_red_blk_node *x; - rb_red_blk_node *newNode; - - x = (rb_red_blk_node *) SafeMalloc(sizeof(rb_red_blk_node)); - x->key = key; - x->info = info; - - TreeInsertHelp(tree, x); - newNode = x; - x->red = 1; - while (x->parent->red) { /* use sentinel instead of checking for root */ - if (x->parent == x->parent->parent->left) { - y = x->parent->parent->right; - if (y->red) { - x->parent->red = 0; - y->red = 0; - x->parent->parent->red = 1; - x = x->parent->parent; - } else { - if (x == x->parent->right) { - x = x->parent; - LeftRotate(tree, x); - } - x->parent->red = 0; - x->parent->parent->red = 1; - RightRotate(tree, x->parent->parent); - } - } else { /* case for x->parent == x->parent->parent->right */ - y = x->parent->parent->left; - if (y->red) { - x->parent->red = 0; - y->red = 0; - x->parent->parent->red = 1; - x = x->parent->parent; - } else { - if (x == x->parent->left) { - x = x->parent; - RightRotate(tree, x); - } - x->parent->red = 0; - x->parent->parent->red = 1; - LeftRotate(tree, x->parent->parent); - } - } - } - tree->root->left->red = 0; - return (newNode); - -#ifdef DEBUG_ASSERT - Assert(!tree->nil->red, "nil not red in RBTreeInsert"); - Assert(!tree->root->red, "root not red in RBTreeInsert"); -#endif -} - -/***********************************************************************/ -/* FUNCTION: TreeSuccessor */ - /**/ -/* INPUTS: tree is the tree in question, and x is the node we want the */ -/* the successor of. */ - /**/ -/* OUTPUT: This function returns the successor of x or NULL if no */ -/* successor exists. */ - /**/ -/* Modifies Input: none */ - /**/ -/* Note: uses the algorithm in _Introduction_To_Algorithms_ */ -/***********************************************************************/ - rb_red_blk_node * TreeSuccessor(rb_red_blk_tree * tree, - rb_red_blk_node * x) -{ - rb_red_blk_node *y; - rb_red_blk_node *nil = tree->nil; - rb_red_blk_node *root = tree->root; - - if (nil != (y = x->right)) { /* assignment to y is intentional */ - while (y->left != nil) { /* returns the minium of the right subtree of x */ - y = y->left; - } - return (y); - } else { - y = x->parent; - while (x == y->right) { /* sentinel used instead of checking for nil */ - x = y; - y = y->parent; - } - if (y == root) - return (nil); - return (y); - } -} - -/***********************************************************************/ -/* FUNCTION: Treepredecessor */ - /**/ -/* INPUTS: tree is the tree in question, and x is the node we want the */ -/* the predecessor of. */ - /**/ -/* OUTPUT: This function returns the predecessor of x or NULL if no */ -/* predecessor exists. */ - /**/ -/* Modifies Input: none */ - /**/ -/* Note: uses the algorithm in _Introduction_To_Algorithms_ */ -/***********************************************************************/ - rb_red_blk_node * TreePredecessor(rb_red_blk_tree * tree, - rb_red_blk_node * x) -{ - rb_red_blk_node *y; - rb_red_blk_node *nil = tree->nil; - rb_red_blk_node *root = tree->root; - - if (nil != (y = x->left)) { /* assignment to y is intentional */ - while (y->right != nil) { /* returns the maximum of the left subtree of x */ - y = y->right; - } - return (y); - } else { - y = x->parent; - while (x == y->left) { - if (y == root) - return (nil); - x = y; - y = y->parent; - } - return (y); - } -} - -/***********************************************************************/ -/* FUNCTION: InorderTreePrint */ - /**/ -/* INPUTS: tree is the tree to print and x is the current inorder node */ - /**/ -/* OUTPUT: none */ - /**/ -/* EFFECTS: This function recursively prints the nodes of the tree */ -/* inorder using the PrintKey and PrintInfo functions. */ - /**/ -/* Modifies Input: none */ - /**/ -/* Note: This function should only be called from RBTreePrint */ -/***********************************************************************/ -void InorderTreePrint(rb_red_blk_tree * tree, rb_red_blk_node * x) -{ - rb_red_blk_node *nil = tree->nil; - rb_red_blk_node *root = tree->root; - if (x != tree->nil) { - InorderTreePrint(tree, x->left); - printf("info="); - tree->PrintInfo(x->info); - printf(" key="); - tree->PrintKey(x->key); - printf(" l->key="); - if (x->left == nil) - printf("NULL"); - else - tree->PrintKey(x->left->key); - printf(" r->key="); - if (x->right == nil) - printf("NULL"); - else - tree->PrintKey(x->right->key); - printf(" p->key="); - if (x->parent == root) - printf("NULL"); - else - tree->PrintKey(x->parent->key); - printf(" red=%i\n", x->red); - InorderTreePrint(tree, x->right); - } -} - - -/***********************************************************************/ -/* FUNCTION: TreeDestHelper */ - /**/ -/* INPUTS: tree is the tree to destroy and x is the current node */ - /**/ -/* OUTPUT: none */ - /**/ -/* EFFECTS: This function recursively destroys the nodes of the tree */ -/* postorder using the DestroyKey and DestroyInfo functions. */ - /**/ -/* Modifies Input: tree, x */ - /**/ -/* Note: This function should only be called by RBTreeDestroy */ -/***********************************************************************/ -void TreeDestHelper(rb_red_blk_tree * tree, rb_red_blk_node * x) -{ - rb_red_blk_node *nil = tree->nil; - if (x != nil) { - TreeDestHelper(tree, x->left); - TreeDestHelper(tree, x->right); - tree->DestroyKey(x->key); - tree->DestroyInfo(x->info); - free(x); - } -} - - -/***********************************************************************/ -/* FUNCTION: RBTreeDestroy */ - /**/ -/* INPUTS: tree is the tree to destroy */ - /**/ -/* OUTPUT: none */ - /**/ -/* EFFECT: Destroys the key and frees memory */ - /**/ -/* Modifies Input: tree */ - /**/ -/***********************************************************************/ -void RBTreeDestroy(rb_red_blk_tree * tree) -{ - TreeDestHelper(tree, tree->root->left); - free(tree->root); - free(tree->nil); - free(tree); -} - - -/***********************************************************************/ -/* FUNCTION: RBTreePrint */ - /**/ -/* INPUTS: tree is the tree to print */ - /**/ -/* OUTPUT: none */ - /**/ -/* EFFECT: This function recursively prints the nodes of the tree */ -/* inorder using the PrintKey and PrintInfo functions. */ - /**/ -/* Modifies Input: none */ - /**/ -/***********************************************************************/ -void RBTreePrint(rb_red_blk_tree * tree) -{ - InorderTreePrint(tree, tree->root->left); -} - - -/***********************************************************************/ -/* FUNCTION: RBExactQuery */ - /**/ -/* INPUTS: tree is the tree to print and q is a pointer to the key */ -/* we are searching for */ - /**/ -/* OUTPUT: returns the a node with key equal to q. If there are */ -/* multiple nodes with key equal to q this function returns */ -/* the one highest in the tree */ - /**/ -/* Modifies Input: none */ - /**/ -/***********************************************************************/ - rb_red_blk_node * RBExactQuery(rb_red_blk_tree * tree, void *q) -{ - rb_red_blk_node *x = tree->root->left; - rb_red_blk_node *nil = tree->nil; - int compVal; - if (x == nil) - return (0); - compVal = tree->Compare(x->key, (int *) q); - while (0 != compVal) { /*assignemnt */ - if (1 == compVal) { /* x->key > q */ - x = x->left; - } else { - x = x->right; - } - if (x == nil) - return (0); - compVal = tree->Compare(x->key, (int *) q); - } - return (x); -} - - -/***********************************************************************/ -/* FUNCTION: RBDeleteFixUp */ - /**/ -/* INPUTS: tree is the tree to fix and x is the child of the spliced */ -/* out node in RBTreeDelete. */ - /**/ -/* OUTPUT: none */ - /**/ -/* EFFECT: Performs rotations and changes colors to restore red-black */ -/* properties after a node is deleted */ - /**/ -/* Modifies Input: tree, x */ - /**/ -/* The algorithm from this function is from _Introduction_To_Algorithms_ */ -/***********************************************************************/ -void RBDeleteFixUp(rb_red_blk_tree * tree, rb_red_blk_node * x) -{ - rb_red_blk_node *root = tree->root->left; - rb_red_blk_node *w; - - while ((!x->red) && (root != x)) { - if (x == x->parent->left) { - w = x->parent->right; - if (w->red) { - w->red = 0; - x->parent->red = 1; - LeftRotate(tree, x->parent); - w = x->parent->right; - } - if ((!w->right->red) && (!w->left->red)) { - w->red = 1; - x = x->parent; - } else { - if (!w->right->red) { - w->left->red = 0; - w->red = 1; - RightRotate(tree, w); - w = x->parent->right; - } - w->red = x->parent->red; - x->parent->red = 0; - w->right->red = 0; - LeftRotate(tree, x->parent); - x = root; /* this is to exit while loop */ - } - } else { /* the code below is has left and right switched from above */ - w = x->parent->left; - if (w->red) { - w->red = 0; - x->parent->red = 1; - RightRotate(tree, x->parent); - w = x->parent->left; - } - if ((!w->right->red) && (!w->left->red)) { - w->red = 1; - x = x->parent; - } else { - if (!w->left->red) { - w->right->red = 0; - w->red = 1; - LeftRotate(tree, w); - w = x->parent->left; - } - w->red = x->parent->red; - x->parent->red = 0; - w->left->red = 0; - RightRotate(tree, x->parent); - x = root; /* this is to exit while loop */ - } - } - } - x->red = 0; - -#ifdef DEBUG_ASSERT - Assert(!tree->nil->red, "nil not black in RBDeleteFixUp"); -#endif -} - - -/***********************************************************************/ -/* FUNCTION: RBDelete */ - /**/ -/* INPUTS: tree is the tree to delete node z from */ - /**/ -/* OUTPUT: none */ - /**/ -/* EFFECT: Deletes z from tree and frees the key and info of z */ -/* using DestoryKey and DestoryInfo. Then calls */ -/* RBDeleteFixUp to restore red-black properties */ - /**/ -/* Modifies Input: tree, z */ - /**/ -/* The algorithm from this function is from _Introduction_To_Algorithms_ */ -/***********************************************************************/ -void RBDelete(rb_red_blk_tree * tree, rb_red_blk_node * z) -{ - rb_red_blk_node *y; - rb_red_blk_node *x; - rb_red_blk_node *nil = tree->nil; - rb_red_blk_node *root = tree->root; - - y = ((z->left == nil) - || (z->right == nil)) ? z : TreeSuccessor(tree, z); - x = (y->left == nil) ? y->right : y->left; - if (root == (x->parent = y->parent)) { /* assignment of y->p to x->p is intentional */ - root->left = x; - } else { - if (y == y->parent->left) { - y->parent->left = x; - } else { - y->parent->right = x; - } - } - if (y != z) { /* y should not be nil in this case */ - -#ifdef DEBUG_ASSERT - Assert((y != tree->nil), "y is nil in RBDelete\n"); -#endif - /* y is the node to splice out and x is its child */ - - if (!(y->red)) - RBDeleteFixUp(tree, x); - - tree->DestroyKey(z->key); - tree->DestroyInfo(z->info); - y->left = z->left; - y->right = z->right; - y->parent = z->parent; - y->red = z->red; - z->left->parent = z->right->parent = y; - if (z == z->parent->left) { - z->parent->left = y; - } else { - z->parent->right = y; - } - free(z); - } else { - tree->DestroyKey(y->key); - tree->DestroyInfo(y->info); - if (!(y->red)) - RBDeleteFixUp(tree, x); - free(y); - } - -#ifdef DEBUG_ASSERT - Assert(!tree->nil->red, "nil not black in RBDelete"); -#endif -} - - -/***********************************************************************/ -/* FUNCTION: RBDEnumerate */ - /**/ -/* INPUTS: tree is the tree to look for keys >= low */ -/* and <= high with respect to the Compare function */ - /**/ -/* OUTPUT: stack containing pointers to the nodes between [low,high] */ - /**/ -/* Modifies Input: none */ -/***********************************************************************/ - stk_stack * RBEnumerate(rb_red_blk_tree * tree, void *low, void *high) -{ - stk_stack *enumResultStack; - rb_red_blk_node *nil = tree->nil; - rb_red_blk_node *x = tree->root->left; - rb_red_blk_node *lastBest = nil; - - enumResultStack = StackCreate(); - while (nil != x) { - if (1 == (tree->Compare(x->key, high))) { /* x->key > high */ - x = x->left; - } else { - lastBest = x; - x = x->right; - } - } - while ((lastBest != nil) && (1 != tree->Compare(low, lastBest->key))) { - StackPush(enumResultStack, lastBest); - lastBest = TreePredecessor(tree, lastBest); - } - return (enumResultStack); -} - - - - - - - - -/***********************************************************************/ -/* FUNCTION: void Assert(int assertion, char* error) */ - /**/ -/* INPUTS: assertion should be a predicated that the programmer */ -/* assumes to be true. If this assumption is not true the message */ -/* error is printed and the program exits. */ - /**/ -/* OUTPUT: None. */ - /**/ -/* Modifies input: none */ - /**/ -/* Note: If DEBUG_ASSERT is not defined then assertions should not */ -/* be in use as they will slow down the code. Therefore the */ -/* compiler will complain if an assertion is used when */ -/* DEBUG_ASSERT is undefined. */ -/***********************************************************************/ -void Assert(int assertion, char *error) -{ - if (!assertion) { - printf("Assertion Failed: %s\n", error); - exit(-1); - } -} - - -/*------------------------ misc.c -------------*/ - -/***********************************************************************/ -/* FUNCTION: SafeMalloc */ - /**/ -/* INPUTS: size is the size to malloc */ - /**/ -/* OUTPUT: returns pointer to allocated memory if succesful */ - /**/ -/* EFFECT: mallocs new memory. If malloc fails, prints error message */ -/* and terminates program. */ - /**/ -/* Modifies Input: none */ - /**/ -/***********************************************************************/ -void *SafeMalloc(size_t size) -{ - void *result; - - if ((result = malloc(size))) { /* assignment intentional */ - return (result); - } else { - printf("memory overflow: malloc failed in SafeMalloc."); - printf(" Exiting Program.\n"); - exit(-1); - return (0); - } -} - -/* NullFunction does nothing it is included so that it can be passed */ -/* as a function to RBTreeCreate when no other suitable function has */ -/* been defined */ - - -void NullFunction(void *junk) -{; -} - - - -/*----------------- stack.c ------------ */ - -int StackNotEmpty(stk_stack * theStack) -{ - return (theStack ? (theStack->top != NULL) : 0); -} - -stk_stack *StackJoin(stk_stack * stack1, stk_stack * stack2) -{ - if (!stack1->tail) { - free(stack1); - return (stack2); - } else { - stack1->tail->next = stack2->top; - stack1->tail = stack2->tail; - free(stack2); - return (stack1); - } -} - -stk_stack *StackCreate() -{ - stk_stack *newStack; - - newStack = (stk_stack *) SafeMalloc(sizeof(stk_stack)); - newStack->top = newStack->tail = NULL; - return (newStack); -} - - -void StackPush(stk_stack * theStack, DATA_TYPE newInfoPointer) -{ - stk_stack_node *newNode; - - if (!theStack->top) { - newNode = (stk_stack_node *) SafeMalloc(sizeof(stk_stack_node)); - newNode->info = newInfoPointer; - newNode->next = theStack->top; - theStack->top = newNode; - theStack->tail = newNode; - } else { - newNode = (stk_stack_node *) SafeMalloc(sizeof(stk_stack_node)); - newNode->info = newInfoPointer; - newNode->next = theStack->top; - theStack->top = newNode; - } - -} - -DATA_TYPE StackPop(stk_stack * theStack) -{ - DATA_TYPE popInfo; - stk_stack_node *oldNode; - - if (theStack->top) { - popInfo = theStack->top->info; - oldNode = theStack->top; - theStack->top = theStack->top->next; - free(oldNode); - if (!theStack->top) - theStack->tail = NULL; - } else { - popInfo = NULL; - } - return (popInfo); -} - -void StackDestroy(stk_stack * theStack, void DestFunc(void *a)) -{ - stk_stack_node *x = theStack->top; - stk_stack_node *y; - - if (theStack) { - while (x) { - y = x->next; - DestFunc(x->info); - free(x); - x = y; - } - free(theStack); - } -} diff --git a/lib/sfdpgen/red_black_tree.h b/lib/sfdpgen/red_black_tree.h deleted file mode 100644 index 85c83c9d2..000000000 --- a/lib/sfdpgen/red_black_tree.h +++ /dev/null @@ -1,76 +0,0 @@ -/* 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 RED_BLACK_TREE_H -#define RED_BLACK_TREE_H - -/* CONVENTIONS: All data structures for red-black trees have the prefix */ -/* "rb_" to prevent name conflicts. */ -/* */ -/* Function names: Each word in a function name begins with */ -/* a capital letter. An example funcntion name is */ -/* CreateRedTree(a,b,c). Furthermore, each function name */ -/* should begin with a capital letter to easily distinguish */ -/* them from variables. */ -/* */ -/* Variable names: Each word in a variable name begins with */ -/* a capital letter EXCEPT the first letter of the variable */ -/* name. For example, int newLongInt. Global variables have */ -/* names beginning with "g". An example of a global */ -/* variable name is gNewtonsConstant. */ - -typedef struct rb_red_blk_node { - void *key; - void *info; - int red; /* if red=0 then the node is black */ - struct rb_red_blk_node *left; - struct rb_red_blk_node *right; - struct rb_red_blk_node *parent; -} rb_red_blk_node; - - -/* Compare(a,b) should return 1 if *a > *b, -1 if *a < *b, and 0 otherwise */ -/* Destroy(a) takes a pointer to whatever key might be and frees it accordingly */ -typedef struct rb_red_blk_tree { - int (*Compare) (const void *a, const void *b); - void (*DestroyKey) (void *a); - void (*DestroyInfo) (void *a); - void (*PrintKey) (const void *a); - void (*PrintInfo) (void *a); - /* A sentinel is used for root and for nil. These sentinels are */ - /* created when RBTreeCreate is caled. root->left should always */ - /* point to the node which is the root of the tree. nil points to a */ - /* node which should always be black but has aribtrary children and */ - /* parent and no key or info. The point of using these sentinels is so */ - /* that the root and nil nodes do not require special cases in the code */ - rb_red_blk_node *root; - rb_red_blk_node *nil; -} rb_red_blk_tree; - -rb_red_blk_tree *RBTreeCreate(int (*CompFunc) (const void *, const void *), - void (*DestFunc) (void *), - void (*InfoDestFunc) (void *), - void (*PrintFunc) (const void *), - void (*PrintInfo) (void *)); -rb_red_blk_node *RBTreeInsert(rb_red_blk_tree *, void *key, void *info); -void RBTreePrint(rb_red_blk_tree *); -void RBDelete(rb_red_blk_tree *, rb_red_blk_node *); -void RBTreeDestroy(rb_red_blk_tree *); -rb_red_blk_node *TreePredecessor(rb_red_blk_tree *, rb_red_blk_node *); -rb_red_blk_node *TreeSuccessor(rb_red_blk_tree *, rb_red_blk_node *); -rb_red_blk_node *RBExactQuery(rb_red_blk_tree *, void *); -void NullFunction(void *); - -#endif diff --git a/lib/sfdpgen/sfdp.h b/lib/sfdpgen/sfdp.h deleted file mode 100644 index 7adb735c2..000000000 --- a/lib/sfdpgen/sfdp.h +++ /dev/null @@ -1,29 +0,0 @@ -/* 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 SFDP_H -#define SFDP_H - -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif - -#include "render.h" - -void sfdp_layout (graph_t * g); -void sfdp_cleanup (graph_t * g); -int fdpAdjust (graph_t * g); - -#endif diff --git a/lib/sfdpgen/sfdpinit.c b/lib/sfdpgen/sfdpinit.c deleted file mode 100644 index 10b06e600..000000000 --- a/lib/sfdpgen/sfdpinit.c +++ /dev/null @@ -1,382 +0,0 @@ -/* 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 * -**********************************************************/ - - -#ifdef HAVE_CONFIG_H -#include "config.h" -#endif -#ifdef HAVE_LIMITS_H -#include -#else -#ifdef HAVE_VALUES_H -#include -#endif -#endif -#include -#include -#include -#include -#include -#include -#include -#include - -#ifdef DEBUG -double _statistics[10]; -#endif - -static void sfdp_init_node(node_t * n) -{ - common_init_node(n); - ND_pos(n) = N_NEW(GD_ndim(n->graph), double); - gv_nodesize(n, GD_flip(n->graph)); -} - -static void sfdp_init_edge(edge_t * e) -{ - common_init_edge(e); -} - -static void sfdp_init_node_edge(graph_t * g) -{ - node_t *n; - edge_t *e; - int nnodes = agnnodes(g); - attrsym_t *N_pos = agfindattr(g->proto->n, "pos"); - - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - sfdp_init_node(n); - user_pos(N_pos, NULL, n, nnodes); - } - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - for (e = agfstout(g, n); e; e = agnxtout(g, e)) - sfdp_init_edge(e); - } -} - -static void sfdp_init_graph(Agraph_t * g) -{ - setEdgeType(g, ET_LINE); - g->u.ndim = late_int(g, agfindattr(g, "dim"), 2, 2); - Ndim = g->u.ndim = MIN(g->u.ndim, MAXDIM); - - sfdp_init_node_edge(g); -} - -/* getPos: - */ -static real *getPos(Agraph_t * g) -{ - Agnode_t *n; - real *pos = N_NEW(Ndim * agnnodes(g), real); - int ix, i; - - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - i = ND_id(n); - if (hasPos(n)) { - for (ix = 0; ix < Ndim; ix++) { - pos[i * Ndim + ix] = ND_pos(n)[ix]; - } - } - } - - return pos; -} - -/* getSizes: - * Set up array of half sizes in inches. - */ -static real *getSizes(Agraph_t * g, pointf pad) -{ - Agnode_t *n; - real *sizes = N_GNEW(2 * agnnodes(g), real); - int i; - - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - i = ND_id(n); - sizes[i * 2] = ND_width(n) * .5 + pad.x; - sizes[i * 2 + 1] = ND_height(n) * .5 + pad.y; - } - - return sizes; -} - -/* makeMatrix: - * Assumes g is connected and simple, i.e., we can have a->b and b->a - * but not a->b and a->b - */ -static SparseMatrix makeMatrix(Agraph_t * g, int dim) -{ - SparseMatrix A = 0; - Agnode_t *n; - Agedge_t *e; - Agsym_t *sym; - int nnodes; - int nedges; - int i, row; - int *I; - int *J; - real *val; - real v; - int type = MATRIX_TYPE_REAL; - - if (!g) - return NULL; - nnodes = agnnodes(g); - nedges = agnedges(g); - - /* Assign node ids */ - i = 0; - for (n = agfstnode(g); n; n = agnxtnode(g, n)) - ND_id(n) = i++; - - I = N_GNEW(nedges, int); - J = N_GNEW(nedges, int); - val = N_GNEW(nedges, real); - - sym = agfindattr(g->proto->e, "wt"); - i = 0; - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - row = ND_id(n); - for (e = agfstout(g, n); e; e = agnxtout(g, e)) { - I[i] = row; - J[i] = ND_id(e->head); - if (!sym || (sscanf(agxget(e, sym->index), "%lf", &v) != 1)) - v = 1; - val[i] = v; - i++; - } - } - - A = SparseMatrix_from_coordinate_arrays(nedges, nnodes, nnodes, I, J, - val, type); - - free(I); - free(J); - free(val); - - return A; -} - -int -fdpAdjust (graph_t* g) -{ - SparseMatrix A = makeMatrix(g, Ndim); - real *sizes; - real *pos = N_NEW(Ndim * agnnodes(g), real); - Agnode_t *n; - int flag, i; - expand_t sep = sepFactor(g); - pointf pad; - - if (sep.doAdd) { - pad.x = PS2INCH(sep.x); - pad.y = PS2INCH(sep.y); - } else { - pad.x = PS2INCH(DFLT_MARGIN); - pad.y = PS2INCH(DFLT_MARGIN); - } - sizes = getSizes(g, pad); - - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - real* npos = pos + (Ndim * ND_id(n)); - for (i = 0; i < Ndim; i++) { - npos[i] = ND_pos(n)[i]; - } - } - - flag = remove_overlap(Ndim, A, pos, sizes, 3); - - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - real *npos = pos + (Ndim * ND_id(n)); - for (i = 0; i < Ndim; i++) { - ND_pos(n)[i] = npos[i]; - } - } - - free(sizes); - free(pos); - SparseMatrix_delete (A); - - return flag; -} - -static void sfdpLayout(graph_t * g, spring_electrical_control ctrl, - pointf pad) -{ - real *sizes; - real *pos; - Agnode_t *n; - int flag, i; - - SparseMatrix A = makeMatrix(g, Ndim); - if (ctrl->overlap) - sizes = getSizes(g, pad); - else - sizes = NULL; - pos = getPos(g); - - flag = multilevel_spring_electrical_embedding(Ndim, A, ctrl, NULL, sizes, pos); - - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - real *npos = pos + (Ndim * ND_id(n)); - for (i = 0; i < Ndim; i++) { - ND_pos(n)[i] = npos[i]; - } - } - - free(sizes); - free(pos); - SparseMatrix_delete (A); -} - -static smooth_t -late_smooth (graph_t* g, Agsym_t* sym, smooth_t dflt) -{ - char* s; - int v; - smooth_t rv; - - if (!sym) return dflt; - s = agxget (g, sym->index); - if (isdigit(*s)) { -#if HAVE_GTS - if ((v = atoi (s)) <= SMOOTHING_RNG) -#else - if ((v = atoi (s)) <= SMOOTHING_SPRING) -#endif - rv = (smooth_t)v; - else - rv = dflt; - } - else if (isalpha(*s)) { - if (!strcasecmp(s, "avg_dist")) - rv = SMOOTHING_STRESS_MAJORIZATION_AVG_DIST; - else if (!strcasecmp(s, "graph_dist")) - rv = SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST; - else if (!strcasecmp(s, "none")) - rv = SMOOTHING_NONE; - else if (!strcasecmp(s, "power_dist")) - rv = SMOOTHING_STRESS_MAJORIZATION_POWER_DIST; -#if HAVE_GTS - else if (!strcasecmp(s, "rng")) - rv = SMOOTHING_RNG; -#endif - else if (!strcasecmp(s, "spring")) - rv = SMOOTHING_SPRING; -#if HAVE_GTS - else if (!strcasecmp(s, "triangle")) - rv = SMOOTHING_TRIANGLE; -#endif - else - rv = dflt; - } - else - rv = dflt; - return rv; -} - -/* tuneControl: - * Use user values to reset control - * - * Possible parameters: - * ctrl->use_node_weights - */ -static void -tuneControl (graph_t* g, spring_electrical_control ctrl) -{ - ctrl->p = -1.0*late_double(g, agfindattr(g, "K"), -AUTOP, 0.0); - ctrl->multilevels = late_int(g, agfindattr(g, "maxiter"), INT_MAX, 0); - ctrl->smoothing = late_smooth(g, agfindattr(g, "smooth"), SMOOTHING_NONE); -} - -void sfdp_layout(graph_t * g) -{ - sfdp_init_graph(g); - - if (agnnodes(g)) { - Agraph_t **ccs; - Agraph_t *sg; - int ncc; - int i; - expand_t sep; - pointf pad; - spring_electrical_control ctrl = spring_electrical_control_new(); - - tuneControl (g, ctrl); - -#if HAVE_GTS - if (!agget(g, "overlap") || (graphAdjustMode(g)->mode == AM_FDP)) { - ctrl->overlap = 3; /* FIX: User control? */ - sep = sepFactor(g); - if (sep.doAdd) { - pad.x = PS2INCH(sep.x); - pad.y = PS2INCH(sep.y); - } else { - pad.x = PS2INCH(DFLT_MARGIN); - pad.y = PS2INCH(DFLT_MARGIN); - } - } -#endif - - ccs = ccomps(g, &ncc, 0); - if (ncc == 1) { - sfdpLayout(g, ctrl, pad); - adjustNodes(g); - spline_edges(g); - } else { - pack_info pinfo; - pack_mode pmode = getPackMode(g, l_node); - - for (i = 0; i < ncc; i++) { - sg = ccs[i]; - nodeInduce(sg); - sfdpLayout(g, ctrl, pad); - adjustNodes(sg); - } - spline_edges(g); - pinfo.margin = getPack(g, CL_OFFSET, CL_OFFSET); - pinfo.doSplines = 1; - pinfo.mode = pmode; - pinfo.fixed = 0; - packSubgraphs(ncc, ccs, g, &pinfo); - } - for (i = 0; i < ncc; i++) { - agdelete(g, ccs[i]); - } - free(ccs); - spring_electrical_control_delete(ctrl); - } - - dotneato_postprocess(g); -} - -static void sfdp_cleanup_graph(graph_t * g) -{ -} - -void sfdp_cleanup(graph_t * g) -{ - node_t *n; - edge_t *e; - - for (n = agfstnode(g); n; n = agnxtnode(g, n)) { - for (e = agfstout(g, n); e; e = agnxtout(g, e)) { - gv_cleanup_edge(e); - } - gv_cleanup_node(n); - } - sfdp_cleanup_graph(g); -} diff --git a/lib/sfdpgen/sfdpinternal.h b/lib/sfdpgen/sfdpinternal.h deleted file mode 100644 index 1fc359349..000000000 --- a/lib/sfdpgen/sfdpinternal.h +++ /dev/null @@ -1,27 +0,0 @@ -/* 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 SFDPINTERNAL_H -#define SFDPINTERNAL_H - -#include - -#ifdef DEBUG -extern double _statistics[10]; -#endif - -typedef double real; - -#endif diff --git a/lib/sfdpgen/sparse_solve.c b/lib/sfdpgen/sparse_solve.c deleted file mode 100644 index 854374c12..000000000 --- a/lib/sfdpgen/sparse_solve.c +++ /dev/null @@ -1,230 +0,0 @@ -/* 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 "sparse_solve.h" -#include "memory.h" -#include "logic.h" -#include "math.h" -#include "assert.h" -#include "string.h" - -static real *vector_subtract_to(int n, real * x, real * y) -{ - /* y = x-y */ - int i; - for (i = 0; i < n; i++) - y[i] = x[i] - y[i]; - return y; -} - -static real *vector_saxpy(int n, real * x, real * y, real beta) -{ - /* y = x+beta*y */ - int i; - for (i = 0; i < n; i++) - y[i] = x[i] + beta * y[i]; - return y; -} - -static real *vector_saxpy2(int n, real * x, real * y, real beta) -{ - /* x = x+beta*y */ - int i; - for (i = 0; i < n; i++) - x[i] = x[i] + beta * y[i]; - return x; -} - -static real vector_product(int n, real *x, real *y){ - real res = 0; - int i; - for (i = 0; i < n; i++) res += x[i]*y[i]; - return res; -} - -real *Operator_matmul_apply(Operator o, real * x, real * y) -{ - SparseMatrix A = (SparseMatrix) o->data; - SparseMatrix_multiply_vector(A, x, &y, FALSE); - return y; -} - -Operator Operator_matmul_new(SparseMatrix A) -{ - Operator o = GNEW(struct Operator_s); - o->data = (void *) A; - o->Operator_apply = Operator_matmul_apply; - return o; -} - -void Operator_matmul_delete(Operator o) -{ -} - -real *Operator_diag_precon_apply(Operator o, real * x, real * y) -{ - int i, m; - real *diag = (real *) o->data; - m = (int) diag[0]; - diag++; - for (i = 0; i < m; i++) - y[i] = x[i] * diag[i]; - return y; -} - -Operator Operator_diag_precon_new(SparseMatrix A) -{ - Operator o; - real *diag; - int i, j, m = A->m, *ia = A->ia, *ja = A->ja; - real *a = (real *) A->a; - - assert(A->type == MATRIX_TYPE_REAL); - - assert(a); - - o = GNEW(struct Operator_s); - o->data = N_GNEW(A->m + 1, real); - diag = (real *) o->data; - - diag[0] = m; - diag++; - for (i = 0; i < m; i++) { - diag[i] = 1.; - for (j = ia[i]; j < ia[i + 1]; j++) { - if (i == ja[j] && ABS(a[j]) > 0) - diag[i] = 1. / a[j]; - } - } - - o->Operator_apply = Operator_diag_precon_apply; - - return o; -} - -void Operator_diag_precon_delete(Operator o) -{ - free(o->data); -} - -static real conjugate_gradient(Operator A, Operator precon, int n, real * x, - real * rhs, real tol, int maxit, int *flag) -{ - real *z, *r, *p, *q, res = 10 * tol, alpha; - real rho = 1.0e20, rho_old = 1, res0, beta; - op_apply_fn Ax = A->Operator_apply; - op_apply_fn Minvx = precon->Operator_apply; - int iter = 0; - - z = N_GNEW(n, real); - r = N_GNEW(n, real); - p = N_GNEW(n, real); - q = N_GNEW(n, real); - - r = Ax(A, x, r); - r = vector_subtract_to(n, rhs, r); - - res0 = res = sqrt(vector_product(n, r, r)) / n; -#ifdef DEBUG_PRINT - if (Verbose && 0) { - fprintf(stderr, " cg iter = %d, residual = %g\n", iter, res); - } -#endif - - while ((iter++) < maxit && res > tol * res0) { - z = Minvx(precon, r, z); - rho = vector_product(n, r, z); - - if (iter > 1) { - beta = rho / rho_old; - p = vector_saxpy(n, z, p, beta); - } else { - memcpy(p, z, sizeof(real) * n); - } - - q = Ax(A, p, q); - - alpha = rho / vector_product(n, p, q); - - x = vector_saxpy2(n, x, p, alpha); - r = vector_saxpy2(n, r, q, -alpha); - - res = sqrt(vector_product(n, r, r)) / n; - -#ifdef DEBUG_PRINT - if (Verbose && 0) { - fprintf(stderr, " cg iter = %d, residual = %g\n", iter, res); - } -#endif - - - - rho_old = rho; - } -#ifdef DEBUG - _statistics[0] += iter - 1; -#endif - -#ifdef DEBUG_PRINT - if (Verbose && 0) { - fprintf(stderr, " cg iter = %d, residual = %g\n", iter, res); - } -#endif - return res; -} - -real SparseMatrix_solve(SparseMatrix A, int dim, real * x0, real * rhs, - real tol, int maxit, int method, int *flag) -{ - Operator Ax; - Operator precond; - real *x, *b, res = 0; - int n = A->m, k, i; - - *flag = 0; - - switch (method) { - case SOLVE_METHOD_CG: - Ax = Operator_matmul_new(A); - precond = Operator_diag_precon_new(A); - - x = N_GNEW(n, real); - b = N_GNEW(n, real); - for (k = 0; k < dim; k++) { - for (i = 0; i < n; i++) { - x[i] = x0[i * dim + k]; - b[i] = rhs[i * dim + k]; - } - - res += - conjugate_gradient(Ax, precond, n, x, b, tol, maxit, flag); - for (i = 0; i < n; i++) { - rhs[i * dim + k] = x[i]; - } - } - Operator_matmul_delete(Ax); - Operator_diag_precon_delete(precond); - free(x); - free(b); - - break; - default: - assert(0); - break; - - } - return res; -} - diff --git a/lib/sfdpgen/sparse_solve.h b/lib/sfdpgen/sparse_solve.h deleted file mode 100644 index dfe502aee..000000000 --- a/lib/sfdpgen/sparse_solve.h +++ /dev/null @@ -1,36 +0,0 @@ -/* 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 SPARSE_SOLVER_H -#define SPARSE_SOLVER_H - -#include "SparseMatrix.h" - -enum { SOLVE_METHOD_CG }; - -typedef struct Operator_s* Operator; - -typedef real *(*op_apply_fn) (Operator o, real * in, real * out); - -struct Operator_s { - void *data; - op_apply_fn Operator_apply; -}; - -real SparseMatrix_solve(SparseMatrix A, int dim, real * x0, real * rhs, - real tol, int maxit, int method, int *flag); - -#endif diff --git a/lib/sfdpgen/spring_electrical.c b/lib/sfdpgen/spring_electrical.c deleted file mode 100644 index ae2621c8a..000000000 --- a/lib/sfdpgen/spring_electrical.c +++ /dev/null @@ -1,1340 +0,0 @@ -/* 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 "spring_electrical.h" -#include "QuadTree.h" -#include "Multilevel.h" -#include "post_process.h" -#include "overlap.h" -#include "types.h" -#include "memory.h" -#include "arith.h" -#include "logic.h" -#include "math.h" -#include "globals.h" -#include -#include - -#define drand() (rand()/(real) RAND_MAX) - -spring_electrical_control spring_electrical_control_new() -{ - spring_electrical_control ctrl; - ctrl = GNEW(struct spring_electrical_control_s); - ctrl->p = AUTOP; /*a negativve number default to -1. repulsive force = dist^p */ - ctrl->random_start = TRUE; /* whether to apply SE from a random layout, or from exisiting layout */ - ctrl->K = -1; /* the natural distance. If K < 0, K will be set to the average distance of an edge */ - ctrl->C = 0.2; /* another parameter. f_a(i,j) = C*dist(i,j)^2/K * d_ij, f_r(i,j) = K^(3-p)/dist(i,j)^(-p). By default C = 0.2. */ - ctrl->multilevels = FALSE; /* if <=1, single level */ - ctrl->quadtree_size = 45; /* cut off size above which quadtree approximation is used */ - ctrl->max_qtree_level = 10; /* max level of quadtree */ - ctrl->bh = 0.6; /* Barnes-Hutt constant, if width(snode)/dist[i,snode] < bh, treat snode as a supernode. */ - ctrl->tol = 0.001; /* minimum different between two subsequence config before terminating. ||x-xold||_infinity < tol/K */ - ctrl->maxiter = 500; - ctrl->cool = 0.90; /* default 0.9 */ - ctrl->step = 0.1; - ctrl->adaptive_cooling = TRUE; - ctrl->random_seed = 123; - ctrl->beautify_leaves = FALSE; - ctrl->use_node_weights = FALSE; - ctrl->smoothing = SMOOTHING_NONE; - ctrl->overlap = 0; - - return ctrl; -} - -void spring_electrical_control_delete(spring_electrical_control ctrl) -{ - free(ctrl); -} - - - -enum { MAX_I = 20, OPT_UP = 1, OPT_DOWN = -1, OPT_INIT = 0 }; -typedef struct { - int i; - real work[MAX_I + 1]; - int direction; -} oned_optimizer; - -static void oned_optimizer_delete(oned_optimizer * opt) -{ - free(opt); -} - -static oned_optimizer *oned_optimizer_new(int i) -{ - oned_optimizer *opt; - opt = GNEW(oned_optimizer); - opt->i = i; - opt->direction = OPT_INIT; - return opt; -} - -static void oned_optimizer_train(oned_optimizer * opt, real work) -{ - int i = opt->i; - - assert(i >= 0); - opt->work[i] = work; - if (opt->direction == OPT_INIT) { - if (opt->i == MAX_I) { - opt->direction = OPT_DOWN; - opt->i = opt->i - 1; - } else { - opt->direction = OPT_UP; - opt->i = MIN(MAX_I, opt->i + 1); - } - } else if (opt->direction == OPT_UP) { - /* fprintf(stderr, "{current_level, prev_level} = {%d,%d}, {work, work_prev} = {%f,%f}",i,i-1,opt->work[i], opt->work[i-1]); */ - assert(i >= 1); - if (opt->work[i] < opt->work[i - 1] && opt->i < MAX_I) { - /* fprintf(stderr, "keep going up to level %d\n",opt->i+1); */ - opt->i = MIN(MAX_I, opt->i + 1); - } else { - /* fprintf(stderr, "going down to level %d\n",opt->i-1); */ - (opt->i)--; - opt->direction = OPT_DOWN; - } - } else { - assert(i < MAX_I); - /* fprintf(stderr, "{current_level, prev_level} = {%d,%d}, {work, work_prev} = {%f,%f}",i,i+1,opt->work[i], opt->work[i+1]); */ - if (opt->work[i] < opt->work[i + 1] && opt->i > 0) { - /* fprintf(stderr, "keep going down to level %d\n",opt->i-1); */ - opt->i = MAX(0, opt->i - 1); - } else { - /* fprintf(stderr, "keep up to level %d\n",opt->i+1); */ - (opt->i)++; - opt->direction = OPT_UP; - } - } -} - -static int oned_optimizer_get(oned_optimizer * opt) -{ - return opt->i; -} - - -real average_edge_length(SparseMatrix A, int dim, real * coord) -{ - real dist = 0, d; - int *ia = A->ia, *ja = A->ja, i, j, k; - assert(SparseMatrix_is_symmetric(A, TRUE)); - - if (ia[A->m] == 0) - return 1; - for (i = 0; i < A->m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - d = 0; - for (k = 0; k < dim; k++) { - d += (coord[dim * i + k] - - coord[dim * ja[j]]) * (coord[dim * i + k] - - coord[dim * ja[j]]); - } - dist += sqrt(d); - } - } - return dist / ia[A->m]; -} - -real distance_cropped(real * x, int dim, int i, int j) -{ - int k; - real dist = 0.; - for (k = 0; k < dim; k++) - dist += - (x[i * dim + k] - x[j * dim + k]) * (x[i * dim + k] - - x[j * dim + k]); - dist = sqrt(dist); - return MAX(dist, MINDIST); -} - -real distance(real * x, int dim, int i, int j) -{ - int k; - real dist = 0.; - for (k = 0; k < dim; k++) - dist += - (x[i * dim + k] - x[j * dim + k]) * (x[i * dim + k] - - x[j * dim + k]); - dist = sqrt(dist); - return dist; -} - -#ifdef ENERGY -static real spring_electrical_energy(int dim, SparseMatrix A, real * x, - real p, real CRK, real KP) -{ - /* 1. Grad[||x-y||^k,x] = k||x-y||^(k-1)*0.5*(x-y)/||x-y|| = k/2*||x-y||^(k-2) (x-y) - which should equal to -force (force = -gradient), - hence energy for force ||x-y||^m (y-x) is ||x-y||^(m+2)*2/(m+2) where m != 2 - 2. Grad[Log[||x-y||],x] = 1/||x-y||*0.5*(x-y)/||x-y|| = 0.5*(x-y)/||x-y||^2, - hence the energy to give force ||x-y||^-2 (x-y) is -2*Log[||x-y||] - - */ - int i, j, k, *ia = A->ia, *ja = A->ja, n = A->m; - real energy = 0, dist; - - for (i = 0; i < n; i++) { - /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */ - for (j = ia[i]; j < ia[i + 1]; j++) { - if (ja[j] == i) - continue; - dist = distance(x, dim, i, ja[j]); - energy += CRK * pow(dist, 3.) * 2. / 3.; - } - - /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */ - for (j = 0; j < n; j++) { - if (j == i) - continue; - dist = distance_cropped(x, dim, i, j); - for (k = 0; k < dim; k++) { - if (p == -1) { - energy += -KP * 2 * log(dist); - } else { - energy += -KP * pow(dist, p + 1) * 2 / (p + 1); - } - } - } - } - return energy; -} - -#endif - -void export_embedding(FILE * fp, int dim, SparseMatrix A, real * x, - real * width) -{ - int i, j, k, *ia = A->ia, *ja = A->ja; - int ne = 0; - if (dim == 2) { - fprintf(fp, "Graphics[{GrayLevel[0.5],Line[{"); - } else { - fprintf(fp, "Graphics3D[{GrayLevel[0.5],Line[{"); - } - for (i = 0; i < A->m; i++) { - for (j = ia[i]; j < ia[i + 1]; j++) { - if (ja[j] == i) - continue; - ne++; - if (ne > 1) - fprintf(fp, ","); - fprintf(fp, "{{"); - for (k = 0; k < dim; k++) { - if (k > 0) - fprintf(fp, ","); - fprintf(fp, "%f", x[i * dim + k]); - } - fprintf(fp, "},{"); - for (k = 0; k < dim; k++) { - if (k > 0) - fprintf(fp, ","); - fprintf(fp, "%f", x[ja[j] * dim + k]); - } - fprintf(fp, "}}"); - } - } - - fprintf(fp, "}],Hue[%f],", /*drand() */ 1.); - - if (A->m < 100) { - for (i = 0; i < A->m; i++) { - if (i > 0) - fprintf(fp, ","); - fprintf(fp, "Text[%d,{", i + 1); - for (k = 0; k < dim; k++) { - if (k > 0) - fprintf(fp, ","); - fprintf(fp, "%f", x[i * dim + k]); - } - fprintf(fp, "}]"); - } - } else if (A->m < 500000) { - fprintf(fp, "Point[{"); - for (i = 0; i < A->m; i++) { - if (i > 0) - fprintf(fp, ","); - fprintf(fp, "{"); - for (k = 0; k < dim; k++) { - if (k > 0) - fprintf(fp, ","); - fprintf(fp, "%f", x[i * dim + k]); - } - fprintf(fp, "}"); - } - fprintf(fp, "}]"); - } else { - fprintf(fp, "{}"); - } - - if (width && dim == 2) { - fprintf(fp, ","); - for (i = 0; i < A->m; i++) { - if (i > 0) - fprintf(fp, ","); - fprintf(fp, - "(*%f,%f*){GrayLevel[.5,.5],Rectangle[{%f,%f},{%f,%f}]}", - width[i * dim], width[i * dim + 1], - x[i * dim] - width[i * dim], - x[i * dim + 1] - width[i * dim + 1], - x[i * dim] + width[i * dim], - x[i * dim + 1] + width[i * dim + 1]); - } - } - - fprintf(fp, "},ImageSize->600]\n"); - -} - -static real update_step(int adaptive_cooling, real step, real Fnorm, - real Fnorm0, real cool) -{ - - if (!adaptive_cooling) { - return cool * step; - } - if (Fnorm >= Fnorm0) { - step = cool * step; - } else if (Fnorm > 0.95 * Fnorm0) { - step = step; - } else { - step = 0.99 * step / cool; - } - return step; -} - - -#define node_degree(i) (ia[(i)+1] - ia[(i)]) - -int check_real_array_size(real ** a, int len, int lenmax) -{ - if (len >= lenmax) { - lenmax = len + MAX((int) 0.2 * len, 10); - *a = RALLOC(lenmax, *a, real); - } - return lenmax; - -} -int check_int_array_size(int **a, int len, int lenmax) -{ - if (len >= lenmax) { - lenmax = len + MAX((int) 0.2 * len, 10); - *a = RALLOC(lenmax, *a, int); - } - return lenmax; - -} - -real get_angle(real * x, int dim, int i, int j) -{ - /* between [0, 2Pi) */ - int k; - real y[2], res; - real eps = 0.00001; - for (k = 0; k < 2; k++) { - y[k] = x[j * dim + k] - x[i * dim + k]; - } - if (ABS(y[0]) <= ABS(y[1]) * eps) { - if (y[1] > 0) - return 0.5 * M_PI; - return 1.5 * M_PI; - } - res = atan(y[1] / y[0]); - if (y[0] > 0) { - if (y[1] < 0) - res = 2 * M_PI + res; - } else if (y[0] < 0) { - res = res + M_PI; - } - return res; -} - -int comp_real(const void *x, const void *y) -{ - real *xx = (real *) x; - real *yy = (real *) y; - - if (*xx > *yy) { - return 1; - } else if (*xx < *yy) { - return -1; - } - return 0; -} -static void sort_real(int n, real * a) -{ - qsort(a, n, sizeof(real), comp_real); -} - - -static void set_leaves(real * x, int dim, real dist, real ang, int i, - int j) -{ - x[dim * j] = cos(ang) * dist + x[dim * i]; - x[dim * j + 1] = sin(ang) * dist + x[dim * i + 1]; -} - -static void beautify_leaves(int dim, SparseMatrix A, real * x) -{ - int m = A->m, i, j, *ia = A->ia, *ja = A->ja, k; - int *checked, p; - real dist; - int nleaves, nleaves_max = 10; - real *angles, maxang, ang1 = 0, ang2 = 0, pad, step; - int *leaves, nangles_max = 10, nangles; - - assert(!SparseMatrix_has_diagonal(A)); - - checked = N_NEW(m, int); - angles = N_GNEW(nangles_max, real); - leaves = N_GNEW(nleaves_max, int); - - for (i = 0; i < m; i++) { - if (ia[i + 1] - ia[i] != 1) - continue; - if (checked[i]) - continue; - p = ja[ia[i]]; - if (!checked[p]) { - checked[p] = TRUE; - dist = 0; - nleaves = 0; - nangles = 0; - for (j = ia[p]; j < ia[p + 1]; j++) { - if (node_degree(ja[j]) == 1) { - checked[ja[j]] = TRUE; - nleaves_max = - check_int_array_size(&leaves, nleaves, - nleaves_max); - dist += distance(x, dim, p, ja[j]); - leaves[nleaves] = ja[j]; - nleaves++; - } else { - nangles_max = - check_real_array_size(&angles, nangles, - nangles_max); - angles[nangles++] = get_angle(x, dim, p, ja[j]); - } - } - assert(nleaves > 0); - dist /= nleaves; - if (nangles > 0) { - sort_real(nangles, angles); - maxang = 0; - for (k = 0; k < nangles - 1; k++) { - if (angles[k + 1] - angles[k] > maxang) { - maxang = angles[k + 1] - angles[k]; - ang1 = angles[k]; - ang2 = angles[k + 1]; - } - } - if (2 * M_PI + angles[0] - angles[nangles - 1] > maxang) { - maxang = 2 * M_PI + angles[0] - angles[nangles - 1]; - ang1 = angles[nangles - 1]; - ang2 = 2 * M_PI + angles[0]; - } - } else { - ang1 = 0; - ang2 = 2 * M_PI; - maxang = 2 * M_PI; - } - pad = MAX(maxang - M_PI * 0.166667 * (nleaves - 1), 0) * 0.5; - ang1 += pad * 0.95; - ang2 -= pad * 0.95; - assert(ang2 >= ang1); - step = 0.; - if (nleaves > 1) - step = (ang2 - ang1) / (nleaves - 1); - for (i = 0; i < nleaves; i++) { - set_leaves(x, dim, dist, ang1, p, leaves[i]); - ang1 += step; - } - } - } - - - free(checked); - free(angles); - free(leaves); -} - -void spring_electrical_spring_embedding(int dim, SparseMatrix A0, - SparseMatrix D, - spring_electrical_control ctrl, - real * node_weights, real * x, - int *flag) -{ - /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. Same as the spring-electrical except we also - introduce force due to spring length - */ - SparseMatrix A = A0; - int m, n; - int i, j, k; - real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = - ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = - ctrl->step, KP; - int *ia = NULL, *ja = NULL; - int *id = NULL, *jd = NULL; - real *d; - real *xold = NULL; - real *f = NULL, dist, F, Fnorm = 0, Fnorm0; - int iter = 0; - int adaptive_cooling = ctrl->adaptive_cooling; - QuadTree qt = NULL; - int USE_QT = FALSE; - int nsuper = 0, nsupermax = 10; - real *center = NULL, *supernode_wgts = NULL, *distances = - NULL, nsuper_avg, counts = 0; - int max_qtree_level = 10; - - if (!A) - return; - m = A->m, n = A->n; - if (n <= 0 || dim <= 0) - return; - - if (n >= ctrl->quadtree_size) { - USE_QT = TRUE; - center = N_GNEW(nsupermax * dim, real); - supernode_wgts = N_GNEW(nsupermax, real); - distances = N_GNEW(nsupermax, real); - } - *flag = 0; - if (m != n) { - *flag = ERROR_NOT_SQUARE_MATRIX; - goto RETURN; - } - assert(A->format == FORMAT_CSR); - A = SparseMatrix_symmetrize(A, TRUE); - ia = A->ia; - ja = A->ja; - id = D->ia; - jd = D->ja; - d = (real *) D->a; - - if (ctrl->random_start) { - srand(ctrl->random_seed); - for (i = 0; i < dim * n; i++) - x[i] = drand(); - } - if (K < 0) { - ctrl->K = K = average_edge_length(A, dim, x); - } - if (C < 0) - ctrl->C = C = 0.2; - if (p >= 0) - ctrl->p = p = -1; - KP = pow(K, 1 - p); - CRK = pow(C, (2. - p) / 3.) / K; - -#ifdef DEBUG_0 - { - FILE *f; - char fname[10000]; - strcpy(fname, "/tmp/graph_layout_0_"); - sprintf(&(fname[strlen(fname)]), "%d", n); - f = fopen(fname, "w"); - export_embedding(f, dim, A, x, NULL); - fclose(f); - } -#endif - - f = N_GNEW(dim, real); - xold = N_GNEW(dim * n, real); - do { - iter++; - xold = memcpy(xold, x, sizeof(real) * dim * n); - Fnorm0 = Fnorm; - Fnorm = 0.; - nsuper_avg = -0; - - if (USE_QT) { - if (ctrl->use_node_weights) { - qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, - x, node_weights); - } else { - qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, - x, NULL); - } - } - - for (i = 0; i < n; i++) { - for (k = 0; k < dim; k++) - f[k] = 0.; - /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */ - - for (j = ia[i]; j < ia[i + 1]; j++) { - if (ja[j] == i) - continue; - dist = distance(x, dim, i, ja[j]); - for (k = 0; k < dim; k++) { - f[k] -= - CRK * (x[i * dim + k] - x[ja[j] * dim + k]) * dist; - } - } - - for (j = id[i]; j < id[i + 1]; j++) { - if (jd[j] == i) - continue; - dist = distance_cropped(x, dim, i, jd[j]); - for (k = 0; k < dim; k++) { - if (dist < d[j]) { - f[k] += - 0.2 * CRK * (x[i * dim + k] - - x[jd[j] * dim + k]) * (dist - - d[j]) * - (dist - d[j]) / dist; - } else { - f[k] -= - 0.2 * CRK * (x[i * dim + k] - - x[jd[j] * dim + k]) * (dist - - d[j]) * - (dist - d[j]) / dist; - } - /* f[k] -= 0.2*CRK*(x[i*dim+k] - x[jd[j]*dim+k])*(dist - d[j]); */ - } - } - - /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */ - if (USE_QT) { - QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim * i]), i, - &nsuper, &nsupermax, ¢er, - &supernode_wgts, &distances, - &counts, flag); - nsuper_avg += nsuper; - if (*flag) - goto RETURN; - for (j = 0; j < nsuper; j++) { - dist = MAX(distances[j], MINDIST); - for (k = 0; k < dim; k++) { - if (p == -1) { - f[k] += - supernode_wgts[j] * KP * (x[i * dim + k] - - center[j * dim + - k]) / - (dist * dist); - } else { - f[k] += - supernode_wgts[j] * KP * (x[i * dim + k] - - center[j * dim + - k]) / - pow(dist, 1. - p); - } - } - } - } else { - if (ctrl->use_node_weights && node_weights) { - for (j = 0; j < n; j++) { - if (j == i) - continue; - dist = distance_cropped(x, dim, i, j); - for (k = 0; k < dim; k++) { - if (p == -1) { - f[k] += - node_weights[j] * KP * - (x[i * dim + k] - - x[j * dim + k]) / (dist * dist); - } else { - f[k] += - node_weights[j] * KP * - (x[i * dim + k] - - x[j * dim + k]) / pow(dist, 1. - p); - } - } - } - } else { - for (j = 0; j < n; j++) { - if (j == i) - continue; - dist = distance_cropped(x, dim, i, j); - for (k = 0; k < dim; k++) { - if (p == -1) { - f[k] += - KP * (x[i * dim + k] - - x[j * dim + k]) / (dist * dist); - } else { - f[k] += - KP * (x[i * dim + k] - - x[j * dim + k]) / pow(dist, - 1. - p); - } - } - } - } - } - - /* normalize force */ - F = 0.; - for (k = 0; k < dim; k++) - F += f[k] * f[k]; - F = sqrt(F); - Fnorm += F; - - if (F > 0) - for (k = 0; k < dim; k++) - f[k] /= F; - - for (k = 0; k < dim; k++) - x[i * dim + k] += step * f[k]; - - } /* done vertex i */ - - if (qt) - QuadTree_delete(qt); - nsuper_avg /= n; -#ifdef DEBUG_PRINT - if (Verbose && 0) { - fprintf(stderr, - "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ", - iter, step, Fnorm, (int) nsuper_avg, A->nz, K); -#ifdef ENERGY - fprintf(stderr, "energy = %f\n", - spring_electrical_energy(dim, A, x, p, CRK, KP)); -#endif - } -#endif - - step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool); - } while (step > tol && iter < maxiter); - -#ifdef DEBUG_PRINT - if (Verbose && 0) - fputs("\n", stderr); -#endif - -#ifdef DEBUG_PRINT_0 - { - FILE *f; - char fname[10000]; - strcpy(fname, "/tmp/graph_layout"); - sprintf(&(fname[strlen(fname)]), "%d", n); - f = fopen(fname, "w"); - export_embedding(f, dim, A, x, NULL); - fclose(f); - } -#endif - - if (ctrl->beautify_leaves) - beautify_leaves(dim, A, x); - - RETURN: - if (xold) - free(xold); - if (A != A0) - SparseMatrix_delete(A); - if (f) - free(f); - if (center) - free(center); - if (supernode_wgts) - free(supernode_wgts); - if (distances) - free(distances); - -} - -int spring_electrical_embedding(int dim, SparseMatrix A0, - spring_electrical_control ctrl, - real * node_weights, real * x) -{ - /* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. */ - SparseMatrix A = A0; - int m, n; - int i, j, k; - real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = - ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = - ctrl->step, KP; - int *ia = NULL, *ja = NULL; - real *xold = NULL; - real *f = NULL, dist, F, Fnorm = 0, Fnorm0; - int iter = 0; - int adaptive_cooling = ctrl->adaptive_cooling; - QuadTree qt = NULL; - int USE_QT = FALSE; - int nsuper = 0, nsupermax = 10; - real *center = NULL, *supernode_wgts = NULL, *distances = - NULL, nsuper_avg, counts = 0, counts_avg = 0; -#ifdef DEBUG - clock_t start, end, start0, start2; - real qtree_cpu = 0, qtree_cpu0 = 0; - real total_cpu = 0; - start0 = clock(); -#endif - int max_qtree_level = ctrl->max_qtree_level; - oned_optimizer *qtree_level_optimizer = NULL; - int flag = 0; - - if (!A) - return flag; - - m = A->m, n = A->n; - if (n <= 0 || dim <= 0) - return flag; - - if (n >= ctrl->quadtree_size) { - USE_QT = TRUE; - qtree_level_optimizer = oned_optimizer_new(max_qtree_level); - center = N_GNEW(nsupermax * dim, real); - supernode_wgts = N_GNEW(nsupermax, real); - distances = N_GNEW(nsupermax, real); - } - if (m != n) { - flag = ERROR_NOT_SQUARE_MATRIX; - goto RETURN; - } - assert(A->format == FORMAT_CSR); - A = SparseMatrix_symmetrize(A, TRUE); - ia = A->ia; - ja = A->ja; - - if (ctrl->random_start) { - srand(ctrl->random_seed); - for (i = 0; i < dim * n; i++) - x[i] = drand(); - } - if (K < 0) { - ctrl->K = K = average_edge_length(A, dim, x); - } - if (C < 0) - ctrl->C = C = 0.2; - if (p >= 0) - ctrl->p = p = -1; - KP = pow(K, 1 - p); - CRK = pow(C, (2. - p) / 3.) / K; - -#ifdef DEBUG_0 - { - FILE *f; - char fname[10000]; - strcpy(fname, "/tmp/graph_layout_0_"); - sprintf(&(fname[strlen(fname)]), "%d", n); - f = fopen(fname, "w"); - export_embedding(f, dim, A, x, NULL); - fclose(f); - } -#endif - - f = N_GNEW(dim, real); - xold = N_GNEW(dim * n, real); - do { - iter++; - xold = memcpy(xold, x, sizeof(real) * dim * n); - Fnorm0 = Fnorm; - Fnorm = 0.; - nsuper_avg = -0; - - if (USE_QT) { - max_qtree_level = oned_optimizer_get(qtree_level_optimizer); - if (ctrl->use_node_weights) { - qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, - x, node_weights); - } else { - qt = QuadTree_new_from_point_list(dim, n, max_qtree_level, - x, NULL); - } - } -#ifdef DEBUG - start2 = clock(); -#endif - - for (i = 0; i < n; i++) { - for (k = 0; k < dim; k++) - f[k] = 0.; - /* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */ - for (j = ia[i]; j < ia[i + 1]; j++) { - if (ja[j] == i) - continue; - dist = distance(x, dim, i, ja[j]); - for (k = 0; k < dim; k++) { - f[k] -= - CRK * (x[i * dim + k] - x[ja[j] * dim + k]) * dist; - } - } - - /* repulsive force K^(1 - p)/||x_i-x_j||^(1 - p) (x_i - x_j) */ - if (USE_QT) { -#ifdef DEBUG - start = clock(); -#endif - QuadTree_get_supernodes(qt, ctrl->bh, &(x[dim * i]), i, - &nsuper, &nsupermax, ¢er, - &supernode_wgts, &distances, - &counts, &flag); -#ifdef DEBUG - end = clock(); - qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC; -#endif - counts_avg += counts; - nsuper_avg += nsuper; - if (flag) - goto RETURN; - for (j = 0; j < nsuper; j++) { - dist = MAX(distances[j], MINDIST); - for (k = 0; k < dim; k++) { - if (p == -1) { - f[k] += - supernode_wgts[j] * KP * (x[i * dim + k] - - center[j * dim + - k]) / - (dist * dist); - } else { - f[k] += - supernode_wgts[j] * KP * (x[i * dim + k] - - center[j * dim + - k]) / - pow(dist, 1. - p); - } - } - } - } else { - if (ctrl->use_node_weights && node_weights) { - for (j = 0; j < n; j++) { - if (j == i) - continue; - dist = distance_cropped(x, dim, i, j); - for (k = 0; k < dim; k++) { - if (p == -1) { - f[k] += - node_weights[j] * KP * - (x[i * dim + k] - - x[j * dim + k]) / (dist * dist); - } else { - f[k] += - node_weights[j] * KP * - (x[i * dim + k] - - x[j * dim + k]) / pow(dist, 1. - p); - } - } - } - } else { - for (j = 0; j < n; j++) { - if (j == i) - continue; - dist = distance_cropped(x, dim, i, j); - for (k = 0; k < dim; k++) { - if (p == -1) { - f[k] += - KP * (x[i * dim + k] - - x[j * dim + k]) / (dist * dist); - } else { - f[k] += - KP * (x[i * dim + k] - - x[j * dim + k]) / pow(dist, - 1. - p); - } - } - } - } - } - - /* normalize force */ - F = 0.; - for (k = 0; k < dim; k++) - F += f[k] * f[k]; - F = sqrt(F); - Fnorm += F; - - if (F > 0) - for (k = 0; k < dim; k++) - f[k] /= F; - - for (k = 0; k < dim; k++) - x[i * dim + k] += step * f[k]; - - } /* done vertex i */ - - if (qt) { - QuadTree_delete(qt); - nsuper_avg /= n; - counts_avg /= n; -#ifdef DEBUG - qtree_cpu0 = qtree_cpu - qtree_cpu0; - if (Verbose && 0) - fprintf(stderr, - "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n", - ((real) (clock() - start2)) / CLOCKS_PER_SEC, - qtree_cpu0, - ((real) (clock() - start2)) / CLOCKS_PER_SEC - - qtree_cpu0); - qtree_cpu0 = qtree_cpu; -#endif - if (Verbose && 0) - fprintf(stderr, - "nsuper_avg=%f, counts_avg = %f 2*nsuper+counts=%f\n", - nsuper_avg, counts_avg, - 2 * nsuper_avg + counts_avg); - oned_optimizer_train(qtree_level_optimizer, - 5 * nsuper_avg + counts_avg); - } -#ifdef DEBUG_PRINT - if (Verbose && 0) { - fprintf(stderr, - "\r iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ", - iter, step, Fnorm, (int) nsuper_avg, A->nz, K); -#ifdef ENERGY - fprintf(stderr, "energy = %f\n", - spring_electrical_energy(dim, A, x, p, CRK, KP)); -#endif - } -#endif - - - step = update_step(adaptive_cooling, step, Fnorm, Fnorm0, cool); - } while (step > tol && iter < maxiter); - -#ifdef DEBUG_PRINT - if (Verbose && 0) - fputs("\n", stderr); -#endif - -#ifdef DEBUG_PRINT_0 - { - FILE *f; - char fname[10000]; - strcpy(fname, "/tmp/graph_layout"); - sprintf(&(fname[strlen(fname)]), "%d", n); - f = fopen(fname, "w"); - export_embedding(f, dim, A, x, NULL); - fclose(f); - } -#endif - - -#ifdef DEBUG_PRINT - if (Verbose) { - if (USE_QT) { - fprintf(stderr, - "iter = %d, step = %f Fnorm = %f qt_level = %d nsuper = %d nz = %d K = %f ", - iter, step, Fnorm, max_qtree_level, (int) nsuper_avg, - A->nz, K); - } else { - fprintf(stderr, - "iter = %d, step = %f Fnorm = %f nsuper = %d nz = %d K = %f ", - iter, step, Fnorm, (int) nsuper_avg, A->nz, K); - } - } -#endif - - if (ctrl->beautify_leaves) - beautify_leaves(dim, A, x); - -#ifdef DEBUG - total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC; - if (Verbose) - fprintf(stderr, "time for supernode = %f, total cpu = %f\n", - qtree_cpu, total_cpu); -#endif - - RETURN: - if (USE_QT) { - oned_optimizer_delete(qtree_level_optimizer); - ctrl->max_qtree_level = max_qtree_level; - } - if (xold) - free(xold); - if (A != A0) - SparseMatrix_delete(A); - if (f) - free(f); - if (center) - free(center); - if (supernode_wgts) - free(supernode_wgts); - if (distances) - free(distances); - - return flag; -} - -void print_matrix(real * x, int n, int dim) -{ - int i, k; - printf("{"); - for (i = 0; i < n; i++) { - if (i != 0) - printf(","); - printf("{"); - for (k = 0; k < dim; k++) { - if (k != 0) - printf(","); - printf("%f", x[i * dim + k]); - } - printf("}"); - } - printf("}\n"); -} - - -static void interpolate(int dim, SparseMatrix A, real * x) -{ - int i, j, k, *ia = A->ia, *ja = A->ja, nz; - real alpha = 0.5, beta, *y; - - y = N_GNEW(dim, real); - for (i = 0; i < A->m; i++) { - for (k = 0; k < dim; k++) - y[k] = 0; - nz = 0; - for (j = ia[i]; j < ia[i + 1]; j++) { - if (ja[j] == i) - continue; - nz++; - for (k = 0; k < dim; k++) { - y[k] += x[ja[j] * dim + k]; - } - } - if (nz > 0) { - beta = (1 - alpha) / nz; - for (k = 0; k < dim; k++) - x[i * dim + k] = alpha * x[i * dim + k] + beta * y[k]; - } - } - - free(y); -} - -static void prolongate(int dim, SparseMatrix A, SparseMatrix P, - SparseMatrix R, real * x, real * y, - int coarsen_scheme_used, real delta) -{ - int nc, *ia, *ja, i, j, k; - SparseMatrix_multiply_dense(P, FALSE, x, FALSE, &y, FALSE, dim); - - /* xu yao rao dong */ - if (coarsen_scheme_used > EDGE_BASED_STA - && coarsen_scheme_used < EDGE_BASED_STO) { - interpolate(dim, A, y); - nc = R->m; - ia = R->ia; - ja = R->ja; - for (i = 0; i < nc; i++) { - for (j = ia[i] + 1; j < ia[i + 1]; j++) { - for (k = 0; k < dim; k++) { - y[ja[j] * dim + k] += delta * (drand() - 0.5); - } - } - } - } -} - - - -int power_law_graph(SparseMatrix A) -{ - int *mask, m, max = 0, i, *ia = A->ia, *ja = A->ja, j, deg; - int res = FALSE; - m = A->m; - mask = N_GNEW(m + 1, int); - - for (i = 0; i < m + 1; i++) { - mask[i] = 0; - } - - for (i = 0; i < m; i++) { - deg = 0; - for (j = ia[i]; j < ia[i + 1]; j++) { - if (i == ja[j]) - continue; - deg++; - } - mask[deg]++; - max = MAX(max, mask[deg]); - } - if (mask[1] > 0.8 * max && mask[1] > 0.3 * m) - res = TRUE; - free(mask); - return res; -} - -void pcp_rotate(int n, int dim, real * x) -{ - int i, k, l; - real y[4], axis[2], center[2], dist, x0, x1; - - assert(dim == 2); - for (i = 0; i < dim * dim; i++) - y[i] = 0; - for (i = 0; i < dim; i++) - center[i] = 0; - for (i = 0; i < n; i++) { - for (k = 0; k < dim; k++) { - center[k] += x[i * dim + k]; - } - } - for (i = 0; i < dim; i++) - center[i] /= n; - for (i = 0; i < n; i++) { - for (k = 0; k < dim; k++) { - x[dim * i + k] = x[dim * i + k] - center[k]; - } - } - - for (i = 0; i < n; i++) { - for (k = 0; k < dim; k++) { - for (l = 0; l < dim; l++) { - y[dim * k + l] += x[i * dim + k] * x[i * dim + l]; - } - } - } - if (y[1] == 0) { - axis[0] = 0; - axis[1] = 1; - } else { - /* Eigensystem[{{x0, x1}, {x1, x3}}] = - {{(x0 + x3 - Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/2, - (x0 + x3 + Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/2}, - {{-(-x0 + x3 + Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/(2*x1), 1}, - {-(-x0 + x3 - Sqrt[x0^2 + 4*x1^2 - 2*x0*x3 + x3^2])/(2*x1), 1}}} - */ - axis[0] = - -(-y[0] + y[3] - - sqrt(y[0] * y[0] + 4 * y[1] * y[1] - 2 * y[0] * y[3] + - y[3] * y[3])) / (2 * y[1]); - axis[1] = 1; - } - dist = sqrt(1 + axis[0] * axis[0]); - axis[0] = axis[0] / dist; - axis[1] = axis[1] / dist; - for (i = 0; i < n; i++) { - x0 = x[dim * i] * axis[0] + x[dim * i + 1] * axis[1]; - x1 = -x[dim * i] * axis[1] + x[dim * i + 1] * axis[0]; - x[dim * i] = x0; - x[dim * i + 1] = x1; - - } - - -} - - -int multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, - spring_electrical_control - ctrl, real * node_weights, - real * label_sizes, real * x) -{ - Multilevel_control mctrl = NULL; - int n, plg, coarsen_scheme_used; - SparseMatrix A = A0; - SparseMatrix P = NULL; - Multilevel grid; - Multilevel grid0; - real *xc = NULL, *xf = NULL; - int flag = 0; -#ifdef DEBUG - clock_t cpu = clock(); -#endif - - if (!A) - return flag; - n = A->n; - if (n <= 0 || dim <= 0) - return flag; - - if (!SparseMatrix_is_symmetric(A, FALSE) - || A->type != MATRIX_TYPE_REAL) { - A = SparseMatrix_get_real_adjacency_matrix_symmetrized(A); - } else { - A = SparseMatrix_remove_diagonal(A); - } - - mctrl = Multilevel_control_new(); - mctrl->maxlevel = ctrl->multilevels; - grid0 = Multilevel_new(A, node_weights, mctrl); - - grid = Multilevel_get_coarsest(grid0); - if (Multilevel_is_finest(grid)) { - xc = x; - } else { - xc = N_GNEW(grid->n * dim, real); - } - - plg = power_law_graph(A); - if (ctrl->p == AUTOP) { - ctrl->p = -1; - if (plg) - ctrl->p = -1.8; - } - - do { -#ifdef DEBUG_PRINT - if (Verbose) { - print_padding(grid->level); - if (Multilevel_is_coarsest(grid)) { - fprintf(stderr, "coarsest level -- %d, n = %d\n", - grid->level, grid->n); - } else { - fprintf(stderr, "level -- %d, n = %d\n", grid->level, - grid->n); - } - } -#endif - flag = spring_electrical_embedding(dim, grid->A, ctrl, grid->node_weights, - xc); - if (Multilevel_is_finest(grid)) - break; - if (flag) { - free(xc); - goto RETURN; - } - P = grid->P; - coarsen_scheme_used = grid->coarsen_scheme_used; - grid = grid->prev; - if (Multilevel_is_finest(grid)) { - xf = x; - } else { - xf = N_GNEW(grid->n * dim, real); - } - prolongate(dim, grid->A, P, grid->R, xc, xf, coarsen_scheme_used, - (ctrl->K) * 0.001); - free(xc); - xc = xf; - ctrl->random_start = FALSE; - ctrl->K = ctrl->K * 0.75; - ctrl->adaptive_cooling = FALSE; - if (grid->next->coarsen_scheme_used > VERTEX_BASED_STA && - grid->next->coarsen_scheme_used < VERTEX_BASED_STO) { - ctrl->step = 1; - } else { - ctrl->step = .1; - } - } while (grid); - -#ifdef DEBUG - if (Verbose) - fprintf(stderr, "layout time %f\n", - ((real) (clock() - cpu)) / CLOCKS_PER_SEC); - cpu = clock(); -#endif - - flag = post_process_smoothing(dim, A, ctrl, node_weights, x); - - if (dim == 2) { - pcp_rotate(n, dim, x); - } - - if (Verbose) - fprintf(stderr, "ctrl->overlap=%d\n", ctrl->overlap); - - flag = remove_overlap(dim, A, x, label_sizes, ctrl->overlap); - - RETURN: - if (A != A0) - SparseMatrix_delete(A); - Multilevel_control_delete(mctrl); - Multilevel_delete(grid0); - return flag; -} diff --git a/lib/sfdpgen/spring_electrical.h b/lib/sfdpgen/spring_electrical.h deleted file mode 100644 index 646a55f33..000000000 --- a/lib/sfdpgen/spring_electrical.h +++ /dev/null @@ -1,86 +0,0 @@ -/* 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 SPRING_ELECTRICAL_H -#define SPRING_ELECTRICAL_H - -#include - -enum { ERROR_NOT_SQUARE_MATRIX = -100 }; - -/* a flag to indicate that p should be set auto */ -#define AUTOP -1.0001234 - -#define MINDIST 1.e-15 - -typedef enum { - SMOOTHING_NONE, - SMOOTHING_STRESS_MAJORIZATION_GRAPH_DIST, - SMOOTHING_STRESS_MAJORIZATION_AVG_DIST, - SMOOTHING_STRESS_MAJORIZATION_POWER_DIST, - SMOOTHING_SPRING, - SMOOTHING_TRIANGLE, - SMOOTHING_RNG -} smooth_t; - -struct spring_electrical_control_s { - real p; /*a negative number default to -1. repulsive force = dist^p */ - int random_start; /* whether to apply SE from a random layout, or from exisiting layout */ - real K; /* the natural distance. If K < 0, K will be set to the average distance of an edge */ - real C; /* another parameter. f_a(i,j) = C*dist(i,j)^2/K * d_ij, f_r(i,j) = K^(3-p)/dist(i,j)^(-p). By default C = 0.2. */ - int multilevels; /* if <=1, single level */ - int quadtree_size; /* cut off size above which quadtree approximation is used */ - int max_qtree_level; /* max level of quadtree */ - real bh; /* Barnes-Hutt constant, if width(snode)/dist[i,snode] < bh, treat snode as a supernode. default 0.2 */ - real tol; /* minimum different between two subsequence config before terminating. ||x-xold|| < tol */ - int maxiter; - real cool; /* default 0.9 */ - real step; /* initial step size */ - int adaptive_cooling; - int random_seed; - int beautify_leaves; - int use_node_weights; - int smoothing; - int overlap; -}; - -typedef struct spring_electrical_control_s* spring_electrical_control; - -spring_electrical_control spring_electrical_control_new(); -int spring_electrical_embedding(int dim, SparseMatrix A0, - spring_electrical_control ctrl, - real * node_weights, real * x); - -int multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, - spring_electrical_control - ctrl0, real * node_weights, - real * label_sizes, real * x); - -void export_embedding(FILE * fp, int dim, SparseMatrix A, real * x, - real * width); -void spring_electrical_control_delete(spring_electrical_control); -void print_matrix(real * x, int n, int dim); - -real distance(real * x, int dim, int i, int j); -real distance_cropped(real * x, int dim, int i, int j); -real average_edge_length(SparseMatrix A, int dim, real * coord); - -void spring_electrical_spring_embedding(int dim, SparseMatrix A, - SparseMatrix D, - spring_electrical_control ctrl, - real * node_weights, real * x, - int *flag); - -#endif -- 2.40.0