+++ /dev/null
-/* 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;
-}
+++ /dev/null
-/* 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
+++ /dev/null
-# $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
+++ /dev/null
-/* 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;
-}
+++ /dev/null
-/* 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
+++ /dev/null
-/* 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 <math.h>
-#include <assert.h>
-
-typedef struct {
- real node_weight;
- real *coord;
- real id;
-} node_data;
-
-real point_distance(real * p1, real * p2, int dim)
-{
- int i;
- real dist;
- dist = 0;
- for (i = 0; i < dim; i++)
- dist += (p1[i] - p2[i]) * (p1[i] - p2[i]);
- return sqrt(dist);
-}
-
-static node_data *node_data_new(int dim, real weight, real * coord, int id)
-{
- node_data *nd;
- int i;
- nd = 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");
- }
-}
+++ /dev/null
-/* vim:set shiftwidth=4 ts=8: */
-
-/**********************************************************
-* This software is part of the graphviz package *
-* http://www.graphviz.org/ *
-* *
-* Copyright (c) 1994-2004 AT&T Corp. *
-* and is licensed under the *
-* Common Public License, Version 1.0 *
-* by AT&T Corp. *
-* *
-* Information and Software Systems Research *
-* AT&T Research, Florham Park NJ *
-**********************************************************/
-
-#ifndef QUAD_TREE_H
-#define QUAD_TREE_H
-
-#include "LinkedList.h"
-#include "sfdpinternal.h"
-#include <stdio.h>
-
-typedef struct QuadTree_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
+++ /dev/null
-/* 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);
-
-}
+++ /dev/null
-/* 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 <stdio.h>
-#include <sfdpinternal.h>
-
-#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
+++ /dev/null
-/* 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;
-
-}
+++ /dev/null
-/* 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
+++ /dev/null
-#include <stdio.h>
-#include <stdlib.h>
-#include <getopt.h>
-#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;
-}
-
+++ /dev/null
-/* 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 <time.h>
-
-#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;
-}
+++ /dev/null
-
-/* 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
+++ /dev/null
-/* vim:set shiftwidth=4 ts=8: */
-
-/**********************************************************
-* This software is part of the graphviz package *
-* http://www.graphviz.org/ *
-* *
-* Copyright (c) 1994-2004 AT&T Corp. *
-* and is licensed under the *
-* Common Public License, Version 1.0 *
-* by AT&T Corp. *
-* *
-* Information and Software Systems Research *
-* AT&T Research, Florham Park NJ *
-**********************************************************/
-
-
-#include "post_process.h"
-#include "spring_electrical.h"
-#include "sparse_solve.h"
-#include "call_tri.h"
-#include "types.h"
-#include "memory.h"
-#include "globals.h"
-#include <time.h>
-#include <string.h>
-#include <math.h>
-
-#define node_degree(i) (ia[(i)+1] - ia[(i)])
-
-SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, real * x)
-{
- /* find the ideal distance between edges, either 1, or |N[i] \Union N[j]| - |N[i] \Intersection N[j]|
- */
- SparseMatrix D;
- int *ia, *ja, i, j, k, l, nz;
- real *d;
- int *mask = NULL;
- real len, di, sum, sumd;
-
- assert(SparseMatrix_is_symmetric(A, FALSE));
-
- D = SparseMatrix_copy(A);
- ia = D->ia;
- ja = D->ja;
- if (D->type != MATRIX_TYPE_REAL) {
- free(D->a);
- D->type = MATRIX_TYPE_REAL;
- D->a = 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;
-}
+++ /dev/null
-/* vim:set shiftwidth=4 ts=8: */
-
-/**********************************************************
-* This software is part of the graphviz package *
-* http://www.graphviz.org/ *
-* *
-* Copyright (c) 1994-2004 AT&T Corp. *
-* and is licensed under the *
-* Common Public License, Version 1.0 *
-* by AT&T Corp. *
-* *
-* Information and Software Systems Research *
-* AT&T Research, Florham Park NJ *
-**********************************************************/
-
-#ifndef POST_PROCESS_H
-#define POST_PROCESS_H
-
-#include <spring_electrical.h>
-
-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
+++ /dev/null
-/* 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 <stdlib.h>
-#include <stdio.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. */
-
-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);
- }
-}
+++ /dev/null
-/* 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
+++ /dev/null
-/* 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
+++ /dev/null
-/* 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 <limits.h>
-#else
-#ifdef HAVE_VALUES_H
-#include <values.h>
-#endif
-#endif
-#include <sfdp.h>
-#include <neatoprocs.h>
-#include <adjust.h>
-#include <pack.h>
-#include <assert.h>
-#include <ctype.h>
-#include <spring_electrical.h>
-#include <overlap.h>
-
-#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);
-}
+++ /dev/null
-/* 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 <sfdp.h>
-
-#ifdef DEBUG
-extern double _statistics[10];
-#endif
-
-typedef double real;
-
-#endif
+++ /dev/null
-/* 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;
-}
-
+++ /dev/null
-/* 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
+++ /dev/null
-/* 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 <string.h>
-#include <time.h>
-
-#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;
-}
+++ /dev/null
-/* 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 <SparseMatrix.h>
-
-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
projects / graphviz / commitdiff