From d1d1ae30865c14bb7e6a75f8772543c10419e48b Mon Sep 17 00:00:00 2001 From: erg Date: Mon, 28 Feb 2011 16:08:27 +0000 Subject: [PATCH] Add new libraries and features needed for additions to sfdp and gvmap --- lib/sparse/BinaryHeap.c | 291 ++++++++++++ lib/sparse/BinaryHeap.h | 77 +++ lib/sparse/IntStack.c | 95 ++++ lib/sparse/IntStack.h | 39 ++ lib/sparse/Makefile.am | 4 +- lib/sparse/SparseMatrix.c | 963 ++++++++++++++++++++++++++++++++++---- lib/sparse/SparseMatrix.h | 38 +- lib/sparse/general.c | 300 ++++++++++++ lib/sparse/general.h | 141 ++++++ 9 files changed, 1836 insertions(+), 112 deletions(-) create mode 100644 lib/sparse/BinaryHeap.c create mode 100644 lib/sparse/BinaryHeap.h create mode 100644 lib/sparse/IntStack.c create mode 100644 lib/sparse/IntStack.h create mode 100644 lib/sparse/general.c create mode 100644 lib/sparse/general.h diff --git a/lib/sparse/BinaryHeap.c b/lib/sparse/BinaryHeap.c new file mode 100644 index 000000000..0321efda6 --- /dev/null +++ b/lib/sparse/BinaryHeap.c @@ -0,0 +1,291 @@ +/* $Id$Revision: */ +/* vim:set shiftwidth=4 ts=8: */ + +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#include "BinaryHeap.h" + +BinaryHeap BinaryHeap_new(int (*cmp)(void*item1, void*item2)){ + BinaryHeap h; + int max_len = 1<<8, i; + + h = MALLOC(sizeof(struct BinaryHeap_struct)); + h->max_len = max_len; + h->len = 0; + h->heap = MALLOC(sizeof(void*)*max_len); + h->id_to_pos = MALLOC(sizeof(int)*max_len); + for (i = 0; i < max_len; i++) (h->id_to_pos)[i] = -1; + + h->pos_to_id = MALLOC(sizeof(int)*max_len); + h->id_stack = IntStack_new(); + h->cmp = cmp; + return h; +} + + +void BinaryHeap_delete(BinaryHeap h, void (*del)(void* item)){ + int i; + if (!h) return; + FREE(h->id_to_pos); + FREE(h->pos_to_id); + IntStack_delete(h->id_stack); + if (del) for (i = 0; i < h->len; i++) del((h->heap)[i]); + FREE(h->heap); + FREE(h); +} + +static BinaryHeap BinaryHeap_realloc(BinaryHeap h){ + int max_len0 = h->max_len, max_len = h->max_len, i; + max_len = max_len + MAX(0.2*max_len, 10); + h->max_len = max_len; + + h->heap = REALLOC(h->heap, sizeof(void*)*max_len); + if (!(h->heap)) return NULL; + + h->id_to_pos = REALLOC(h->id_to_pos, sizeof(int)*max_len); + if (!(h->id_to_pos)) return NULL; + + h->pos_to_id = REALLOC(h->pos_to_id, sizeof(int)*max_len); + if (!(h->pos_to_id)) return NULL; + + for (i = max_len0; i < max_len; i++) (h->id_to_pos)[i] = -1; + return h; + +} + +#define ParentPos(pos) (pos - 1)/2 +#define ChildrenPos1(pos) (2*(pos)+1) +#define ChildrenPos2(pos) (2*(pos)+2) + +static void swap(BinaryHeap h, int parentPos, int nodePos){ + int parentID, nodeID; + void *tmp; + void **heap = h->heap; + int *id_to_pos = h->id_to_pos, *pos_to_id = h->pos_to_id; + + assert(parentPos < h->len); + assert(nodePos < h->len); + + parentID = pos_to_id[parentPos]; + nodeID = pos_to_id[nodePos]; + + tmp = heap[parentPos]; + heap[parentPos] = heap[nodePos]; + heap[nodePos] = tmp; + + pos_to_id[parentPos] = nodeID; + id_to_pos[nodeID] = parentPos; + + pos_to_id[nodePos] = parentID; + id_to_pos[parentID] = nodePos; + +} + +static int siftUp(BinaryHeap h, int nodePos){ + int parentPos; + + void **heap = h->heap; + + + if (nodePos != 0) { + parentPos = ParentPos(nodePos); + + if ((h->cmp)(heap[parentPos], heap[nodePos]) == 1) {/* if smaller than parent, swap */ + swap(h, parentPos, nodePos); + nodePos = siftUp(h, parentPos); + } + } + return nodePos; +} + +static int siftDown(BinaryHeap h, int nodePos){ + int childPos, childPos1, childPos2; + + void **heap = h->heap; + + + childPos1 = ChildrenPos1(nodePos); + childPos2 = ChildrenPos2(nodePos); + if (childPos1 > h->len - 1) return nodePos;/* no child */ + if (childPos1 == h->len - 1) { + childPos = childPos1;/* one child */ + } else {/* two child */ + if ((h->cmp)(heap[childPos1], heap[childPos2]) == 1){ /* pick the smaller of the two child */ + childPos = childPos2; + } else { + childPos = childPos1; + } + } + + if ((h->cmp)(heap[nodePos], heap[childPos]) == 1) { + /* if larger than child, swap */ + swap(h, nodePos, childPos); + nodePos = siftDown(h, childPos); + } + + return nodePos; +} + +int BinaryHeap_insert(BinaryHeap h, void *item){ + int len = h->len, id = len, flag, pos; + + /* insert an item, and return its ID. This ID can be used later to extract the item */ + if (len > h->max_len - 1) { + if (BinaryHeap_realloc(h) == NULL) return BinaryHeap_error_malloc; + } + + /* check if we have IDs in the stack to reuse. If no, then assign the last pos as the ID */ + id = IntStack_pop(h->id_stack, &flag); + if (flag) id = len; + + h->heap[len] = item; + h->id_to_pos[id] = len; + h->pos_to_id[len] = id; + + (h->len)++; + + pos = siftUp(h, len); + assert(h->id_to_pos[id] == pos); + assert(h->pos_to_id[pos] == id); + + + return id; +} + + +void* BinaryHeap_get_min(BinaryHeap h){ + /* return the min item */ + return h->heap[0]; +} + +void* BinaryHeap_extract_min(BinaryHeap h){ + /* return and remove the min item */ + if (h->len == 0) return NULL; + return BinaryHeap_extract_item(h, (h->pos_to_id)[0]); +} + +void* BinaryHeap_extract_item(BinaryHeap h, int id){ + /* extract an item with ID out and delete it */ + void *item; + int *id_to_pos = h->id_to_pos; + int pos; + + if (id >= h->max_len) return NULL; + + pos = id_to_pos[id]; + + if (pos < 0) return NULL; + + assert(pos < h->len); + + item = (h->heap)[pos]; + + IntStack_push(h->id_stack, id); + + if (pos < h->len - 1){/* move the last item to occupy the position of extracted item */ + swap(h, pos, h->len - 1); + (h->len)--; + pos = siftUp(h, pos); + pos = siftDown(h, pos); + } else { + (h->len)--; + } + + (h->id_to_pos)[id] = -1; + + return item; + +} + +int BinaryHeap_reset(BinaryHeap h, int id, void *item){ + /* reset value of an item with specified id */ + int pos; + + if (id >= h->max_len) return -1; + pos = h->id_to_pos[id]; + if (pos < 0) return -1; + + h->heap[pos] = item; + + pos = siftUp(h, pos); + + pos = siftDown(h, pos); + + return pos; + +} +void* BinaryHeap_get_item(BinaryHeap h, int id){ + /* get an item with ID out, without deleting */ + int pos; + + if (id >= h->max_len) return NULL; + + pos = h->id_to_pos[id]; + + if (pos < 0) return NULL; + return (h->heap)[pos]; +} + +void BinaryHeap_sanity_check(BinaryHeap h){ + int i, key_spare, parentPos, *id_to_pos = h->id_to_pos, *pos_to_id = h->pos_to_id; + void **heap = h->heap; + int *mask; + + /* check that this is a binary heap: children is smaller than parent */ + for (i = 1; i < h->len; i++){ + parentPos = ParentPos(i); + assert((h->cmp)(heap[i], heap[parentPos]) >= 0); + } + + mask = MALLOC(sizeof(int)*(h->len + IntStack_get_length(h->id_stack))); + for (i = 0; i < h->len + IntStack_get_length(h->id_stack); i++) mask[i] = -1; + + /* check that spare keys has negative id_to_pos mapping */ + for (i = 0; i <= h->id_stack->last; i++){ + key_spare = h->id_stack->stack[i]; + assert(h->id_to_pos[key_spare] < 0); + mask[key_spare] = 1;/* mask spare ID */ + } + + /* check that + pos_to_id[id_to_pos[i]] = i, for i not in the id_stack & i < length(id_stack)+len + id_to_pos[pos_to_id[i]] = i, 0 <= i < len + */ + for (i = 1; i < h->len; i++){ + assert(mask[pos_to_id[i]] == -1);/* that id is in use so can't be spare */ + mask[pos_to_id[i]] = 1; + assert(id_to_pos[pos_to_id[i]] == i); + } + + /* all IDs, spare or in use, are ccounted for and give a contiguous set */ + for (i = 0; i < h->len + IntStack_get_length(h->id_stack); i++) assert(mask[i] =- 1); + + FREE(mask); +} +void BinaryHeap_print(BinaryHeap h, void (*pnt)(void*)){ + int i, k = 2; + + for (i = 0; i < h->len; i++){ + pnt(h->heap[i]); + fprintf(stderr, "(%d) ",(h->pos_to_id)[i]); + if (i == k-2) { + fprintf(stderr, "\n"); + k *= 2; + } + } + fprintf(stderr, "\nSpare keys ="); + for (i = 0; i <= h->id_stack->last; i++){ + fprintf(stderr,"%d(%d) ",h->id_stack->stack[i], h->id_to_pos[h->id_stack->stack[i]]); + } + fprintf(stderr, "\n"); +} + + diff --git a/lib/sparse/BinaryHeap.h b/lib/sparse/BinaryHeap.h new file mode 100644 index 000000000..946442511 --- /dev/null +++ b/lib/sparse/BinaryHeap.h @@ -0,0 +1,77 @@ +/* $Id$Revision: */ +/* vim:set shiftwidth=4 ts=8: */ + +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifndef BinaryHeap_H +#define BinaryHeap_H + +#include "general.h" +#include "IntStack.h" + +/* binary heap code. + Caution: items inserted should be kept untouched, e.g., the value of the item should be kepted unchanged while the heap is still in use! + To change the valud of am item, use BinaryHeap_reset +*/ + +struct BinaryHeap_struct { + int max_len;/* storage allocated for the heap */ + int len;/*number of elements in the heap so far. <= maxlen */ + void **heap; + int *id_to_pos;/* this array store the position of an item with ID. For ID that was not in use, + i.e., no in pos_to_id[1:len], + id_to_pos[id] = -1 + */ + int *pos_to_id;/* this array stores the ID of item at position pos. + pos_to_id[id_to_pos[i]] = i, for i not in the id_stack & i < length(id_stack)+len + id_to_pos[pos_to_id[i]] = i, 0 <= i < len + */ + IntStack id_stack;/* a stack that store IDs that is no longer used, to + be assigned to newly inserted elements. + For sanity check, the union of items in + the id_stack, and that is pos_to_id[1:len], + should form a set of contiguous numbers + {1, 2, ..., len, ...} + */ + int (*cmp)(void*item1, void*item2);/* comparison function. Return 1,0,-1 + if item1 >, =, < item2 */ +}; + +enum {BinaryHeap_error_malloc = -10}; + +typedef struct BinaryHeap_struct* BinaryHeap; + +BinaryHeap BinaryHeap_new(int (*cmp)(void*item1, void*item2)); + + +void BinaryHeap_delete(BinaryHeap h, void (*del)(void* item));/* delete not just the heap data structure, but also the data items + through a user supplied del function. */ + +int BinaryHeap_insert(BinaryHeap h, void *item);/* insert an item, and return its ID. + This ID can be used later to extract the item. ID + are always nonnegative. If the return value is + negative, it is an error message */ + +void* BinaryHeap_get_min(BinaryHeap h);/* return the min item */ + +void* BinaryHeap_extract_min(BinaryHeap h);/* return and remove the min item */ + +void* BinaryHeap_extract_item(BinaryHeap h, int id);/* extract an item with ID out and delete it */ + +void* BinaryHeap_get_item(BinaryHeap h, int id);/* get an item with ID out, without deleting */ + +int BinaryHeap_reset(BinaryHeap h, int id, void *item);/* reset value of an item with specified id. Return new pos. If ID is invalue, return -1 */ + +void BinaryHeap_print(BinaryHeap h, void (*pnt)(void*)); + +void BinaryHeap_sanity_check(BinaryHeap h); + +#endif diff --git a/lib/sparse/IntStack.c b/lib/sparse/IntStack.c new file mode 100644 index 000000000..8961174a9 --- /dev/null +++ b/lib/sparse/IntStack.c @@ -0,0 +1,95 @@ +/* $Id$Revision: */ +/* vim:set shiftwidth=4 ts=8: */ + +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#include "general.h" +#include "IntStack.h" + +IntStack IntStack_new(void){ + IntStack s; + int max_len = 1<<5; + + s = MALLOC(sizeof(struct IntStack_struct)); + s->max_len = max_len; + s->last = -1; + s->stack = MALLOC(sizeof(int)*max_len); + return s; +} + +void IntStack_delete(IntStack s){ + if (s){ + FREE(s->stack); + FREE(s); + } +} + +static IntStack IntStack_realloc(IntStack s){ + int max_len = s->max_len; + + max_len = max_len + MAX(10,0.2*max_len); + s->max_len = max_len; + s->stack = REALLOC(s->stack, sizeof(int)*max_len); + if (!s->stack) return NULL; + return s; +} + +int IntStack_push(IntStack s, int i){ + /* add an item and return the pos. Return negative value of malloc failed */ + if (s->last >= s->max_len - 1){ + if (!(IntStack_realloc(s))) return -1; + } + s->stack[++(s->last)] = i; + return s->last; +} +int IntStack_pop(IntStack s, int *flag){ + /* remove the last item. If none exist, return -1 */ + *flag = 0; + if (s->last < 0){ + *flag = -1; + return -1; + } + return s->stack[(s->last)--]; +} +void IntStack_print(IntStack s){ + /* remove the last item. If none exist, return -1 */ + int i; + for (i = 0; i <= s->last; i++) fprintf(stderr,"%d,",s->stack[i]); + fprintf(stderr,"\n"); +} + +/* +main(){ + + IntStack s; + int i, len = 1, pos, flag; + + for (;;){ + s = IntStack_new(); + fprintf(stderr,"=============== stack with %d elements ============\n",len); + for (i = 0; i < len; i++){ + pos = IntStack_push(s, i); + if (pos < 0){ + fprintf(stderr," fail to push element %d, quit\n", i); + exit(1); + } + } + for (i = 0; i < len+1; i++){ + IntStack_pop(s, &flag); + if (flag) { + fprintf(stderr, "no more element\n"); + } + } + IntStack_delete(s); + len *= 2; + } +} +*/ diff --git a/lib/sparse/IntStack.h b/lib/sparse/IntStack.h new file mode 100644 index 000000000..a1bfbe8f4 --- /dev/null +++ b/lib/sparse/IntStack.h @@ -0,0 +1,39 @@ +/* $Id$Revision: */ +/* vim:set shiftwidth=4 ts=8: */ + +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifndef IntStack_H +#define IntStack_H + +/* last in first out integer stack */ +struct IntStack_struct{ + int last;/* position of the last element, If empty, last = -1 */ + int max_len; + int *stack; +}; + +typedef struct IntStack_struct* IntStack; + +IntStack IntStack_new(void); + +void IntStack_delete(IntStack s); + +#define IntStack_get_length(s) (1+(s)->last) + +int IntStack_push(IntStack s, int i);/* add an item and return the pos (>=0). + Return negative value of malloc failed */ + +int IntStack_pop(IntStack s, int *flag);/* remove the last item. If none exist, flag = -1, and return -1. */ + +void IntStack_print(IntStack s); + +#endif diff --git a/lib/sparse/Makefile.am b/lib/sparse/Makefile.am index 680e7abe6..f84c189b2 100644 --- a/lib/sparse/Makefile.am +++ b/lib/sparse/Makefile.am @@ -5,9 +5,9 @@ AM_CPPFLAGS = \ -I$(top_srcdir) \ -I$(top_srcdir)/lib/common -noinst_HEADERS = SparseMatrix.h +noinst_HEADERS = SparseMatrix.h general.h BinaryHeap.h IntStack.h noinst_LTLIBRARIES = libsparse_C.la -libsparse_C_la_SOURCES = SparseMatrix.c +libsparse_C_la_SOURCES = SparseMatrix.c general.c BinaryHeap.c IntStack.c EXTRA_DIST = Makefile.old gvsparse.vcproj diff --git a/lib/sparse/SparseMatrix.c b/lib/sparse/SparseMatrix.c index 9a001a090..40011f77f 100644 --- a/lib/sparse/SparseMatrix.c +++ b/lib/sparse/SparseMatrix.c @@ -1,3 +1,4 @@ +/* $Id$Revision: */ /* vim:set shiftwidth=4 ts=8: */ /************************************************************************* @@ -9,22 +10,16 @@ * * Contributors: See CVS logs. Details at http://www.graphviz.org/ *************************************************************************/ + #include #include #include #include - -#include "SparseMatrix.h" - #include "logic.h" #include "memory.h" #include "arith.h" - -#define MALLOC gmalloc -#define REALLOC grealloc -#define FREE free -#define MEMCPY memcpy - +#include "SparseMatrix.h" +#include "BinaryHeap.h" static size_t size_of_matrix_type(int type){ int size = 0; @@ -52,6 +47,14 @@ static size_t size_of_matrix_type(int type){ return size; } +SparseMatrix SparseMatrix_sort(SparseMatrix A){ + SparseMatrix B; + B = SparseMatrix_transpose(A); + SparseMatrix_delete(A); + A = SparseMatrix_transpose(B); + SparseMatrix_delete(B); + return A; +} SparseMatrix SparseMatrix_make_undirected(SparseMatrix A){ /* make it strictly low diag only, and set flag to undirected */ SparseMatrix B; @@ -152,6 +155,21 @@ SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, int pattern_symmetric_only) return A; } +SparseMatrix SparseMatrix_symmetrize_nodiag(SparseMatrix A, int pattern_symmetric_only){ + SparseMatrix B; + if (SparseMatrix_is_symmetric(A, pattern_symmetric_only)) { + B = SparseMatrix_copy(A); + return SparseMatrix_remove_diagonal(B); + } + 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 SparseMatrix_remove_diagonal(A); +} + int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only){ /* assume no repeated entries! */ SparseMatrix B; @@ -177,7 +195,7 @@ int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only){ jb = B->ja; m = A->m; - mask = N_GNEW(m,int); + mask = MALLOC(sizeof(int)*m); for (i = 0; i < m; i++) mask[i] = -1; type = A->type; @@ -235,6 +253,7 @@ int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only){ if (bi[j] != ai[mask[jb[j]]]) goto RETURN; } } + res = TRUE; break; } case MATRIX_TYPE_PATTERN: @@ -246,6 +265,7 @@ int SparseMatrix_is_symmetric(SparseMatrix A, int test_pattern_symmetry_only){ if (mask[jb[j]] < ia[i]) goto RETURN; } } + res = TRUE; break; case MATRIX_TYPE_UNKNOWN: goto RETURN; @@ -272,7 +292,7 @@ static SparseMatrix SparseMatrix_init(int m, int n, int type, int format){ SparseMatrix A; - A = N_GNEW(1,struct SparseMatrix_struct); + A = MALLOC(sizeof(struct SparseMatrix_struct)); A->m = m; A->n = n; A->nz = 0; @@ -285,7 +305,7 @@ static SparseMatrix SparseMatrix_init(int m, int n, int type, int format){ case FORMAT_CSC: case FORMAT_CSR: default: - A->ia = N_GNEW((m+1),int); + A->ia = MALLOC(sizeof(int)*(m+1)); } A->ja = NULL; A->a = NULL; @@ -305,14 +325,14 @@ static SparseMatrix SparseMatrix_alloc(SparseMatrix A, int nz){ A->a = NULL; switch (format){ case FORMAT_COORD: - A->ia = N_GNEW(nz,int); - A->ja = N_GNEW(nz,int); + A->ia = MALLOC(sizeof(int)*nz); + A->ja = MALLOC(sizeof(int)*nz); A->a = MALLOC(size_of_matrix_type(type)*nz); break; case FORMAT_CSR: case FORMAT_CSC: default: - A->ja = N_GNEW(nz,int); + A->ja = MALLOC(sizeof(int)*nz); if (size_of_matrix_type(type) && nz > 0) A->a = MALLOC(size_of_matrix_type(type)*nz); break; } @@ -324,8 +344,8 @@ 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); + A->ia = REALLOC(A->ia, sizeof(int)*nz); + A->ja = REALLOC(A->ja, sizeof(int)*nz); if (size_of_matrix_type(type) > 0) { if (A->a){ A->a = REALLOC(A->a, size_of_matrix_type(type)*nz); @@ -337,7 +357,7 @@ static SparseMatrix SparseMatrix_realloc(SparseMatrix A, int nz){ case FORMAT_CSR: case FORMAT_CSC: default: - A->ja = RALLOC(nz,A->ja,int); + A->ja = REALLOC(A->ja, sizeof(int)*nz); if (size_of_matrix_type(type) > 0) { if (A->a){ A->a = REALLOC(A->a, size_of_matrix_type(type)*nz); @@ -377,6 +397,7 @@ void SparseMatrix_print_csr(char *c, SparseMatrix A){ int *ai; int i, j, m = A->m; + assert (A->format == FORMAT_CSR); printf("%s\n SparseArray[{",c); ia = A->ia; ja = A->ja; @@ -437,6 +458,7 @@ void SparseMatrix_print_coord(char *c, SparseMatrix A){ int *ai; int i, m = A->m; + assert (A->format == FORMAT_COORD); printf("%s\n SparseArray[{",c); ia = A->ia; ja = A->ja; @@ -577,7 +599,6 @@ static void SparseMatrix_export_csr(FILE *f, SparseMatrix A){ void SparseMatrix_export_binary(char *name, SparseMatrix A, int *flag){ FILE *f; - int rc; *flag = 0; f = fopen(name, "wb"); @@ -585,62 +606,52 @@ void SparseMatrix_export_binary(char *name, SparseMatrix A, int *flag){ *flag = 1; return; } - rc = fwrite(&(A->m), sizeof(int), 1, f); - rc += fwrite(&(A->n), sizeof(int), 1, f); - rc += fwrite(&(A->nz), sizeof(int), 1, f); - rc += fwrite(&(A->nzmax), sizeof(int), 1, f); - rc += fwrite(&(A->type), sizeof(int), 1, f); - rc += fwrite(&(A->format), sizeof(int), 1, f); - rc += fwrite(&(A->property), sizeof(int), 1, f); - if (rc != 7) 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){ - rc = fwrite(A->ia, sizeof(int), A->nz, f); - if (rc != A->nz) return; + fwrite(A->ia, sizeof(int), A->nz, f); } else { - rc = fwrite(A->ia, sizeof(int), A->m + 1, f); - if (rc != A->m + 1) return; - } - rc = fwrite(A->ja, sizeof(int), A->nz, f); - if (rc != A->nz) return; - if (size_of_matrix_type(A->type) > 0) { - rc = fwrite(A->a, size_of_matrix_type(A->type), A->nz, f); - if (rc != A->nz) return; + 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, rc; + int m, n, nz, nzmax, type, format, property, iread; FILE *f; f = fopen(name, "rb"); if (!f) return NULL; - rc = fread(&m, sizeof(int), 1, f); - rc += fread(&n, sizeof(int), 1, f); - rc += fread(&nz, sizeof(int), 1, f); - rc += fread(&nzmax, sizeof(int), 1, f); - rc += fread(&type, sizeof(int), 1, f); - rc += fread(&format, sizeof(int), 1, f); - rc += fread(&property, sizeof(int), 1, f); - if (rc != 7) return NULL; + iread = fread(&m, sizeof(int), 1, f); + iread = fread(&n, sizeof(int), 1, f); + iread = fread(&nz, sizeof(int), 1, f); + iread = fread(&nzmax, sizeof(int), 1, f); + iread = fread(&type, sizeof(int), 1, f); + iread = fread(&format, sizeof(int), 1, f); + iread = fread(&property, sizeof(int), 1, f); + A = SparseMatrix_new(m, n, nz, type, format); A->nz = nz; A->property = property; - + if (format == FORMAT_COORD){ - rc = fread(A->ia, sizeof(int), A->nz, f); - if (rc != A->nz) return NULL; + iread = fread(A->ia, sizeof(int), A->nz, f); } else { - rc = fread(A->ia, sizeof(int), A->m + 1, f); - if (rc != A->m + 1) return NULL; + iread = fread(A->ia, sizeof(int), A->m + 1, f); } - rc = fread(A->ja, sizeof(int), A->nz, f); - if (rc != A->nz) return NULL; + iread = fread(A->ja, sizeof(int), A->nz, f); if (size_of_matrix_type(A->type) > 0) { - rc = fread(A->a, size_of_matrix_type(A->type), A->nz, f); - if (rc != size_of_matrix_type(A->type)) return NULL; + iread = fread(A->a, size_of_matrix_type(A->type), A->nz, f); } fclose(f); return A; @@ -865,7 +876,11 @@ SparseMatrix SparseMatrix_from_coordinate_arrays_internal(int nz, int m, int n, return NULL; } A->nz = nz; + + + if(sum_repeated) A = SparseMatrix_sum_repeat_entries(A, sum_repeated); + return A; } @@ -897,7 +912,7 @@ SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B){ ic = C->ia; jc = C->ja; - mask = N_GNEW(n,int); + mask = MALLOC(sizeof(int)*n); for (i = 0; i < n; i++) mask[i] = -1; @@ -1014,7 +1029,7 @@ 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); + u = MALLOC(sizeof(real)*m*n); MEMCPY(u,v, sizeof(real)*m*n); for (i = 0; i < m; i++){ @@ -1042,7 +1057,7 @@ static void SparseMatrix_multiply_dense1(SparseMatrix A, real *v, real **res, in u = *res; if (!transposed){ - if (!u) u = N_GNEW(m*dim,real); + if (!u) u = MALLOC(sizeof(real)*m*dim); 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++){ @@ -1051,7 +1066,7 @@ static void SparseMatrix_multiply_dense1(SparseMatrix A, real *v, real **res, in } if (res_transposed) dense_transpose(u, m, dim); } else { - if (!u) u = N_GNEW(n*dim,real); + if (!u) u = MALLOC(sizeof(real)*n*dim); 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++){ @@ -1080,14 +1095,14 @@ static void SparseMatrix_multiply_dense2(SparseMatrix A, real *v, real **res, in n = A->n; if (!transposed){ - if (!u) u = N_GNEW(m*dim,real); + if (!u) u = MALLOC(sizeof(real)*m*dim); 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); + if (!u) u = MALLOC(sizeof(real)*n*dim); for (i = 0; i < dim; i++){ rr = &(u[n*i]); SparseMatrix_multiply_vector(A, &(v[m*i]), &rr, transposed); @@ -1143,7 +1158,7 @@ void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int trans if (v){ if (!transposed){ - if (!u) u = N_GNEW(m,real); + if (!u) u = MALLOC(sizeof(real)*m); for (i = 0; i < m; i++){ u[i] = 0.; for (j = ia[i]; j < ia[i+1]; j++){ @@ -1151,7 +1166,7 @@ void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int trans } } } else { - if (!u) u = N_GNEW(n,real); + if (!u) u = MALLOC(sizeof(real)*n); for (i = 0; i < n; i++) u[i] = 0.; for (i = 0; i < m; i++){ for (j = ia[i]; j < ia[i+1]; j++){ @@ -1162,7 +1177,7 @@ void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int trans } else { /* v is assumed to be all 1's */ if (!transposed){ - if (!u) u = N_GNEW(m,real); + if (!u) u = MALLOC(sizeof(real)*m); for (i = 0; i < m; i++){ u[i] = 0.; for (j = ia[i]; j < ia[i+1]; j++){ @@ -1170,7 +1185,7 @@ void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int trans } } } else { - if (!u) u = N_GNEW(n,real); + if (!u) u = MALLOC(sizeof(real)*n); for (i = 0; i < n; i++) u[i] = 0.; for (i = 0; i < m; i++){ for (j = ia[i]; j < ia[i+1]; j++){ @@ -1244,7 +1259,7 @@ SparseMatrix SparseMatrix_multiply_by_scaler(SparseMatrix A, real s){ SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B){ - int m, n; + int m; SparseMatrix C = NULL; int *mask = NULL; int *ia = A->ia, *ja = A->ja, *ib = B->ia, *jb = B->ja, *ic, *jc; @@ -1253,8 +1268,7 @@ SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B){ 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->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"); @@ -1263,10 +1277,10 @@ SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B){ } type = A->type; - mask = N_GNEW(n,int); + mask = MALLOC(sizeof(int)*(B->n)); if (!mask) return NULL; - for (i = 0; i < n; i++) mask[i] = -1; + for (i = 0; i < B->n; i++) mask[i] = -1; nz = 0; for (i = 0; i < m; i++){ @@ -1363,7 +1377,7 @@ SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B){ if (mask[jb[k]] < ic[i]){ mask[jb[k]] = nz; jc[nz] = jb[k]; - c[nz] += a[j]*b[k]; + c[nz] = a[j]*b[k]; nz++; } else { assert(jc[mask[jb[k]]] == jb[k]); @@ -1408,6 +1422,193 @@ SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B){ } + + +SparseMatrix SparseMatrix_multiply3(SparseMatrix A, SparseMatrix B, SparseMatrix C){ + int m; + SparseMatrix D = NULL; + int *mask = NULL; + int *ia = A->ia, *ja = A->ja, *ib = B->ia, *jb = B->ja, *ic = C->ia, *jc = C->ja, *id, *jd; + int i, j, k, l, ll, jj, type, nz; + + assert(A->format == B->format && A->format == FORMAT_CSR);/* other format not yet supported */ + + m = A->m; + if (A->n != B->m) return NULL; + if (B->n != C->m) return NULL; + + if (A->type != B->type || B->type != C->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 = MALLOC(sizeof(int)*(C->n)); + if (!mask) return NULL; + + for (i = 0; i < C->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 (l = ib[jj]; l < ib[jj+1]; l++){ + ll = jb[l]; + for (k = ic[ll]; k < ic[ll+1]; k++){ + if (mask[jc[k]] != -i - 2){ + if ((nz+1) <= nz) { +#ifdef DEBUG_PRINT + fprintf(stderr,"overflow in SparseMatrix_multiply !!!\n"); +#endif + return NULL; + } + nz++; + mask[jc[k]] = -i - 2; + } + } + } + } + } + + D = SparseMatrix_new(m, C->n, nz, type, FORMAT_CSR); + if (!D) goto RETURN; + id = D->ia; + jd = D->ja; + + nz = 0; + + switch (type){ + case MATRIX_TYPE_REAL: + { + real *a = (real*) A->a; + real *b = (real*) B->a; + real *c = (real*) C->a; + real *d = (real*) D->a; + id[0] = 0; + for (i = 0; i < m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + for (l = ib[jj]; l < ib[jj+1]; l++){ + ll = jb[l]; + for (k = ic[ll]; k < ic[ll+1]; k++){ + if (mask[jc[k]] < id[i]){ + mask[jc[k]] = nz; + jd[nz] = jc[k]; + d[nz] = a[j]*b[l]*c[k]; + nz++; + } else { + assert(jd[mask[jc[k]]] == jc[k]); + d[mask[jc[k]]] += a[j]*b[l]*c[k]; + } + } + } + } + id[i+1] = nz; + } + } + break; + case MATRIX_TYPE_COMPLEX: + { + real *a = (real*) A->a; + real *b = (real*) B->a; + real *c = (real*) C->a; + real *d = (real*) D->a; + id[0] = 0; + for (i = 0; i < m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + for (l = ib[jj]; l < ib[jj+1]; l++){ + ll = jb[l]; + for (k = ic[ll]; k < ic[ll+1]; k++){ + if (mask[jc[k]] < id[i]){ + mask[jc[k]] = nz; + jd[nz] = jc[k]; + d[2*nz] = (a[2*j]*b[2*l] - a[2*j+1]*b[2*l+1])*c[2*k] + - (a[2*j]*b[2*l+1] + a[2*j+1]*b[2*l])*c[2*k+1];/*real part */ + d[2*nz+1] = (a[2*j]*b[2*l+1] + a[2*j+1]*b[2*l])*c[2*k] + + (a[2*j]*b[2*l] - a[2*j+1]*b[2*l+1])*c[2*k+1];/*img part */ + nz++; + } else { + assert(jd[mask[jc[k]]] == jc[k]); + d[2*mask[jc[k]]] += (a[2*j]*b[2*l] - a[2*j+1]*b[2*l+1])*c[2*k] + - (a[2*j]*b[2*l+1] + a[2*j+1]*b[2*l])*c[2*k+1];/*real part */ + d[2*mask[jc[k]]+1] += (a[2*j]*b[2*l+1] + a[2*j+1]*b[2*l])*c[2*k] + + (a[2*j]*b[2*l] - a[2*j+1]*b[2*l+1])*c[2*k+1];/*img part */ + } + } + } + } + id[i+1] = nz; + } + } + break; + case MATRIX_TYPE_INTEGER: + { + int *a = (int*) A->a; + int *b = (int*) B->a; + int *c = (int*) C->a; + int *d = (int*) D->a; + id[0] = 0; + for (i = 0; i < m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + for (l = ib[jj]; l < ib[jj+1]; l++){ + ll = jb[l]; + for (k = ic[ll]; k < ic[ll+1]; k++){ + if (mask[jc[k]] < id[i]){ + mask[jc[k]] = nz; + jd[nz] = jc[k]; + d[nz] += a[j]*b[l]*c[k]; + nz++; + } else { + assert(jd[mask[jc[k]]] == jc[k]); + d[mask[jc[k]]] += a[j]*b[l]*c[k]; + } + } + } + } + id[i+1] = nz; + } + } + break; + case MATRIX_TYPE_PATTERN: + id[0] = 0; + for (i = 0; i < m; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + for (l = ib[jj]; l < ib[jj+1]; l++){ + ll = jb[l]; + for (k = ic[ll]; k < ic[ll+1]; k++){ + if (mask[jc[k]] < id[i]){ + mask[jc[k]] = nz; + jd[nz] = jc[k]; + nz++; + } else { + assert(jd[mask[jc[k]]] == jc[k]); + } + } + } + } + id[i+1] = nz; + } + break; + case MATRIX_TYPE_UNKNOWN: + default: + SparseMatrix_delete(D); + D = NULL; goto RETURN; + break; + } + + D->nz = nz; + + RETURN: + FREE(mask); + return D; + +} + /* For complex matrix: if what_to_sum = SUM_REPEATED_REAL_PART, we find entries {i,j,x + i y} and sum the x's if {i,j,Round(y)} are the same if what_to_sum = SUM_REPEATED_REAL_PART, we find entries {i,j,x + i y} and sum the y's if {i,j,Round(x)} are the same @@ -1432,7 +1633,7 @@ SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A, int what_to_sum){ if (what_to_sum == SUM_REPEATED_NONE) return A; - mask = N_GNEW(n,int); + mask = MALLOC(sizeof(int)*n); for (i = 0; i < n; i++) mask[i] = -1; switch (type){ @@ -1491,7 +1692,7 @@ SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A, int what_to_sum){ } } FREE(mask); - mask = N_GNEW(n*(ymax-ymin+1),int); + mask = MALLOC(sizeof(int)*n*(ymax-ymin+1)); for (i = 0; i < n*(ymax-ymin+1); i++) mask[i] = -1; nz = 0; @@ -1527,7 +1728,7 @@ SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A, int what_to_sum){ } } FREE(mask); - mask = N_GNEW(n*(xmax-xmin+1),int); + mask = MALLOC(sizeof(int)*n*(xmax-xmin+1)); for (i = 0; i < n*(xmax-xmin+1); i++) mask[i] = -1; nz = 0; @@ -1579,7 +1780,7 @@ SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A, int what_to_sum){ nz = 0; sta = ia[0]; for (i = 0; i < A->m; i++){ - for (j = ia[i]; j < ia[i+1]; j++){ + for (j = sta; j < ia[i+1]; j++){ if (mask[ja[j]] < ia[i]){ ja[nz] = ja[j]; mask[ja[j]] = nz++; @@ -1588,7 +1789,7 @@ SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A, int what_to_sum){ } } sta = ia[i+1]; - ia[i+1] = nz; + ia[i+1] = nz; } } break; @@ -1607,7 +1808,6 @@ SparseMatrix SparseMatrix_sum_repeat_entries(SparseMatrix A, int what_to_sum){ SparseMatrix SparseMatrix_coordinate_form_add_entries(SparseMatrix A, int nentries, int *irn, int *jcn, void *val){ int nz, type = A->type, nzmax, i; - assert(A); assert(A->format == FORMAT_COORD); if (nentries <= 0) return A; nz = A->nz; @@ -1615,7 +1815,7 @@ SparseMatrix SparseMatrix_coordinate_form_add_entries(SparseMatrix A, int nentri if (nz + nentries >= A->nzmax){ nzmax = nz + nentries; - A->nzmax = nzmax = MAX(10, (int) 0.2*nzmax) + 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); @@ -1870,7 +2070,7 @@ SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A){ A = SparseMatrix_symmetrize(B, TRUE); SparseMatrix_delete(B); A = SparseMatrix_remove_diagonal(A); - A->a = N_GNEW((A->nz),real); + A->a = MALLOC(sizeof(real)*(A->nz)); a = (real*) A->a; for (i = 0; i < A->nz; i++) a[i] = 1.; A->type = MATRIX_TYPE_REAL; @@ -1936,10 +2136,12 @@ SparseMatrix SparseMatrix_normalize_by_row(SparseMatrix 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 @@ -1948,6 +2150,7 @@ SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, double (*fun)(double x)){ return A; } + a = (real*) A->a; for (i = 0; i < A->m; i++){ for (j = A->ia[i]; j < A->ia[i+1]; j++){ @@ -1957,6 +2160,29 @@ SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, double (*fun)(double x)){ return A; } +SparseMatrix SparseMatrix_apply_fun_general(SparseMatrix A, void (*fun)(int i, int j, int n, double *x)){ + int i, j; + real *a; + int len = 1; + + if (!A) return A; + if (A->format != FORMAT_CSR || (A->type != MATRIX_TYPE_REAL&&A->type != MATRIX_TYPE_COMPLEX)) { +#ifdef DEBUG + printf("SparseMatrix_apply_fun: only CSR and real/complex matrix supported.\n"); +#endif + return A; + } + if (A->type == MATRIX_TYPE_COMPLEX) len = 2; + + a = (real*) A->a; + for (i = 0; i < A->m; i++){ + for (j = A->ia[i]; j < A->ia[i+1]; j++){ + fun(i, A->ja[j], len, &a[len*j]); + } + } + return A; +} + SparseMatrix SparseMatrix_crop(SparseMatrix A, real epsilon){ int i, j, *ia, *ja, nz, sta; @@ -2056,17 +2282,14 @@ void SparseMatrix_level_sets(SparseMatrix A, int root, int *nlevel, int **levels 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 (!(*levelset_ptr)) *levelset_ptr = MALLOC(sizeof(int)*(m+2)); + if (!(*levelset)) *levelset = MALLOC(sizeof(int)*m); 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; + *nlevel = 0; assert(root >= 0 && root < m); (*levelset_ptr)[0] = 0; (*levelset_ptr)[1] = 1; @@ -2102,7 +2325,7 @@ void SparseMatrix_weakly_connected_components(SparseMatrix A0, int *ncomp, int * if (!SparseMatrix_is_symmetric(A, TRUE)){ A = SparseMatrix_symmetrize(A, TRUE); } - if (!(*comps_ptr)) *comps_ptr = N_GNEW((m+1),int); + if (!(*comps_ptr)) *comps_ptr = MALLOC(sizeof(int)*(m+1)); *ncomp = 0; (*comps_ptr)[0] = 0; @@ -2118,11 +2341,203 @@ void SparseMatrix_weakly_connected_components(SparseMatrix A0, int *ncomp, int * } if (A != A0) SparseMatrix_delete(A); + if (levelset_ptr) FREE(levelset_ptr); + 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 */ + + +struct nodedata_struct { + real dist;/* distance to root */ + int id;/*node id */ +}; +typedef struct nodedata_struct* nodedata; + +static int cmp(void*i, void*j){ + nodedata d1, d2; + + d1 = (nodedata) i; + d2 = (nodedata) j; + if (d1->dist > d2->dist){ + return 1; + } else if (d1->dist == d2->dist){ + return 0; + } else { + return -1; + } +} + +static int Dijkstra(SparseMatrix A, int root, real *dist, int *nlist, int *list, real *dist_max){ + /* A: the nxn distance matrix. Entries are assumed to be nonnegative. Absolute value will be taken if + . entry value is negative. + dist: length n. On on exit contain the distance from root to every other node. dist[root] = 0 + . if the graph is disconnetced, unreachable node have a distance -1. + nlist: number of nodes visited + list: length n. the list of node in order of their extraction from the heap. The distance from root to last in the list should be the maximum + dist_max: the maximum distance, should be realized at node list[nlist-1]. + return: 0 if every node is reachable. -1 if not */ + + int m = A->m, i, j, jj, *ia = A->ia, *ja = A->ja, heap_id; + BinaryHeap h; + real *a = NULL, *aa; + int *ai; + nodedata ndata, ndata_min; + enum {UNVISITED = -2, FINISHED = -1}; + int *heap_ids; /* node ID to heap ID array. Initialised to UNVISITED. + Set to FINISHED after extracted as min from heap */ + int found = 0; + + assert(SparseMatrix_is_symmetric(A, TRUE)); + + assert(m == A->n); + + switch (A->type){ + case MATRIX_TYPE_COMPLEX: + aa = (real*) A->a; + a = MALLOC(sizeof(real)*(A->nz)); + for (i = 0; i < A->nz; i++) a[i] = aa[i*2]; + break; + case MATRIX_TYPE_REAL: + a = (real*) A->a; + break; + case MATRIX_TYPE_INTEGER: + ai = (int*) A->a; + a = MALLOC(sizeof(real)*(A->nz)); + for (i = 0; i < A->nz; i++) a[i] = (real) ai[i]; + break; + case MATRIX_TYPE_PATTERN: + a = MALLOC(sizeof(real)*A->nz); + for (i = 0; i < A->nz; i++) a[i] = 1.; + break; + default: + assert(0);/* no such matrix type */ + } + + heap_ids = MALLOC(sizeof(int)*m); + for (i = 0; i < m; i++) { + dist[i] = -1; + heap_ids[i] = UNVISITED; + } + + h = BinaryHeap_new(cmp); + assert(h); + + /* add root as the first item in the heap */ + ndata = MALLOC(sizeof(struct nodedata_struct)); + ndata->dist = 0; + ndata->id = root; + heap_ids[root] = BinaryHeap_insert(h, ndata); + + assert(heap_ids[root] >= 0);/* by design heap ID from BinaryHeap_insert >=0*/ + + while ((ndata_min = BinaryHeap_extract_min(h))){ + i = ndata_min->id; + dist[i] = ndata_min->dist; + list[found++] = i; + heap_ids[i] = FINISHED; + //fprintf(stderr," =================\n min extracted is id=%d, dist=%f\n",i, ndata_min->dist); + for (j = ia[i]; j < ia[i+1]; j++){ + jj = ja[j]; + heap_id = heap_ids[jj]; + if (jj == i || heap_id == FINISHED) continue; + if (heap_id == UNVISITED){ + ndata = MALLOC(sizeof(struct nodedata_struct)); + ndata->dist = ABS(a[j]) + ndata_min->dist; + ndata->id = jj; + heap_ids[jj] = BinaryHeap_insert(h, ndata); + //fprintf(stderr," set neighbor id=%d, dist=%f, hid = %d\n",jj, ndata->dist, heap_ids[jj]); + } else { + ndata = BinaryHeap_get_item(h, heap_id); + ndata->dist = MIN(ndata->dist, ABS(a[j]) + ndata_min->dist); + assert(ndata->id == jj); + BinaryHeap_reset(h, heap_id, ndata); + //fprintf(stderr," reset neighbor id=%d, dist=%f, hid = %d\n",jj, ndata->dist,heap_id); + } + } + FREE(ndata_min); + } + *nlist = found; + *dist_max = dist[i]; + + + BinaryHeap_delete(h, FREE); + FREE(heap_ids); + if (a && a != A->a) FREE(a); + if (found == m){ + return 0; + } else { + return -1; + } +} + + +real SparseMatrix_pseudo_diameter_weighted(SparseMatrix A0, int root, int aggressive, int *end1, int *end2, int *connectedQ){ + /* weighted graph. But still assume to be undirected. unsymmetric matrix ill be symmetrized */ + SparseMatrix A = A0; + int m = A->m, i, *list = NULL, nlist; + int flag; + real *dist = NULL, dist_max = -1, dist0 = -1; + int roots[5], iroots, end11, end22; + + if (!SparseMatrix_is_symmetric(A, TRUE)){ + A = SparseMatrix_symmetrize(A, TRUE); + } + assert(m == A->n); + + dist = MALLOC(sizeof(real)*m); + list = MALLOC(sizeof(int)*m); + nlist = 1; + list[nlist-1] = root; + + assert(SparseMatrix_is_symmetric(A, TRUE)); + + do { + dist0 = dist_max; + root = list[nlist - 1]; + flag = Dijkstra(A, root, dist, &nlist, list, &dist_max); + //fprintf(stderr,"after Dijkstra, {%d,%d}-%f\n",root, list[nlist-1], dist_max); + assert(dist[list[nlist-1]] == dist_max); + assert(root == list[0]); + assert(nlist > 0); + } while (dist_max > dist0); + + *connectedQ = (!flag); + assert((dist_max - dist0)/MAX(1, MAX(ABS(dist0), ABS(dist_max))) < 1.e-10); + + *end1 = root; + *end2 = list[nlist-1]; + //fprintf(stderr,"after search for diameter, diam = %f, ends = {%d,%d}\n", dist_max, *end1, *end2); + + if (aggressive){ + iroots = 0; + for (i = MAX(0, nlist - 6); i < nlist - 1; i++){ + roots[iroots++] = list[i]; + } + for (i = 0; i < iroots; i++){ + root = roots[i]; + dist0 = dist_max; + fprintf(stderr,"search for diameter again from root=%d\n", root); + dist_max = SparseMatrix_pseudo_diameter_weighted(A, root, FALSE, &end11, &end22, connectedQ); + if (dist_max > dist0){ + *end1 = end11; *end2 = end22; + } else { + dist_max = dist0; + } + } + fprintf(stderr,"after aggressive search for diameter, diam = %f, ends = {%d,%d}\n", dist_max, *end1, *end2); + } + + FREE(dist); + FREE(list); + + if (A != A0) SparseMatrix_delete(A); + return dist_max; + +} + +real SparseMatrix_pseudo_diameter_unweighted(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; @@ -2137,11 +2552,14 @@ int SparseMatrix_pseudo_diameter(SparseMatrix A0, int root, int aggressive, int assert(SparseMatrix_is_symmetric(A, TRUE)); SparseMatrix_level_sets(A, root, &nlevel, &levelset_ptr, &levelset, &mask, TRUE); + // fprintf(stderr,"after level set, {%d,%d}=%d\n",levelset[0], levelset[levelset_ptr[nlevel]-1], nlevel); + *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); + //fprintf(stderr,"after level set, {%d,%d}=%d\n",levelset[0], levelset[levelset_ptr[nlevel]-1], nlevel); } *end1 = levelset[0]; *end2 = levelset[levelset_ptr[nlevel]-1]; @@ -2149,13 +2567,13 @@ int SparseMatrix_pseudo_diameter(SparseMatrix A0, int root, int aggressive, int if (aggressive){ nlevel0 = nlevel; iroots = 0; - for (i = levelset_ptr[nlevel]; i < MIN(levelset_ptr[nlevel + 1], levelset_ptr[nlevel]+5); i++){ + for (i = levelset_ptr[nlevel-1]; i < MIN(levelset_ptr[nlevel], levelset_ptr[nlevel-1]+5); i++){ iroots++; - roots[i - levelset_ptr[nlevel]] = levelset[levelset_ptr[i]]; + roots[i - levelset_ptr[nlevel-1]] = levelset[i]; } for (i = 0; i < iroots; i++){ root = roots[i]; - nlevel = SparseMatrix_pseudo_diameter(A, root, FALSE, &enda, &endb, connectedQ); + nlevel = (int) SparseMatrix_pseudo_diameter_unweighted(A, root, FALSE, &enda, &endb, connectedQ); if (nlevel > nlevel0) { nlevel0 = nlevel; *end1 = enda; @@ -2168,12 +2586,12 @@ int SparseMatrix_pseudo_diameter(SparseMatrix A0, int root, int aggressive, int FREE(levelset); FREE(mask); if (A != A0) SparseMatrix_delete(A); - return nlevel0; + return (real) nlevel0 - 1; } -int SparseMatrix_pseudo_diameter_only(SparseMatrix A){ +real SparseMatrix_pseudo_diameter_only(SparseMatrix A){ int end1, end2, connectedQ; - return SparseMatrix_pseudo_diameter(A, 0, FALSE, &end1, &end2, &connectedQ); + return SparseMatrix_pseudo_diameter_unweighted(A, 0, FALSE, &end1, &end2, &connectedQ); } int SparseMatrix_connectedQ(SparseMatrix A0){ @@ -2204,10 +2622,10 @@ void SparseMatrix_decompose_to_supervariables(SparseMatrix A, int *ncluster, int 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); + super = MALLOC(sizeof(int)*(n)); + nsuper = MALLOC(sizeof(int)*(n+1)); + mask = MALLOC(sizeof(int)*n); + newmap = MALLOC(sizeof(int)*n); nsuper++; isup = 0; @@ -2305,8 +2723,8 @@ SparseMatrix SparseMatrix_get_augmented(SparseMatrix A){ SparseMatrix B = NULL; if (!A) return NULL; if (nz > 0){ - irn = N_GNEW(nz*2,int); - jcn = N_GNEW(nz*2,int); + irn = MALLOC(sizeof(int)*nz*2); + jcn = MALLOC(sizeof(int)*nz*2); } if (A->a){ @@ -2340,10 +2758,351 @@ SparseMatrix SparseMatrix_get_augmented(SparseMatrix A){ } -SparseMatrix SparseMatrix_to_square_matrix(SparseMatrix A){ +SparseMatrix SparseMatrix_to_square_matrix(SparseMatrix A, int bipartite_options){ SparseMatrix B; - if (A->m == A->n) return A; + switch (bipartite_options){ + case BIPARTITE_RECT: + if (A->m == A->n) return A; + break; + case BIPARTITE_PATTERN_UNSYM: + if (A->m == A->n && SparseMatrix_is_symmetric(A, TRUE)) return A; + break; + case BIPARTITE_UNSYM: + if (A->m == A->n && SparseMatrix_is_symmetric(A, FALSE)) return A; + break; + case BIPARTITE_ALWAYS: + break; + default: + assert(0); + } B = SparseMatrix_get_augmented(A); SparseMatrix_delete(A); return B; } + +SparseMatrix SparseMatrix_get_submatrix(SparseMatrix A, int nrow, int ncol, int *rindices, int *cindices){ + /* get the submatrix from row/columns indices[0,...,l-1]. + row rindices[i] will be the new row i + column cindices[i] will be the new column i. + if rindices = NULL, it is assume that 1 -- nrow is needed. Same for cindices/ncol. + */ + int nz = 0, i, j, *irn, *jcn, *ia = A->ia, *ja = A->ja, m = A->m, n = A->n; + int *cmask, *rmask; + void *v = NULL; + SparseMatrix B = NULL; + int irow = 0, icol = 0; + + if (nrow <= 0 || ncol <= 0) return NULL; + + + + rmask = MALLOC(sizeof(int)*m); + cmask = MALLOC(sizeof(int)*n); + for (i = 0; i < m; i++) rmask[i] = -1; + for (i = 0; i < n; i++) cmask[i] = -1; + + if (rindices){ + for (i = 0; i < nrow; i++) { + if (rindices[i] >= 0 && rindices[i] < m){ + rmask[rindices[i]] = irow++; + } + } + } else { + for (i = 0; i < nrow; i++) { + rmask[i] = irow++; + } + } + + if (cindices){ + for (i = 0; i < ncol; i++) { + if (cindices[i] >= 0 && cindices[i] < n){ + cmask[cindices[i]] = icol++; + } + } + } else { + for (i = 0; i < ncol; i++) { + cmask[i] = icol++; + } + } + + for (i = 0; i < m; i++){ + if (rmask[i] < 0) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (cmask[ja[j]] < 0) continue; + nz++; + } + } + + + switch (A->type){ + case MATRIX_TYPE_REAL:{ + real *a = (real*) A->a; + real *val; + irn = MALLOC(sizeof(int)*nz); + jcn = MALLOC(sizeof(int)*nz); + val = MALLOC(sizeof(real)*nz); + + nz = 0; + for (i = 0; i < m; i++){ + if (rmask[i] < 0) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (cmask[ja[j]] < 0) continue; + irn[nz] = rmask[i]; + jcn[nz] = cmask[ja[j]]; + val[nz++] = a[j]; + } + } + v = (void*) val; + break; + } + case MATRIX_TYPE_COMPLEX:{ + real *a = (real*) A->a; + real *val; + + irn = MALLOC(sizeof(int)*nz); + jcn = MALLOC(sizeof(int)*nz); + val = MALLOC(sizeof(real)*2*nz); + + nz = 0; + for (i = 0; i < m; i++){ + if (rmask[i] < 0) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (cmask[ja[j]] < 0) continue; + irn[nz] = rmask[i]; + jcn[nz] = cmask[ja[j]]; + val[2*nz] = a[2*j]; + val[2*nz+1] = a[2*j+1]; + nz++; + } + } + v = (void*) val; + break; + } + case MATRIX_TYPE_INTEGER:{ + int *a = (int*) A->a; + int *val; + + irn = MALLOC(sizeof(int)*nz); + jcn = MALLOC(sizeof(int)*nz); + val = MALLOC(sizeof(int)*nz); + + nz = 0; + for (i = 0; i < m; i++){ + if (rmask[i] < 0) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (cmask[ja[j]] < 0) continue; + irn[nz] = rmask[i]; + jcn[nz] = cmask[ja[j]]; + val[nz] = a[j]; + nz++; + } + } + v = (void*) val; + break; + } + case MATRIX_TYPE_PATTERN: + irn = MALLOC(sizeof(int)*nz); + jcn = MALLOC(sizeof(int)*nz); + nz = 0; + for (i = 0; i < m; i++){ + if (rmask[i] < 0) continue; + for (j = ia[i]; j < ia[i+1]; j++){ + if (cmask[ja[j]] < 0) continue; + irn[nz] = rmask[i]; + jcn[nz++] = cmask[ja[j]]; + } + } + break; + case MATRIX_TYPE_UNKNOWN: + FREE(rmask); + FREE(cmask); + return NULL; + default: + FREE(rmask); + FREE(cmask); + return NULL; + } + + B = SparseMatrix_from_coordinate_arrays(nz, nrow, ncol, irn, jcn, v, A->type); + FREE(cmask); + FREE(rmask); + FREE(irn); + FREE(jcn); + if (v) FREE(v); + + + return B; + +} + +SparseMatrix SparseMatrix_exclude_submatrix(SparseMatrix A, int nrow, int ncol, int *rindices, int *cindices){ + /* get a submatrix by excluding rows and columns */ + int *r, *c, nr, nc, i; + SparseMatrix B; + + if (nrow <= 0 && ncol <= 0) return A; + + r = MALLOC(sizeof(int)*A->m); + c = MALLOC(sizeof(int)*A->n); + + for (i = 0; i < A->m; i++) r[i] = i; + for (i = 0; i < A->n; i++) c[i] = i; + for (i = 0; i < nrow; i++) { + if (rindices[i] >= 0 && rindices[i] < A->m){ + r[rindices[i]] = -1; + } + } + for (i = 0; i < ncol; i++) { + if (cindices[i] >= 0 && cindices[i] < A->n){ + c[cindices[i]] = -1; + } + } + + nr = nc = 0; + for (i = 0; i < A->m; i++) { + if (r[i] > 0) r[nr++] = r[i]; + } + for (i = 0; i < A->n; i++) { + if (c[i] > 0) c[nc++] = c[i]; + } + + B = SparseMatrix_get_submatrix(A, nr, nc, r, c); + + FREE(r); + FREE(c); + return B; + +} + +SparseMatrix SparseMatrix_largest_component(SparseMatrix A){ + SparseMatrix B; + int ncomp; + int *comps = NULL; + int *comps_ptr = NULL; + int i; + int nmax, imax = 0; + + if (!A) return NULL; + A = SparseMatrix_to_square_matrix(A, BIPARTITE_RECT); + SparseMatrix_weakly_connected_components(A, &ncomp, &comps, &comps_ptr); + if (ncomp == 1) { + B = A; + } else { + nmax = 0; + for (i = 0; i < ncomp; i++){ + if (nmax < comps_ptr[i+1] - comps_ptr[i]){ + nmax = comps_ptr[i+1] - comps_ptr[i]; + imax = i; + } + } + B = SparseMatrix_get_submatrix(A, nmax, nmax, &comps[comps_ptr[imax]], &comps[comps_ptr[imax]]); + } + FREE(comps); + FREE(comps_ptr); + return B; + + +} + +SparseMatrix SparseMatrix_delete_empty_columns(SparseMatrix A, int **new2old, int *nnew, int inplace){ + /* delete empty columns in A. After than number of columns will be nnew, and + the mapping from new matrix column to old matrix column is new2old. + On entry, if new2old is NULL, it is allocated. + */ + SparseMatrix B; + int *ia, *ja; + int *old2new; + int i; + old2new = MALLOC(sizeof(int)*A->n); + for (i = 0; i < A->n; i++) old2new[i] = -1; + + *nnew = 0; + B = SparseMatrix_transpose(A); + ia = B->ia; ja = B->ja; + for (i = 0; i < B->m; i++){ + if (ia[i+1] > ia[i]){ + (*nnew)++; + } + } + if (!(*new2old)) *new2old = MALLOC(sizeof(int)*(*nnew)); + + *nnew = 0; + for (i = 0; i < B->m; i++){ + if (ia[i+1] > ia[i]){ + (*new2old)[*nnew] = i; + old2new[i] = *nnew; + (*nnew)++; + } + } + SparseMatrix_delete(B); + + if (inplace){ + B = A; + } else { + B = SparseMatrix_copy(A); + } + ia = B->ia; ja = B->ja; + for (i = 0; i < ia[B->m]; i++){ + assert(old2new[ja[i]] >= 0); + ja[i] = old2new[ja[i]]; + } + B->n = *nnew; + + FREE(old2new); + return B; + + +} + +SparseMatrix SparseMatrix_set_entries_to_real_one(SparseMatrix A){ + real *a; + int i; + + if (A->a) FREE(A->a); + A->a = MALLOC(sizeof(real)*A->nz); + a = (real*) (A->a); + for (i = 0; i < A->nz; i++) a[i] = 1.; + A->type = MATRIX_TYPE_REAL; + return A; + +} + +SparseMatrix SparseMatrix_complement(SparseMatrix A, int undirected){ + /* find the complement graph A^c, such that {i,h}\in E(A_c) iff {i,j} \notin E(A). Only structural matrix is returned. */ + SparseMatrix B = A; + int *ia, *ja; + int m = A->m, n = A->n; + int *mask, nz = 0; + int *irn, *jcn; + int i, j; + + if (undirected) B = SparseMatrix_symmetrize(A, TRUE); + assert(m == n); + + ia = B->ia; ja = B->ja; + mask = MALLOC(sizeof(int)*n); + irn = MALLOC(sizeof(int)*(n*n - A->nz)); + jcn = MALLOC(sizeof(int)*(n*n - A->nz)); + + for (i = 0; i < n; i++){ + mask[i] = -1; + } + + for (i = 0; i < n; i++){ + for (j = ia[i]; j < ia[i+1]; j++){ + mask[ja[j]] = i; + } + for (j = 0; j < n; j++){ + if (mask[j] != i){ + irn[nz] = i; + jcn[nz++] = j; + } + } + } + + if (B != A) SparseMatrix_delete(B); + B = SparseMatrix_from_coordinate_arrays(nz, m, n, irn, jcn, NULL, MATRIX_TYPE_PATTERN); + FREE(irn); + FREE(jcn); + return B; +} diff --git a/lib/sparse/SparseMatrix.h b/lib/sparse/SparseMatrix.h index 7ea07b65b..e511dbd54 100644 --- a/lib/sparse/SparseMatrix.h +++ b/lib/sparse/SparseMatrix.h @@ -1,3 +1,4 @@ +/* $Id$Revision: */ /* vim:set shiftwidth=4 ts=8: */ /************************************************************************* @@ -11,17 +12,15 @@ *************************************************************************/ #ifndef SPARSEMATRIX_H #define SPARSEMATRIX_H -#define set_flag(a, flag) ((a)=((a)|(flag))) -#define test_flag(a, flag) ((a)&(flag)) -#define clear_flag(a, flag) ((a) &=(~(flag))) +#include #include -typedef double real; #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}; +enum {BIPARTITE_RECT = 0, BIPARTITE_PATTERN_UNSYM, BIPARTITE_UNSYM, BIPARTITE_ALWAYS}; struct SparseMatrix_struct { int m; /* row dimension */ @@ -61,10 +60,11 @@ void SparseMatrix_delete(SparseMatrix A); SparseMatrix SparseMatrix_add(SparseMatrix A, SparseMatrix B); SparseMatrix SparseMatrix_multiply(SparseMatrix A, SparseMatrix B); +SparseMatrix SparseMatrix_multiply3(SparseMatrix A, SparseMatrix B, SparseMatrix C); /* For complex matrix: if what_to_sum = SUM_REPEATED_REAL_PART, we find entries {i,j,x + i y} and sum the x's if {i,j,Round(y)} are the same - if what_to_sum = SUM_REPEATED_REAL_PART, we find entries {i,j,x + i y} and sum the y's if {i,j,Round(x)} are the same + if what_to_sum = SUM_REPEATED_IMAGINARY_PART, we find entries {i,j,x + i y} and sum the y's if {i,j,Round(x)} are the same For other matrix, what_to_sum = SUM_REPEATED_REAL_PART is the same as what_to_sum = SUM_REPEATED_IMAGINARY_PART or what_to_sum = SUM_REPEATED_ALL */ @@ -74,6 +74,7 @@ SparseMatrix SparseMatrix_coordinate_form_add_entries(SparseMatrix A, int nentri 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); +SparseMatrix SparseMatrix_symmetrize_nodiag(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 */ @@ -82,6 +83,7 @@ SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A); 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_apply_fun_general(SparseMatrix A, void (*fun)(int i, int j, int n, double *x));/* for real and complex (n=2) */ SparseMatrix SparseMatrix_copy(SparseMatrix A); int SparseMatrix_has_diagonal(SparseMatrix A); SparseMatrix SparseMatrix_normalize_by_row(SparseMatrix A);/* divide by max of each row */ @@ -90,14 +92,34 @@ SparseMatrix SparseMatrix_scaled_by_vector(SparseMatrix A, real *v, int apply_to 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 */ +real SparseMatrix_pseudo_diameter_only(SparseMatrix A); +real SparseMatrix_pseudo_diameter_weighted(SparseMatrix A0, int root, int aggressive, int *end1, int *end2, int *connectedQ); /* assume real distances, unsymmetric matrix ill be symmetrized */ +real SparseMatrix_pseudo_diameter_unweighted(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); +SparseMatrix SparseMatrix_get_submatrix(SparseMatrix A, int nrow, int ncol, int *rindices, int *cindices); +SparseMatrix SparseMatrix_exclude_submatrix(SparseMatrix A, int nrow, int ncol, int *rindices, int *cindices); SparseMatrix SparseMatrix_get_augmented(SparseMatrix A); -SparseMatrix SparseMatrix_to_square_matrix(SparseMatrix A); + +/* bipartite_options: + BIPARTITE_RECT -- turn rectangular matrix into square), + BIPARTITE_PATTERN_UNSYM -- pattern unsummetric as bipartite + BIPARTITE_UNSYM -- unsymmetric as square + BIPARTITE_ALWAYS -- always as square +*/ +SparseMatrix SparseMatrix_to_square_matrix(SparseMatrix A, int bipartite_options); + +SparseMatrix SparseMatrix_largest_component(SparseMatrix A); + +SparseMatrix SparseMatrix_delete_empty_columns(SparseMatrix A, int **new2old, int *nnew, int inplace); + +SparseMatrix SparseMatrix_sort(SparseMatrix A); + +SparseMatrix SparseMatrix_set_entries_to_real_one(SparseMatrix A); + +SparseMatrix SparseMatrix_complement(SparseMatrix A, int undirected); #define SparseMatrix_set_undirected(A) set_flag((A)->property, MATRIX_UNDIRECTED) #define SparseMatrix_set_symmetric(A) set_flag((A)->property, MATRIX_SYMMETRIC) diff --git a/lib/sparse/general.c b/lib/sparse/general.c new file mode 100644 index 000000000..1acac9821 --- /dev/null +++ b/lib/sparse/general.c @@ -0,0 +1,300 @@ +/* $Id$Revision: */ +/* vim:set shiftwidth=4 ts=8: */ + +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#include "general.h" + +#ifdef DEBUG +double _statistics[10]; +#endif + +real drand(){ + return rand()/(real) RAND_MAX; +} + +int irand(int n){ + /* 0, 1, ..., n-1 */ + assert(n > 1); + /*return (int) MIN(floor(drand()*n),n-1);*/ + return rand()%n; +} + +int *random_permutation(int n){ + int *p; + int i, j, pp, len; + if (n <= 0) return NULL; + p = MALLOC(sizeof(int)*n); + 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; +} + + +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; +} + +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* 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; +} + +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; +} + +void vector_print(char *s, int n, real *x){ + int i; + printf("%s{",s); + for (i = 0; i < n; i++) { + if (i > 0) printf(","); + printf("%f",x[i]); + } + printf("}\n"); +} + +void vector_take(int n, real *v, int m, int *p, real **u){ + /* take m elements v[p[i]]],i=1,...,m and oput in u */ + int i; + + if (!*u) *u = MALLOC(sizeof(real)*m); + + for (i = 0; i < m; i++) { + assert(p[i] < n && p[i] >= 0); + (*u)[i] = v[p[i]]; + } + +} + +void vector_float_take(int n, float *v, int m, int *p, float **u){ + /* take m elements v[p[i]]],i=1,...,m and oput in u */ + int i; + + if (!*u) *u = MALLOC(sizeof(float)*m); + + for (i = 0; i < m; i++) { + assert(p[i] < n && p[i] >= 0); + (*u)[i] = v[p[i]]; + } + +} + +int comp_descend(const void *s1, const void *s2){ + real *ss1, *ss2; + ss1 = (real*) s1; + ss2 = (real*) s2; + + if ((ss1)[0] > (ss2)[0]){ + return -1; + } else if ((ss1)[0] < (ss2)[0]){ + return 1; + } + return 0; +} + +int comp_ascend(const void *s1, const void *s2){ + real *ss1, *ss2; + ss1 = (real*) s1; + ss2 = (real*) s2; + + if ((ss1)[0] > (ss2)[0]){ + return -1; + } else if ((ss1)[0] < (ss2)[0]){ + return 1; + } + return 0; +} +int comp_descend_int(const void *s1, const void *s2){ + int *ss1, *ss2; + ss1 = (int*) s1; + ss2 = (int*) s2; + + if ((ss1)[0] > (ss2)[0]){ + return -1; + } else if ((ss1)[0] < (ss2)[0]){ + return 1; + } + return 0; +} + +int comp_ascend_int(const void *s1, const void *s2){ + int *ss1, *ss2; + ss1 = (int*) s1; + ss2 = (int*) s2; + + if ((ss1)[0] > (ss2)[0]){ + return -1; + } else if ((ss1)[0] < (ss2)[0]){ + return 1; + } + return 0; +} + + +void vector_ordering(int n, real *v, int **p, int ascending){ + /* give the position of the lagest, second largest etc in vector v if ascending = TRUE + + or + + give the position of the smallest, second smallest etc in vector v if ascending = TRUE. + results in p. If *p == NULL, p is asigned. + + */ + + real *u; + int i; + + if (!*p) *p = MALLOC(sizeof(int)*n); + u = MALLOC(sizeof(real)*2*n); + + for (i = 0; i < n; i++) { + u[2*i+1] = i; + u[2*i] = v[i]; + } + + if (ascending){ + qsort(u, n, sizeof(real)*2, comp_ascend); + } else { + qsort(u, n, sizeof(real)*2, comp_descend); + } + + for (i = 0; i < n; i++) (*p)[i] = (int) u[2*i+1]; + FREE(u); + +} + +void vector_sort_real(int n, real *v, int ascending){ + if (ascending){ + qsort(v, n, sizeof(real), comp_ascend); + } else { + qsort(v, n, sizeof(real), comp_descend); + } +} +void vector_sort_int(int n, int *v, int ascending){ + if (ascending){ + qsort(v, n, sizeof(int), comp_ascend_int); + } else { + qsort(v, n, sizeof(int), comp_descend_int); + } +} + +int excute_system_command3(char *s1, char *s2, char *s3){ + char c[1000]; + + strcpy(c, s1); + strcat(c, s2); + strcat(c, s3); + return system(c); +} + +int excute_system_command(char *s1, char *s2){ + char c[1000]; + + strcpy(c, s1); + strcat(c, s2); + return system(c); +} + + +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; +} + +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); +} + +char *strip_dir(char *s){ + int i, first = TRUE; + for (i = strlen(s); i >= 0; i--) { + if (first && s[i] == '.') {/* get rid of .mtx */ + s[i] = '\0'; + first = FALSE; + } + if (s[i] == '/') return (char*) &(s[i+1]); + } + return s; +} + +void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x){ + real min[3], max[3], min0[3], ratio = 1; + int i, k; + + for (i = 0; i < dim; i++) { + min[i] = x[i]; + max[i] = x[i]; + } + + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++) { + min[k] = MIN(x[i*dim+k], min[k]); + max[k] = MAX(x[i*dim+k], max[k]); + } + } + + if (max[0] - min[0] != 0) { + ratio = (xmax-xmin)/(max[0] - min[0]); + } + if (max[1] - min[1] != 0) { + ratio = MIN(ratio, (ymax-ymin)/(max[1] - min[1])); + } + + min0[0] = xmin; + min0[1] = ymin; + min0[2] = 0; + for (i = 0; i < n; i++){ + for (k = 0; k < dim; k++) { + x[i*dim+k] = min0[k] + (x[i*dim+k] - min[k])*ratio; + } + } + + +} diff --git a/lib/sparse/general.h b/lib/sparse/general.h new file mode 100644 index 000000000..6968c0654 --- /dev/null +++ b/lib/sparse/general.h @@ -0,0 +1,141 @@ +/* $Id$Revision: */ +/* vim:set shiftwidth=4 ts=8: */ + +/************************************************************************* + * Copyright (c) 2011 AT&T Intellectual Property + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License v1.0 + * which accompanies this distribution, and is available at + * http://www.eclipse.org/legal/epl-v10.html + * + * Contributors: See CVS logs. Details at http://www.graphviz.org/ + *************************************************************************/ + +#ifndef GENERAL_H +#define GENERAL_H + +#include +#include +#include +#include +#include +/* Applications that do not use the common library can define STANDALONE + * to get definitions/definitions that are normally provided there. + * In particular, note that Verbose is declared but undefined. + */ +#ifndef STANDALONE +#include +#include +#include +#endif + +#define real double + +#define set_flag(a, flag) ((a)=((a)|(flag))) +#define test_flag(a, flag) ((a)&(flag)) +#define clear_flag(a, flag) ((a) &=(~(flag))) + +#ifdef STANDALONE +#define MALLOC malloc +#define REALLOC realloc +#else +#define MALLOC gmalloc +#define REALLOC grealloc +#endif +#define FREE free +#define MEMCPY memcpy + +#ifdef STANDALONE +#define N_NEW(n,t) (t*)malloc((n)*sizeof(t)) +#define NEW(t) (t*)malloc(sizeof(t)) +#endif + +#ifndef DEBUG +#define NDEBUG /* switch off assert*/ +#endif + +#ifdef STANDALONE +#define MAX(a,b) ((a)>(b)?(a):b) +#define MIN(a,b) ((a)<(b)?(a):b) +#define ABS(a) (((a)>0)?(a):(-(a))) + +#ifdef TRUE +#undef TRUE +#endif +#define TRUE 1 + +#ifdef FALSE +#undef FALSE +#endif +#define FALSE 0 + +extern int Verbose; +#endif /* STANDALONE */ + +#ifdef DEBUG +extern double _statistics[10]; +#endif + +int irand(int n); +real drand(); +int *random_permutation(int n);/* random permutation of 0 to n-1 */ + +#ifdef STANDALONE +#define MAXINT 1<<30 +#endif + +#define PI 3.14159 + +real* vector_subtract_to(int n, real *x, real *y); + +real vector_product(int n, real *x, real *y); + +real* vector_saxpy(int n, real *x, real *y, real beta); /* y = x+beta*y */ + + +real* vector_saxpy2(int n, real *x, real *y, real beta);/* x = x+beta*y */ + +/* take m elements v[p[i]]],i=1,...,m and oput in u. u will be assigned if *u = NULL */ +void vector_take(int n, real *v, int m, int *p, real **u); +void vector_float_take(int n, float *v, int m, int *p, float **u); + +/* give the position of the lagest, second largest etc in vector v if ascending = TRUE + or + give the position of the smallest, second smallest etc in vector v if ascending = TRUE. + results in p. If *p == NULL, p is asigned. +*/ +void vector_ordering(int n, real *v, int **p, int ascending); +void vector_sort_real(int n, real *v, int ascending); +void vector_sort_int(int n, int *v, int ascending); + +void vector_print(char *s, int n, real *x); + +#define MACHINEACC 1.0e-16 + +#ifdef STANDALONE +#define POINTS(inch) 72*(inch) + +typedef unsigned int boolean; +#endif /* STANDALONE */ + +int excute_system_command3(char *s1, char *s2, char *s3); +int excute_system_command(char *s1, char *s2); + +#define MINDIST 1.e-15 + +enum {UNMATCHED = -1}; + + +real distance(real *x, int dim, int i, int j); +real distance_cropped(real *x, int dim, int i, int j); + +real point_distance(real *p1, real *p2, int dim); + +char *strip_dir(char *s); + +void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x); + +#endif + + + -- 2.40.0