From b6ede4eb01bcab129be9ca7be102174b4cd4418b Mon Sep 17 00:00:00 2001 From: Allanis Date: Thu, 10 Apr 2014 19:06:30 +0100 Subject: [PATCH] [Add] Added CSparse lib for linear algebra. --- bin/Makefile | 51 ++++-- lib/csparse/Makefile | 37 ++++ lib/csparse/cs.h | 143 +++++++++++++++ lib/csparse/cs_add.c | 28 +++ lib/csparse/cs_amd.c | 364 +++++++++++++++++++++++++++++++++++++ lib/csparse/cs_chol.c | 58 ++++++ lib/csparse/cs_cholsol.c | 26 +++ lib/csparse/cs_compress.c | 22 +++ lib/csparse/cs_counts.c | 61 +++++++ lib/csparse/cs_cumsum.c | 17 ++ lib/csparse/cs_dfs.c | 36 ++++ lib/csparse/cs_dmperm.c | 144 +++++++++++++++ lib/csparse/cs_droptol.c | 9 + lib/csparse/cs_dropzeros.c | 9 + lib/csparse/cs_dupl.c | 34 ++++ lib/csparse/cs_entry.c | 13 ++ lib/csparse/cs_ereach.c | 23 +++ lib/csparse/cs_etree.c | 30 +++ lib/csparse/cs_fkeep.c | 25 +++ lib/csparse/cs_gaxpy.c | 17 ++ lib/csparse/cs_happly.c | 19 ++ lib/csparse/cs_house.c | 23 +++ lib/csparse/cs_ipvec.c | 9 + lib/csparse/cs_leaf.c | 22 +++ lib/csparse/cs_load.c | 15 ++ lib/csparse/cs_lsolve.c | 18 ++ lib/csparse/cs_ltsolve.c | 18 ++ lib/csparse/cs_lu.c | 86 +++++++++ lib/csparse/cs_lusol.c | 26 +++ lib/csparse/cs_malloc.c | 35 ++++ lib/csparse/cs_maxtrans.c | 92 ++++++++++ lib/csparse/cs_multiply.c | 35 ++++ lib/csparse/cs_norm.c | 15 ++ lib/csparse/cs_permute.c | 25 +++ lib/csparse/cs_pinv.c | 11 ++ lib/csparse/cs_post.c | 24 +++ lib/csparse/cs_print.c | 36 ++++ lib/csparse/cs_pvec.c | 9 + lib/csparse/cs_qr.c | 73 ++++++++ lib/csparse/cs_qrsol.c | 53 ++++++ lib/csparse/cs_randperm.c | 22 +++ lib/csparse/cs_reach.c | 19 ++ lib/csparse/cs_scatter.c | 22 +++ lib/csparse/cs_scc.c | 41 +++++ lib/csparse/cs_schol.c | 26 +++ lib/csparse/cs_spsolve.c | 28 +++ lib/csparse/cs_sqr.c | 88 +++++++++ lib/csparse/cs_symperm.c | 39 ++++ lib/csparse/cs_tdfs.c | 24 +++ lib/csparse/cs_transpose.c | 25 +++ lib/csparse/cs_updown.c | 37 ++++ lib/csparse/cs_usolve.c | 18 ++ lib/csparse/cs_util.c | 119 ++++++++++++ lib/csparse/cs_utsolve.c | 18 ++ 54 files changed, 2297 insertions(+), 20 deletions(-) create mode 100644 lib/csparse/Makefile create mode 100644 lib/csparse/cs.h create mode 100644 lib/csparse/cs_add.c create mode 100644 lib/csparse/cs_amd.c create mode 100644 lib/csparse/cs_chol.c create mode 100644 lib/csparse/cs_cholsol.c create mode 100644 lib/csparse/cs_compress.c create mode 100644 lib/csparse/cs_counts.c create mode 100644 lib/csparse/cs_cumsum.c create mode 100644 lib/csparse/cs_dfs.c create mode 100644 lib/csparse/cs_dmperm.c create mode 100644 lib/csparse/cs_droptol.c create mode 100644 lib/csparse/cs_dropzeros.c create mode 100644 lib/csparse/cs_dupl.c create mode 100644 lib/csparse/cs_entry.c create mode 100644 lib/csparse/cs_ereach.c create mode 100644 lib/csparse/cs_etree.c create mode 100644 lib/csparse/cs_fkeep.c create mode 100644 lib/csparse/cs_gaxpy.c create mode 100644 lib/csparse/cs_happly.c create mode 100644 lib/csparse/cs_house.c create mode 100644 lib/csparse/cs_ipvec.c create mode 100644 lib/csparse/cs_leaf.c create mode 100644 lib/csparse/cs_load.c create mode 100644 lib/csparse/cs_lsolve.c create mode 100644 lib/csparse/cs_ltsolve.c create mode 100644 lib/csparse/cs_lu.c create mode 100644 lib/csparse/cs_lusol.c create mode 100644 lib/csparse/cs_malloc.c create mode 100644 lib/csparse/cs_maxtrans.c create mode 100644 lib/csparse/cs_multiply.c create mode 100644 lib/csparse/cs_norm.c create mode 100644 lib/csparse/cs_permute.c create mode 100644 lib/csparse/cs_pinv.c create mode 100644 lib/csparse/cs_post.c create mode 100644 lib/csparse/cs_print.c create mode 100644 lib/csparse/cs_pvec.c create mode 100644 lib/csparse/cs_qr.c create mode 100644 lib/csparse/cs_qrsol.c create mode 100644 lib/csparse/cs_randperm.c create mode 100644 lib/csparse/cs_reach.c create mode 100644 lib/csparse/cs_scatter.c create mode 100644 lib/csparse/cs_scc.c create mode 100644 lib/csparse/cs_schol.c create mode 100644 lib/csparse/cs_spsolve.c create mode 100644 lib/csparse/cs_sqr.c create mode 100644 lib/csparse/cs_symperm.c create mode 100644 lib/csparse/cs_tdfs.c create mode 100644 lib/csparse/cs_transpose.c create mode 100644 lib/csparse/cs_updown.c create mode 100644 lib/csparse/cs_usolve.c create mode 100644 lib/csparse/cs_util.c create mode 100644 lib/csparse/cs_utsolve.c diff --git a/bin/Makefile b/bin/Makefile index ef3679f..4ee1147 100644 --- a/bin/Makefile +++ b/bin/Makefile @@ -25,13 +25,14 @@ APPNAME := Lephisto OBJS := $(patsubst %.c, %.o, $(wildcard ../src/*.c)) # CFLAGS -CLUA := -I../lib/lua -CSDL := $(shell sdl-config --cflags) -DGL_GLEXT_PROTOTYPES -CXML := $(shell xml2-config --cflags) -CTTF := $(shell freetype-config --cflags) -CPNG := #$(shell libpng-config libpng --cflags) -CGL := -CFLAGS := $(CLUA) $(CSDL) $(CXML) $(CTTF) $(CPNG) $(CGL) $(VERSION) -D$(OS) +CLUA := -I../lib/lua +CCSPARSE := -I../lib/csparse +CSDL := $(shell sdl-config --cflags) -DGL_GLEXT_PROTOTYPES +CXML := $(shell xml2-config --cflags) +CTTF := $(shell freetype-config --cflags) +CPNG := #$(shell libpng-config libpng --cflags) +CGL := +CFLAGS := $(CLUA) $(CCSPARSE) $(CSDL) $(CXML) $(CTTF) $(CPNG) $(CGL) $(VERSION) -D$(OS) ifdef DATA_DEF CFLAGS += -DDATA_DEF=$(DATA_DEF) endif @@ -75,13 +76,14 @@ CFLAGS += -O2 -funroll-loops -pipe -std=c99 endif # LDFLAGS. -LDLUA := ../lib/lua/liblua.a -LDSDL := $(shell sdl-config --libs) -lSDL_image -lSDL_mixer -lpng -LDXML := $(shell xml2-config --libs) -LDTTF := $(shell freetype-config --libs) -CPNG := #$(shell libpng-config libpng --libs) -LDGL := -lGL -LDFLAGS := -lm $(LDLUA) $(LDSDL) $(LDXML) $(LDTTF) $(CPNG) $(LDGL) $(LDPNG) +LDLUA := ../lib/lua/liblua.a +LDCSPARSE := ../lib/csparse/libcsparse.a +LDSDL := $(shell sdl-config --libs) -lSDL_image -lSDL_mixer -lpng +LDXML := $(shell xml2-config --libs) +LDTTF := $(shell freetype-config --libs) +CPNG := #$(shell libpng-config libpng --libs) +LDGL := -lGL +LDFLAGS := -lm $(LDLUA) $(LDCSPARSE) $(LDSDL) $(LDXML) $(LDTTF) $(CPNG) $(LDGL) $(LDPNG) # This is just for gstat to run some analysis on performance. ifdef DEBUG @@ -110,17 +112,18 @@ DATA := ldata DATAFILES := $(DATA_AI) $(DATA_GFX) $(DATA_XML) $(DATA_SND) $(DATA_MISN) # TARGETS. -.PHONY: all help lua utils docs clean distclean +.PHONY: all help csparse lua utils docs clean distclean %.o: %.c %.h @$(CC) -c $(CFLAGS) -o $@ $< @echo " CC $@" -all: utils ldata lua Lephisto +all: utils ldata lua csparse Lephisto help: @echo "Possible targets are:": @echo " lua - Builds Lua support." + @echo " csparse - Builds CSparse support." @echo " lephisto - Builds the Lephisto binary." @echo " mkspr - Builds the mkspr utility." @echo " ldata - Creates the ldata file." @@ -130,12 +133,18 @@ help: @echo " distclean - Removes everything done." $(APPNAME): $(OBJS) - @$(CC) $(LDFLAGS) -o $(APPNAME) $(OBJS) ../lib/lua/liblua.a + @$(CC) $(LDFLAGS) -o $(APPNAME) $(OBJS) ../lib/lua/liblua.a ../lib/csparse/libcsparse.a @echo " LD $(APPNAME)" -lua: - @if [ ! -e ../lib/lua/liblua.a ]; then $(MAKE) -C ../lib/lua a; fi +lua: ../lib/lua/liblua.a +../lib/lua/lublua.a: + +@$(MAKE) -C ../lib/lua a + +csparse: ../lib/csparse/libcsparse.a + +../lib/csparse/libcsparse.a: + +@$(MAKE) -C ../lib/csparse pack: ../src/md5.c ../src/pack.c ../utils/pack/main.c @$(MAKE) -C ../utils/pack @@ -169,7 +178,9 @@ distclean: clean @$(MAKE) -C ../utils/pack clean @$(MAKE) -C ../utils/mkspr clean @echo " Cleaning Lua" - @$(MAKE) -C ../lib/lua clean + @$(MAKE) -C ../lib/lua distclean + @echo " Cleaning CSparse." + @$(MAKE) -C ../lib/csparse disclean @echo " Removing build tool binaries." @$(RM) Lephisto mksprite ldata pack gmon.out VERSION diff --git a/lib/csparse/Makefile b/lib/csparse/Makefile new file mode 100644 index 0000000..d8c65c9 --- /dev/null +++ b/lib/csparse/Makefile @@ -0,0 +1,37 @@ +CFLAGS := -O3 + +AR := ar cr +RANLIB := ranlib + +LIB := libcsparse.a + +RELPATH := lib/csparse + +all: libcsparse.a + +%.o: %.c + @$(CC) -c $(CFLAGS) -o $@ $< + @echo " CC $(RELPATH)$@" + +CS = cs_add.o cs_amd.o cs_chol.o cs_cholsol.o cs_counts.o cs_cumsum.o \ + cs_droptol.o cs_dropzeros.o cs_dupl.o cs_entry.o \ + cs_etree.o cs_fkeep.o cs_gaxpy.o cs_happly.o cs_house.o cs_ipvec.o \ + cs_lsolve.o cs_ltsolve.o cs_lu.o cs_lusol.o cs_util.o cs_multiply.o \ + cs_permute.o cs_pinv.o cs_post.o cs_pvec.o cs_qr.o cs_qrsol.o \ + cs_scatter.o cs_schol.o cs_sqr.o cs_symperm.o cs_tdfs.o cs_malloc.o \ + cs_transpose.o cs_compress.o cs_usolve.o cs_utsolve.o cs_scc.o \ + cs_maxtrans.o cs_dmperm.o cs_updown.o cs_print.o cs_norm.o cs_load.o \ + cs_dfs.o cs_reach.o cs_spsolve.o cs_ereach.o cs_leaf.o cs_randperm.o + +libcsparse.a: $(CS) + @$(AR) $(LIB) $(CS) + @$(RANLIB) $(LIB) + @echo " LIB $(RELPATH)$@" + +clean: + @$(RM) *.o + @echo " Cleaning CSparse" + +distclean: clean + @$(RM) *.a + diff --git a/lib/csparse/cs.h b/lib/csparse/cs.h new file mode 100644 index 0000000..e7867b0 --- /dev/null +++ b/lib/csparse/cs.h @@ -0,0 +1,143 @@ +#ifndef _CS_H +#define _CS_H +#include +#include +#include +#include +#ifdef MATLAB_MEX_FILE +#include "mex.h" +#endif +#define CS_VER 2 /* CSparse Version 2.2.2 */ +#define CS_SUBVER 2 +#define CS_SUBSUB 2 +#define CS_DATE "Sept 23, 2008" /* CSparse release date */ +#define CS_COPYRIGHT "Copyright (c) Timothy A. Davis, 2006-2008" + +/* --- primary CSparse routines and data structures ------------------------- */ +typedef struct cs_sparse /* matrix in compressed-column or triplet form */ +{ + int nzmax ; /* maximum number of entries */ + int m ; /* number of rows */ + int n ; /* number of columns */ + int *p ; /* column pointers (size n+1) or col indices (size nzmax) */ + int *i ; /* row indices, size nzmax */ + double *x ; /* numerical values, size nzmax */ + int nz ; /* # of entries in triplet matrix, -1 for compressed-col */ +} cs ; + +cs *cs_add (const cs *A, const cs *B, double alpha, double beta) ; +int cs_cholsol (int order, const cs *A, double *b) ; +cs *cs_compress (const cs *T) ; +int cs_dupl (cs *A) ; +int cs_entry (cs *T, int i, int j, double x) ; +int cs_gaxpy (const cs *A, const double *x, double *y) ; +cs *cs_load (FILE *f) ; +int cs_lusol (int order, const cs *A, double *b, double tol) ; +cs *cs_multiply (const cs *A, const cs *B) ; +double cs_norm (const cs *A) ; +int cs_print (const cs *A, int brief) ; +int cs_qrsol (int order, const cs *A, double *b) ; +cs *cs_transpose (const cs *A, int values) ; +/* utilities */ +void *cs_calloc (int n, size_t size) ; +void *cs_free (void *p) ; +void *cs_realloc (void *p, int n, size_t size, int *ok) ; +cs *cs_spalloc (int m, int n, int nzmax, int values, int triplet) ; +cs *cs_spfree (cs *A) ; +int cs_sprealloc (cs *A, int nzmax) ; +void *cs_malloc (int n, size_t size) ; + +/* --- secondary CSparse routines and data structures ----------------------- */ +typedef struct cs_symbolic /* symbolic Cholesky, LU, or QR analysis */ +{ + int *pinv ; /* inverse row perm. for QR, fill red. perm for Chol */ + int *q ; /* fill-reducing column permutation for LU and QR */ + int *parent ; /* elimination tree for Cholesky and QR */ + int *cp ; /* column pointers for Cholesky, row counts for QR */ + int *leftmost ; /* leftmost[i] = min(find(A(i,:))), for QR */ + int m2 ; /* # of rows for QR, after adding fictitious rows */ + double lnz ; /* # entries in L for LU or Cholesky; in V for QR */ + double unz ; /* # entries in U for LU; in R for QR */ +} css ; + +typedef struct cs_numeric /* numeric Cholesky, LU, or QR factorization */ +{ + cs *L ; /* L for LU and Cholesky, V for QR */ + cs *U ; /* U for LU, R for QR, not used for Cholesky */ + int *pinv ; /* partial pivoting for LU */ + double *B ; /* beta [0..n-1] for QR */ +} csn ; + +typedef struct cs_dmperm_results /* cs_dmperm or cs_scc output */ +{ + int *p ; /* size m, row permutation */ + int *q ; /* size n, column permutation */ + int *r ; /* size nb+1, block k is rows r[k] to r[k+1]-1 in A(p,q) */ + int *s ; /* size nb+1, block k is cols s[k] to s[k+1]-1 in A(p,q) */ + int nb ; /* # of blocks in fine dmperm decomposition */ + int rr [5] ; /* coarse row decomposition */ + int cc [5] ; /* coarse column decomposition */ +} csd ; + +int *cs_amd (int order, const cs *A) ; +csn *cs_chol (const cs *A, const css *S) ; +csd *cs_dmperm (const cs *A, int seed) ; +int cs_droptol (cs *A, double tol) ; +int cs_dropzeros (cs *A) ; +int cs_happly (const cs *V, int i, double beta, double *x) ; +int cs_ipvec (const int *p, const double *b, double *x, int n) ; +int cs_lsolve (const cs *L, double *x) ; +int cs_ltsolve (const cs *L, double *x) ; +csn *cs_lu (const cs *A, const css *S, double tol) ; +cs *cs_permute (const cs *A, const int *pinv, const int *q, int values) ; +int *cs_pinv (const int *p, int n) ; +int cs_pvec (const int *p, const double *b, double *x, int n) ; +csn *cs_qr (const cs *A, const css *S) ; +css *cs_schol (int order, const cs *A) ; +css *cs_sqr (int order, const cs *A, int qr) ; +cs *cs_symperm (const cs *A, const int *pinv, int values) ; +int cs_updown (cs *L, int sigma, const cs *C, const int *parent) ; +int cs_usolve (const cs *U, double *x) ; +int cs_utsolve (const cs *U, double *x) ; +/* utilities */ +css *cs_sfree (css *S) ; +csn *cs_nfree (csn *N) ; +csd *cs_dfree (csd *D) ; + +/* --- tertiary CSparse routines -------------------------------------------- */ +int *cs_counts (const cs *A, const int *parent, const int *post, int ata) ; +double cs_cumsum (int *p, int *c, int n) ; +int cs_dfs (int j, cs *G, int top, int *xi, int *pstack, const int *pinv) ; +int cs_ereach (const cs *A, int k, const int *parent, int *s, int *w) ; +int *cs_etree (const cs *A, int ata) ; +int cs_fkeep (cs *A, int (*fkeep) (int, int, double, void *), void *other) ; +double cs_house (double *x, double *beta, int n) ; +int cs_leaf (int i, int j, const int *first, int *maxfirst, int *prevleaf, + int *ancestor, int *jleaf) ; +int *cs_maxtrans (const cs *A, int seed) ; +int *cs_post (const int *parent, int n) ; +int *cs_randperm (int n, int seed) ; +int cs_reach (cs *G, const cs *B, int k, int *xi, const int *pinv) ; +int cs_scatter (const cs *A, int j, double beta, int *w, double *x, int mark, + cs *C, int nz) ; +csd *cs_scc (cs *A) ; +int cs_spsolve (cs *G, const cs *B, int k, int *xi, double *x, + const int *pinv, int lo) ; +int cs_tdfs (int j, int k, int *head, const int *next, int *post, + int *stack) ; +/* utilities */ +csd *cs_dalloc (int m, int n) ; +csd *cs_ddone (csd *D, cs *C, void *w, int ok) ; +cs *cs_done (cs *C, void *w, void *x, int ok) ; +int *cs_idone (int *p, cs *C, void *w, int ok) ; +csn *cs_ndone (csn *N, cs *C, void *w, void *x, int ok) ; + +#define CS_MAX(a,b) (((a) > (b)) ? (a) : (b)) +#define CS_MIN(a,b) (((a) < (b)) ? (a) : (b)) +#define CS_FLIP(i) (-(i)-2) +#define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i)) +#define CS_MARKED(w,j) (w [j] < 0) +#define CS_MARK(w,j) { w [j] = CS_FLIP (w [j]) ; } +#define CS_CSC(A) (A && (A->nz == -1)) +#define CS_TRIPLET(A) (A && (A->nz >= 0)) +#endif diff --git a/lib/csparse/cs_add.c b/lib/csparse/cs_add.c new file mode 100644 index 0000000..2c3c146 --- /dev/null +++ b/lib/csparse/cs_add.c @@ -0,0 +1,28 @@ +#include "cs.h" +/* C = alpha*A + beta*B */ +cs *cs_add (const cs *A, const cs *B, double alpha, double beta) +{ + int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values ; + double *x, *Bx, *Cx ; + cs *C ; + if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ; /* check inputs */ + if (A->m != B->m || A->n != B->n) return (NULL) ; + m = A->m ; anz = A->p [A->n] ; + n = B->n ; Bp = B->p ; Bx = B->x ; bnz = Bp [n] ; + w = cs_calloc (m, sizeof (int)) ; /* get workspace */ + values = (A->x != NULL) && (Bx != NULL) ; + x = values ? cs_malloc (m, sizeof (double)) : NULL ; /* get workspace */ + C = cs_spalloc (m, n, anz + bnz, values, 0) ; /* allocate result*/ + if (!C || !w || (values && !x)) return (cs_done (C, w, x, 0)) ; + Cp = C->p ; Ci = C->i ; Cx = C->x ; + for (j = 0 ; j < n ; j++) + { + Cp [j] = nz ; /* column j of C starts here */ + nz = cs_scatter (A, j, alpha, w, x, j+1, C, nz) ; /* alpha*A(:,j)*/ + nz = cs_scatter (B, j, beta, w, x, j+1, C, nz) ; /* beta*B(:,j) */ + if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Ci [p]] ; + } + Cp [n] = nz ; /* finalize the last column of C */ + cs_sprealloc (C, 0) ; /* remove extra space from C */ + return (cs_done (C, w, x, 1)) ; /* success; free workspace, return C */ +} diff --git a/lib/csparse/cs_amd.c b/lib/csparse/cs_amd.c new file mode 100644 index 0000000..58b5d0e --- /dev/null +++ b/lib/csparse/cs_amd.c @@ -0,0 +1,364 @@ +#include "cs.h" +/* clear w */ +static int cs_wclear (int mark, int lemax, int *w, int n) +{ + int k ; + if (mark < 2 || (mark + lemax < 0)) + { + for (k = 0 ; k < n ; k++) if (w [k] != 0) w [k] = 1 ; + mark = 2 ; + } + return (mark) ; /* at this point, w [0..n-1] < mark holds */ +} + +/* keep off-diagonal entries; drop diagonal entries */ +static int cs_diag (int i, int j, double aij, void *other) { return (i != j) ; } + +/* p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */ +int *cs_amd (int order, const cs *A) /* order 0:natural, 1:Chol, 2:LU, 3:QR */ +{ + cs *C, *A2, *AT ; + int *Cp, *Ci, *last, *W, *len, *nv, *next, *P, *head, *elen, *degree, *w, + *hhead, *ATp, *ATi, d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1, + k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi, + ok, cnz, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, n, m, t ; + unsigned int h ; + /* --- Construct matrix C ----------------------------------------------- */ + if (!CS_CSC (A) || order <= 0 || order > 3) return (NULL) ; /* check */ + AT = cs_transpose (A, 0) ; /* compute A' */ + if (!AT) return (NULL) ; + m = A->m ; n = A->n ; + dense = CS_MAX (16, 10 * sqrt ((double) n)) ; /* find dense threshold */ + dense = CS_MIN (n-2, dense) ; + if (order == 1 && n == m) + { + C = cs_add (A, AT, 0, 0) ; /* C = A+A' */ + } + else if (order == 2) + { + ATp = AT->p ; /* drop dense columns from AT */ + ATi = AT->i ; + for (p2 = 0, j = 0 ; j < m ; j++) + { + p = ATp [j] ; /* column j of AT starts here */ + ATp [j] = p2 ; /* new column j starts here */ + if (ATp [j+1] - p > dense) continue ; /* skip dense col j */ + for ( ; p < ATp [j+1] ; p++) ATi [p2++] = ATi [p] ; + } + ATp [m] = p2 ; /* finalize AT */ + A2 = cs_transpose (AT, 0) ; /* A2 = AT' */ + C = A2 ? cs_multiply (AT, A2) : NULL ; /* C=A'*A with no dense rows */ + cs_spfree (A2) ; + } + else + { + C = cs_multiply (AT, A) ; /* C=A'*A */ + } + cs_spfree (AT) ; + if (!C) return (NULL) ; + cs_fkeep (C, &cs_diag, NULL) ; /* drop diagonal entries */ + Cp = C->p ; + cnz = Cp [n] ; + P = cs_malloc (n+1, sizeof (int)) ; /* allocate result */ + W = cs_malloc (8*(n+1), sizeof (int)) ; /* get workspace */ + t = cnz + cnz/5 + 2*n ; /* add elbow room to C */ + if (!P || !W || !cs_sprealloc (C, t)) return (cs_idone (P, C, W, 0)) ; + len = W ; nv = W + (n+1) ; next = W + 2*(n+1) ; + head = W + 3*(n+1) ; elen = W + 4*(n+1) ; degree = W + 5*(n+1) ; + w = W + 6*(n+1) ; hhead = W + 7*(n+1) ; + last = P ; /* use P as workspace for last */ + /* --- Initialize quotient graph ---------------------------------------- */ + for (k = 0 ; k < n ; k++) len [k] = Cp [k+1] - Cp [k] ; + len [n] = 0 ; + nzmax = C->nzmax ; + Ci = C->i ; + for (i = 0 ; i <= n ; i++) + { + head [i] = -1 ; /* degree list i is empty */ + last [i] = -1 ; + next [i] = -1 ; + hhead [i] = -1 ; /* hash list i is empty */ + nv [i] = 1 ; /* node i is just one node */ + w [i] = 1 ; /* node i is alive */ + elen [i] = 0 ; /* Ek of node i is empty */ + degree [i] = len [i] ; /* degree of node i */ + } + mark = cs_wclear (0, 0, w, n) ; /* clear w */ + elen [n] = -2 ; /* n is a dead element */ + Cp [n] = -1 ; /* n is a root of assembly tree */ + w [n] = 0 ; /* n is a dead element */ + /* --- Initialize degree lists ------------------------------------------ */ + for (i = 0 ; i < n ; i++) + { + d = degree [i] ; + if (d == 0) /* node i is empty */ + { + elen [i] = -2 ; /* element i is dead */ + nel++ ; + Cp [i] = -1 ; /* i is a root of assembly tree */ + w [i] = 0 ; + } + else if (d > dense) /* node i is dense */ + { + nv [i] = 0 ; /* absorb i into element n */ + elen [i] = -1 ; /* node i is dead */ + nel++ ; + Cp [i] = CS_FLIP (n) ; + nv [n]++ ; + } + else + { + if (head [d] != -1) last [head [d]] = i ; + next [i] = head [d] ; /* put node i in degree list d */ + head [d] = i ; + } + } + while (nel < n) /* while (selecting pivots) do */ + { + /* --- Select node of minimum approximate degree -------------------- */ + for (k = -1 ; mindeg < n && (k = head [mindeg]) == -1 ; mindeg++) ; + if (next [k] != -1) last [next [k]] = -1 ; + head [mindeg] = next [k] ; /* remove k from degree list */ + elenk = elen [k] ; /* elenk = |Ek| */ + nvk = nv [k] ; /* # of nodes k represents */ + nel += nvk ; /* nv[k] nodes of A eliminated */ + /* --- Garbage collection ------------------------------------------- */ + if (elenk > 0 && cnz + mindeg >= nzmax) + { + for (j = 0 ; j < n ; j++) + { + if ((p = Cp [j]) >= 0) /* j is a live node or element */ + { + Cp [j] = Ci [p] ; /* save first entry of object */ + Ci [p] = CS_FLIP (j) ; /* first entry is now CS_FLIP(j) */ + } + } + for (q = 0, p = 0 ; p < cnz ; ) /* scan all of memory */ + { + if ((j = CS_FLIP (Ci [p++])) >= 0) /* found object j */ + { + Ci [q] = Cp [j] ; /* restore first entry of object */ + Cp [j] = q++ ; /* new pointer to object j */ + for (k3 = 0 ; k3 < len [j]-1 ; k3++) Ci [q++] = Ci [p++] ; + } + } + cnz = q ; /* Ci [cnz...nzmax-1] now free */ + } + /* --- Construct new element ---------------------------------------- */ + dk = 0 ; + nv [k] = -nvk ; /* flag k as in Lk */ + p = Cp [k] ; + pk1 = (elenk == 0) ? p : cnz ; /* do in place if elen[k] == 0 */ + pk2 = pk1 ; + for (k1 = 1 ; k1 <= elenk + 1 ; k1++) + { + if (k1 > elenk) + { + e = k ; /* search the nodes in k */ + pj = p ; /* list of nodes starts at Ci[pj]*/ + ln = len [k] - elenk ; /* length of list of nodes in k */ + } + else + { + e = Ci [p++] ; /* search the nodes in e */ + pj = Cp [e] ; + ln = len [e] ; /* length of list of nodes in e */ + } + for (k2 = 1 ; k2 <= ln ; k2++) + { + i = Ci [pj++] ; + if ((nvi = nv [i]) <= 0) continue ; /* node i dead, or seen */ + dk += nvi ; /* degree[Lk] += size of node i */ + nv [i] = -nvi ; /* negate nv[i] to denote i in Lk*/ + Ci [pk2++] = i ; /* place i in Lk */ + if (next [i] != -1) last [next [i]] = last [i] ; + if (last [i] != -1) /* remove i from degree list */ + { + next [last [i]] = next [i] ; + } + else + { + head [degree [i]] = next [i] ; + } + } + if (e != k) + { + Cp [e] = CS_FLIP (k) ; /* absorb e into k */ + w [e] = 0 ; /* e is now a dead element */ + } + } + if (elenk != 0) cnz = pk2 ; /* Ci [cnz...nzmax] is free */ + degree [k] = dk ; /* external degree of k - |Lk\i| */ + Cp [k] = pk1 ; /* element k is in Ci[pk1..pk2-1] */ + len [k] = pk2 - pk1 ; + elen [k] = -2 ; /* k is now an element */ + /* --- Find set differences ----------------------------------------- */ + mark = cs_wclear (mark, lemax, w, n) ; /* clear w if necessary */ + for (pk = pk1 ; pk < pk2 ; pk++) /* scan 1: find |Le\Lk| */ + { + i = Ci [pk] ; + if ((eln = elen [i]) <= 0) continue ;/* skip if elen[i] empty */ + nvi = -nv [i] ; /* nv [i] was negated */ + wnvi = mark - nvi ; + for (p = Cp [i] ; p <= Cp [i] + eln - 1 ; p++) /* scan Ei */ + { + e = Ci [p] ; + if (w [e] >= mark) + { + w [e] -= nvi ; /* decrement |Le\Lk| */ + } + else if (w [e] != 0) /* ensure e is a live element */ + { + w [e] = degree [e] + wnvi ; /* 1st time e seen in scan 1 */ + } + } + } + /* --- Degree update ------------------------------------------------ */ + for (pk = pk1 ; pk < pk2 ; pk++) /* scan2: degree update */ + { + i = Ci [pk] ; /* consider node i in Lk */ + p1 = Cp [i] ; + p2 = p1 + elen [i] - 1 ; + pn = p1 ; + for (h = 0, d = 0, p = p1 ; p <= p2 ; p++) /* scan Ei */ + { + e = Ci [p] ; + if (w [e] != 0) /* e is an unabsorbed element */ + { + dext = w [e] - mark ; /* dext = |Le\Lk| */ + if (dext > 0) + { + d += dext ; /* sum up the set differences */ + Ci [pn++] = e ; /* keep e in Ei */ + h += e ; /* compute the hash of node i */ + } + else + { + Cp [e] = CS_FLIP (k) ; /* aggressive absorb. e->k */ + w [e] = 0 ; /* e is a dead element */ + } + } + } + elen [i] = pn - p1 + 1 ; /* elen[i] = |Ei| */ + p3 = pn ; + p4 = p1 + len [i] ; + for (p = p2 + 1 ; p < p4 ; p++) /* prune edges in Ai */ + { + j = Ci [p] ; + if ((nvj = nv [j]) <= 0) continue ; /* node j dead or in Lk */ + d += nvj ; /* degree(i) += |j| */ + Ci [pn++] = j ; /* place j in node list of i */ + h += j ; /* compute hash for node i */ + } + if (d == 0) /* check for mass elimination */ + { + Cp [i] = CS_FLIP (k) ; /* absorb i into k */ + nvi = -nv [i] ; + dk -= nvi ; /* |Lk| -= |i| */ + nvk += nvi ; /* |k| += nv[i] */ + nel += nvi ; + nv [i] = 0 ; + elen [i] = -1 ; /* node i is dead */ + } + else + { + degree [i] = CS_MIN (degree [i], d) ; /* update degree(i) */ + Ci [pn] = Ci [p3] ; /* move first node to end */ + Ci [p3] = Ci [p1] ; /* move 1st el. to end of Ei */ + Ci [p1] = k ; /* add k as 1st element in of Ei */ + len [i] = pn - p1 + 1 ; /* new len of adj. list of node i */ + h %= n ; /* finalize hash of i */ + next [i] = hhead [h] ; /* place i in hash bucket */ + hhead [h] = i ; + last [i] = h ; /* save hash of i in last[i] */ + } + } /* scan2 is done */ + degree [k] = dk ; /* finalize |Lk| */ + lemax = CS_MAX (lemax, dk) ; + mark = cs_wclear (mark+lemax, lemax, w, n) ; /* clear w */ + /* --- Supernode detection ------------------------------------------ */ + for (pk = pk1 ; pk < pk2 ; pk++) + { + i = Ci [pk] ; + if (nv [i] >= 0) continue ; /* skip if i is dead */ + h = last [i] ; /* scan hash bucket of node i */ + i = hhead [h] ; + hhead [h] = -1 ; /* hash bucket will be empty */ + for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++) + { + ln = len [i] ; + eln = elen [i] ; + for (p = Cp [i]+1 ; p <= Cp [i] + ln-1 ; p++) w [Ci [p]] = mark; + jlast = i ; + for (j = next [i] ; j != -1 ; ) /* compare i with all j */ + { + ok = (len [j] == ln) && (elen [j] == eln) ; + for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++) + { + if (w [Ci [p]] != mark) ok = 0 ; /* compare i and j*/ + } + if (ok) /* i and j are identical */ + { + Cp [j] = CS_FLIP (i) ; /* absorb j into i */ + nv [i] += nv [j] ; + nv [j] = 0 ; + elen [j] = -1 ; /* node j is dead */ + j = next [j] ; /* delete j from hash bucket */ + next [jlast] = j ; + } + else + { + jlast = j ; /* j and i are different */ + j = next [j] ; + } + } + } + } + /* --- Finalize new element------------------------------------------ */ + for (p = pk1, pk = pk1 ; pk < pk2 ; pk++) /* finalize Lk */ + { + i = Ci [pk] ; + if ((nvi = -nv [i]) <= 0) continue ;/* skip if i is dead */ + nv [i] = nvi ; /* restore nv[i] */ + d = degree [i] + dk - nvi ; /* compute external degree(i) */ + d = CS_MIN (d, n - nel - nvi) ; + if (head [d] != -1) last [head [d]] = i ; + next [i] = head [d] ; /* put i back in degree list */ + last [i] = -1 ; + head [d] = i ; + mindeg = CS_MIN (mindeg, d) ; /* find new minimum degree */ + degree [i] = d ; + Ci [p++] = i ; /* place i in Lk */ + } + nv [k] = nvk ; /* # nodes absorbed into k */ + if ((len [k] = p-pk1) == 0) /* length of adj list of element k*/ + { + Cp [k] = -1 ; /* k is a root of the tree */ + w [k] = 0 ; /* k is now a dead element */ + } + if (elenk != 0) cnz = p ; /* free unused space in Lk */ + } + /* --- Postordering ----------------------------------------------------- */ + for (i = 0 ; i < n ; i++) Cp [i] = CS_FLIP (Cp [i]) ;/* fix assembly tree */ + for (j = 0 ; j <= n ; j++) head [j] = -1 ; + for (j = n ; j >= 0 ; j--) /* place unordered nodes in lists */ + { + if (nv [j] > 0) continue ; /* skip if j is an element */ + next [j] = head [Cp [j]] ; /* place j in list of its parent */ + head [Cp [j]] = j ; + } + for (e = n ; e >= 0 ; e--) /* place elements in lists */ + { + if (nv [e] <= 0) continue ; /* skip unless e is an element */ + if (Cp [e] != -1) + { + next [e] = head [Cp [e]] ; /* place e in list of its parent */ + head [Cp [e]] = e ; + } + } + for (k = 0, i = 0 ; i <= n ; i++) /* postorder the assembly tree */ + { + if (Cp [i] == -1) k = cs_tdfs (i, k, head, next, P, w) ; + } + return (cs_idone (P, C, W, 1)) ; +} diff --git a/lib/csparse/cs_chol.c b/lib/csparse/cs_chol.c new file mode 100644 index 0000000..5eb12b3 --- /dev/null +++ b/lib/csparse/cs_chol.c @@ -0,0 +1,58 @@ +#include "cs.h" +/* L = chol (A, [pinv parent cp]), pinv is optional */ +csn *cs_chol (const cs *A, const css *S) +{ + double d, lki, *Lx, *x, *Cx ; + int top, i, p, k, n, *Li, *Lp, *cp, *pinv, *s, *c, *parent, *Cp, *Ci ; + cs *L, *C, *E ; + csn *N ; + if (!CS_CSC (A) || !S || !S->cp || !S->parent) return (NULL) ; + n = A->n ; + N = cs_calloc (1, sizeof (csn)) ; /* allocate result */ + c = cs_malloc (2*n, sizeof (int)) ; /* get int workspace */ + x = cs_malloc (n, sizeof (double)) ; /* get double workspace */ + cp = S->cp ; pinv = S->pinv ; parent = S->parent ; + C = pinv ? cs_symperm (A, pinv, 1) : ((cs *) A) ; + E = pinv ? C : NULL ; /* E is alias for A, or a copy E=A(p,p) */ + if (!N || !c || !x || !C) return (cs_ndone (N, E, c, x, 0)) ; + s = c + n ; + Cp = C->p ; Ci = C->i ; Cx = C->x ; + N->L = L = cs_spalloc (n, n, cp [n], 1, 0) ; /* allocate result */ + if (!L) return (cs_ndone (N, E, c, x, 0)) ; + Lp = L->p ; Li = L->i ; Lx = L->x ; + for (k = 0 ; k < n ; k++) Lp [k] = c [k] = cp [k] ; + for (k = 0 ; k < n ; k++) /* compute L(:,k) for L*L' = C */ + { + /* --- Nonzero pattern of L(k,:) ------------------------------------ */ + top = cs_ereach (C, k, parent, s, c) ; /* find pattern of L(k,:) */ + x [k] = 0 ; /* x (0:k) is now zero */ + for (p = Cp [k] ; p < Cp [k+1] ; p++) /* x = full(triu(C(:,k))) */ + { + if (Ci [p] <= k) x [Ci [p]] = Cx [p] ; + } + d = x [k] ; /* d = C(k,k) */ + x [k] = 0 ; /* clear x for k+1st iteration */ + /* --- Triangular solve --------------------------------------------- */ + for ( ; top < n ; top++) /* solve L(0:k-1,0:k-1) * x = C(:,k) */ + { + i = s [top] ; /* s [top..n-1] is pattern of L(k,:) */ + lki = x [i] / Lx [Lp [i]] ; /* L(k,i) = x (i) / L(i,i) */ + x [i] = 0 ; /* clear x for k+1st iteration */ + for (p = Lp [i] + 1 ; p < c [i] ; p++) + { + x [Li [p]] -= Lx [p] * lki ; + } + d -= lki * lki ; /* d = d - L(k,i)*L(k,i) */ + p = c [i]++ ; + Li [p] = k ; /* store L(k,i) in column i */ + Lx [p] = lki ; + } + /* --- Compute L(k,k) ----------------------------------------------- */ + if (d <= 0) return (cs_ndone (N, E, c, x, 0)) ; /* not pos def */ + p = c [k]++ ; + Li [p] = k ; /* store L(k,k) = sqrt (d) in column k */ + Lx [p] = sqrt (d) ; + } + Lp [n] = cp [n] ; /* finalize L */ + return (cs_ndone (N, E, c, x, 1)) ; /* success: free E,s,x; return N */ +} diff --git a/lib/csparse/cs_cholsol.c b/lib/csparse/cs_cholsol.c new file mode 100644 index 0000000..f3cf414 --- /dev/null +++ b/lib/csparse/cs_cholsol.c @@ -0,0 +1,26 @@ +#include "cs.h" +/* x=A\b where A is symmetric positive definite; b overwritten with solution */ +int cs_cholsol (int order, const cs *A, double *b) +{ + double *x ; + css *S ; + csn *N ; + int n, ok ; + if (!CS_CSC (A) || !b) return (0) ; /* check inputs */ + n = A->n ; + S = cs_schol (order, A) ; /* ordering and symbolic analysis */ + N = cs_chol (A, S) ; /* numeric Cholesky factorization */ + x = cs_malloc (n, sizeof (double)) ; /* get workspace */ + ok = (S && N && x) ; + if (ok) + { + cs_ipvec (S->pinv, b, x, n) ; /* x = P*b */ + cs_lsolve (N->L, x) ; /* x = L\x */ + cs_ltsolve (N->L, x) ; /* x = L'\x */ + cs_pvec (S->pinv, x, b, n) ; /* b = P'*x */ + } + cs_free (x) ; + cs_sfree (S) ; + cs_nfree (N) ; + return (ok) ; +} diff --git a/lib/csparse/cs_compress.c b/lib/csparse/cs_compress.c new file mode 100644 index 0000000..d18e840 --- /dev/null +++ b/lib/csparse/cs_compress.c @@ -0,0 +1,22 @@ +#include "cs.h" +/* C = compressed-column form of a triplet matrix T */ +cs *cs_compress (const cs *T) +{ + int m, n, nz, p, k, *Cp, *Ci, *w, *Ti, *Tj ; + double *Cx, *Tx ; + cs *C ; + if (!CS_TRIPLET (T)) return (NULL) ; /* check inputs */ + m = T->m ; n = T->n ; Ti = T->i ; Tj = T->p ; Tx = T->x ; nz = T->nz ; + C = cs_spalloc (m, n, nz, Tx != NULL, 0) ; /* allocate result */ + w = cs_calloc (n, sizeof (int)) ; /* get workspace */ + if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */ + Cp = C->p ; Ci = C->i ; Cx = C->x ; + for (k = 0 ; k < nz ; k++) w [Tj [k]]++ ; /* column counts */ + cs_cumsum (Cp, w, n) ; /* column pointers */ + for (k = 0 ; k < nz ; k++) + { + Ci [p = w [Tj [k]]++] = Ti [k] ; /* A(i,j) is the pth entry in C */ + if (Cx) Cx [p] = Tx [k] ; + } + return (cs_done (C, w, NULL, 1)) ; /* success; free w and return C */ +} diff --git a/lib/csparse/cs_counts.c b/lib/csparse/cs_counts.c new file mode 100644 index 0000000..ff32b6f --- /dev/null +++ b/lib/csparse/cs_counts.c @@ -0,0 +1,61 @@ +#include "cs.h" +/* column counts of LL'=A or LL'=A'A, given parent & post ordering */ +#define HEAD(k,j) (ata ? head [k] : j) +#define NEXT(J) (ata ? next [J] : -1) +static void init_ata (cs *AT, const int *post, int *w, int **head, int **next) +{ + int i, k, p, m = AT->n, n = AT->m, *ATp = AT->p, *ATi = AT->i ; + *head = w+4*n, *next = w+5*n+1 ; + for (k = 0 ; k < n ; k++) w [post [k]] = k ; /* invert post */ + for (i = 0 ; i < m ; i++) + { + for (k = n, p = ATp[i] ; p < ATp[i+1] ; p++) k = CS_MIN (k, w [ATi[p]]); + (*next) [i] = (*head) [k] ; /* place row i in linked list k */ + (*head) [k] = i ; + } +} +int *cs_counts (const cs *A, const int *parent, const int *post, int ata) +{ + int i, j, k, n, m, J, s, p, q, jleaf, *ATp, *ATi, *maxfirst, *prevleaf, + *ancestor, *head = NULL, *next = NULL, *colcount, *w, *first, *delta ; + cs *AT ; + if (!CS_CSC (A) || !parent || !post) return (NULL) ; /* check inputs */ + m = A->m ; n = A->n ; + s = 4*n + (ata ? (n+m+1) : 0) ; + delta = colcount = cs_malloc (n, sizeof (int)) ; /* allocate result */ + w = cs_malloc (s, sizeof (int)) ; /* get workspace */ + AT = cs_transpose (A, 0) ; /* AT = A' */ + if (!AT || !colcount || !w) return (cs_idone (colcount, AT, w, 0)) ; + ancestor = w ; maxfirst = w+n ; prevleaf = w+2*n ; first = w+3*n ; + for (k = 0 ; k < s ; k++) w [k] = -1 ; /* clear workspace w [0..s-1] */ + for (k = 0 ; k < n ; k++) /* find first [j] */ + { + j = post [k] ; + delta [j] = (first [j] == -1) ? 1 : 0 ; /* delta[j]=1 if j is a leaf */ + for ( ; j != -1 && first [j] == -1 ; j = parent [j]) first [j] = k ; + } + ATp = AT->p ; ATi = AT->i ; + if (ata) init_ata (AT, post, w, &head, &next) ; + for (i = 0 ; i < n ; i++) ancestor [i] = i ; /* each node in its own set */ + for (k = 0 ; k < n ; k++) + { + j = post [k] ; /* j is the kth node in postordered etree */ + if (parent [j] != -1) delta [parent [j]]-- ; /* j is not a root */ + for (J = HEAD (k,j) ; J != -1 ; J = NEXT (J)) /* J=j for LL'=A case */ + { + for (p = ATp [J] ; p < ATp [J+1] ; p++) + { + i = ATi [p] ; + q = cs_leaf (i, j, first, maxfirst, prevleaf, ancestor, &jleaf); + if (jleaf >= 1) delta [j]++ ; /* A(i,j) is in skeleton */ + if (jleaf == 2) delta [q]-- ; /* account for overlap in q */ + } + } + if (parent [j] != -1) ancestor [j] = parent [j] ; + } + for (j = 0 ; j < n ; j++) /* sum up delta's of each child */ + { + if (parent [j] != -1) colcount [parent [j]] += colcount [j] ; + } + return (cs_idone (colcount, AT, w, 1)) ; /* success: free workspace */ +} diff --git a/lib/csparse/cs_cumsum.c b/lib/csparse/cs_cumsum.c new file mode 100644 index 0000000..b8b7778 --- /dev/null +++ b/lib/csparse/cs_cumsum.c @@ -0,0 +1,17 @@ +#include "cs.h" +/* p [0..n] = cumulative sum of c [0..n-1], and then copy p [0..n-1] into c */ +double cs_cumsum (int *p, int *c, int n) +{ + int i, nz = 0 ; + double nz2 = 0 ; + if (!p || !c) return (-1) ; /* check inputs */ + for (i = 0 ; i < n ; i++) + { + p [i] = nz ; + nz += c [i] ; + nz2 += c [i] ; /* also in double to avoid int overflow */ + c [i] = p [i] ; /* also copy p[0..n-1] back into c[0..n-1]*/ + } + p [n] = nz ; + return (nz2) ; /* return sum (c [0..n-1]) */ +} diff --git a/lib/csparse/cs_dfs.c b/lib/csparse/cs_dfs.c new file mode 100644 index 0000000..250dcf3 --- /dev/null +++ b/lib/csparse/cs_dfs.c @@ -0,0 +1,36 @@ +#include "cs.h" +/* depth-first-search of the graph of a matrix, starting at node j */ +int cs_dfs (int j, cs *G, int top, int *xi, int *pstack, const int *pinv) +{ + int i, p, p2, done, jnew, head = 0, *Gp, *Gi ; + if (!CS_CSC (G) || !xi || !pstack) return (-1) ; /* check inputs */ + Gp = G->p ; Gi = G->i ; + xi [0] = j ; /* initialize the recursion stack */ + while (head >= 0) + { + j = xi [head] ; /* get j from the top of the recursion stack */ + jnew = pinv ? (pinv [j]) : j ; + if (!CS_MARKED (Gp, j)) + { + CS_MARK (Gp, j) ; /* mark node j as visited */ + pstack [head] = (jnew < 0) ? 0 : CS_UNFLIP (Gp [jnew]) ; + } + done = 1 ; /* node j done if no unvisited neighbors */ + p2 = (jnew < 0) ? 0 : CS_UNFLIP (Gp [jnew+1]) ; + for (p = pstack [head] ; p < p2 ; p++) /* examine all neighbors of j */ + { + i = Gi [p] ; /* consider neighbor node i */ + if (CS_MARKED (Gp, i)) continue ; /* skip visited node i */ + pstack [head] = p ; /* pause depth-first search of node j */ + xi [++head] = i ; /* start dfs at node i */ + done = 0 ; /* node j is not done */ + break ; /* break, to start dfs (i) */ + } + if (done) /* depth-first search at node j is done */ + { + head-- ; /* remove j from the recursion stack */ + xi [--top] = j ; /* and place in the output stack */ + } + } + return (top) ; +} diff --git a/lib/csparse/cs_dmperm.c b/lib/csparse/cs_dmperm.c new file mode 100644 index 0000000..29f72f6 --- /dev/null +++ b/lib/csparse/cs_dmperm.c @@ -0,0 +1,144 @@ +#include "cs.h" +/* breadth-first search for coarse decomposition (C0,C1,R1 or R0,R3,C3) */ +static int cs_bfs (const cs *A, int n, int *wi, int *wj, int *queue, + const int *imatch, const int *jmatch, int mark) +{ + int *Ap, *Ai, head = 0, tail = 0, j, i, p, j2 ; + cs *C ; + for (j = 0 ; j < n ; j++) /* place all unmatched nodes in queue */ + { + if (imatch [j] >= 0) continue ; /* skip j if matched */ + wj [j] = 0 ; /* j in set C0 (R0 if transpose) */ + queue [tail++] = j ; /* place unmatched col j in queue */ + } + if (tail == 0) return (1) ; /* quick return if no unmatched nodes */ + C = (mark == 1) ? ((cs *) A) : cs_transpose (A, 0) ; + if (!C) return (0) ; /* bfs of C=A' to find R3,C3 from R0 */ + Ap = C->p ; Ai = C->i ; + while (head < tail) /* while queue is not empty */ + { + j = queue [head++] ; /* get the head of the queue */ + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + i = Ai [p] ; + if (wi [i] >= 0) continue ; /* skip if i is marked */ + wi [i] = mark ; /* i in set R1 (C3 if transpose) */ + j2 = jmatch [i] ; /* traverse alternating path to j2 */ + if (wj [j2] >= 0) continue ;/* skip j2 if it is marked */ + wj [j2] = mark ; /* j2 in set C1 (R3 if transpose) */ + queue [tail++] = j2 ; /* add j2 to queue */ + } + } + if (mark != 1) cs_spfree (C) ; /* free A' if it was created */ + return (1) ; +} + +/* collect matched rows and columns into p and q */ +static void cs_matched (int n, const int *wj, const int *imatch, int *p, int *q, + int *cc, int *rr, int set, int mark) +{ + int kc = cc [set], j ; + int kr = rr [set-1] ; + for (j = 0 ; j < n ; j++) + { + if (wj [j] != mark) continue ; /* skip if j is not in C set */ + p [kr++] = imatch [j] ; + q [kc++] = j ; + } + cc [set+1] = kc ; + rr [set] = kr ; +} + +/* collect unmatched rows into the permutation vector p */ +static void cs_unmatched (int m, const int *wi, int *p, int *rr, int set) +{ + int i, kr = rr [set] ; + for (i = 0 ; i < m ; i++) if (wi [i] == 0) p [kr++] = i ; + rr [set+1] = kr ; +} + +/* return 1 if row i is in R2 */ +static int cs_rprune (int i, int j, double aij, void *other) +{ + int *rr = (int *) other ; + return (i >= rr [1] && i < rr [2]) ; +} + +/* Given A, compute coarse and then fine dmperm */ +csd *cs_dmperm (const cs *A, int seed) +{ + int m, n, i, j, k, cnz, nc, *jmatch, *imatch, *wi, *wj, *pinv, *Cp, *Ci, + *ps, *rs, nb1, nb2, *p, *q, *cc, *rr, *r, *s, ok ; + cs *C ; + csd *D, *scc ; + /* --- Maximum matching ------------------------------------------------- */ + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + m = A->m ; n = A->n ; + D = cs_dalloc (m, n) ; /* allocate result */ + if (!D) return (NULL) ; + p = D->p ; q = D->q ; r = D->r ; s = D->s ; cc = D->cc ; rr = D->rr ; + jmatch = cs_maxtrans (A, seed) ; /* max transversal */ + imatch = jmatch + m ; /* imatch = inverse of jmatch */ + if (!jmatch) return (cs_ddone (D, NULL, jmatch, 0)) ; + /* --- Coarse decomposition --------------------------------------------- */ + wi = r ; wj = s ; /* use r and s as workspace */ + for (j = 0 ; j < n ; j++) wj [j] = -1 ; /* unmark all cols for bfs */ + for (i = 0 ; i < m ; i++) wi [i] = -1 ; /* unmark all rows for bfs */ + cs_bfs (A, n, wi, wj, q, imatch, jmatch, 1) ; /* find C1, R1 from C0*/ + ok = cs_bfs (A, m, wj, wi, p, jmatch, imatch, 3) ; /* find R3, C3 from R0*/ + if (!ok) return (cs_ddone (D, NULL, jmatch, 0)) ; + cs_unmatched (n, wj, q, cc, 0) ; /* unmatched set C0 */ + cs_matched (n, wj, imatch, p, q, cc, rr, 1, 1) ; /* set R1 and C1 */ + cs_matched (n, wj, imatch, p, q, cc, rr, 2, -1) ; /* set R2 and C2 */ + cs_matched (n, wj, imatch, p, q, cc, rr, 3, 3) ; /* set R3 and C3 */ + cs_unmatched (m, wi, p, rr, 3) ; /* unmatched set R0 */ + cs_free (jmatch) ; + /* --- Fine decomposition ----------------------------------------------- */ + pinv = cs_pinv (p, m) ; /* pinv=p' */ + if (!pinv) return (cs_ddone (D, NULL, NULL, 0)) ; + C = cs_permute (A, pinv, q, 0) ;/* C=A(p,q) (it will hold A(R2,C2)) */ + cs_free (pinv) ; + if (!C) return (cs_ddone (D, NULL, NULL, 0)) ; + Cp = C->p ; + nc = cc [3] - cc [2] ; /* delete cols C0, C1, and C3 from C */ + if (cc [2] > 0) for (j = cc [2] ; j <= cc [3] ; j++) Cp [j-cc[2]] = Cp [j] ; + C->n = nc ; + if (rr [2] - rr [1] < m) /* delete rows R0, R1, and R3 from C */ + { + cs_fkeep (C, cs_rprune, rr) ; + cnz = Cp [nc] ; + Ci = C->i ; + if (rr [1] > 0) for (k = 0 ; k < cnz ; k++) Ci [k] -= rr [1] ; + } + C->m = nc ; + scc = cs_scc (C) ; /* find strongly connected components of C*/ + if (!scc) return (cs_ddone (D, C, NULL, 0)) ; + /* --- Combine coarse and fine decompositions --------------------------- */ + ps = scc->p ; /* C(ps,ps) is the permuted matrix */ + rs = scc->r ; /* kth block is rs[k]..rs[k+1]-1 */ + nb1 = scc->nb ; /* # of blocks of A(R2,C2) */ + for (k = 0 ; k < nc ; k++) wj [k] = q [ps [k] + cc [2]] ; + for (k = 0 ; k < nc ; k++) q [k + cc [2]] = wj [k] ; + for (k = 0 ; k < nc ; k++) wi [k] = p [ps [k] + rr [1]] ; + for (k = 0 ; k < nc ; k++) p [k + rr [1]] = wi [k] ; + nb2 = 0 ; /* create the fine block partitions */ + r [0] = s [0] = 0 ; + if (cc [2] > 0) nb2++ ; /* leading coarse block A (R1, [C0 C1]) */ + for (k = 0 ; k < nb1 ; k++) /* coarse block A (R2,C2) */ + { + r [nb2] = rs [k] + rr [1] ; /* A (R2,C2) splits into nb1 fine blocks */ + s [nb2] = rs [k] + cc [2] ; + nb2++ ; + } + if (rr [2] < m) + { + r [nb2] = rr [2] ; /* trailing coarse block A ([R3 R0], C3) */ + s [nb2] = cc [3] ; + nb2++ ; + } + r [nb2] = m ; + s [nb2] = n ; + D->nb = nb2 ; + cs_dfree (scc) ; + return (cs_ddone (D, C, NULL, 1)) ; +} diff --git a/lib/csparse/cs_droptol.c b/lib/csparse/cs_droptol.c new file mode 100644 index 0000000..2afb660 --- /dev/null +++ b/lib/csparse/cs_droptol.c @@ -0,0 +1,9 @@ +#include "cs.h" +static int cs_tol (int i, int j, double aij, void *tol) +{ + return (fabs (aij) > *((double *) tol)) ; +} +int cs_droptol (cs *A, double tol) +{ + return (cs_fkeep (A, &cs_tol, &tol)) ; /* keep all large entries */ +} diff --git a/lib/csparse/cs_dropzeros.c b/lib/csparse/cs_dropzeros.c new file mode 100644 index 0000000..441e54d --- /dev/null +++ b/lib/csparse/cs_dropzeros.c @@ -0,0 +1,9 @@ +#include "cs.h" +static int cs_nonzero (int i, int j, double aij, void *other) +{ + return (aij != 0) ; +} +int cs_dropzeros (cs *A) +{ + return (cs_fkeep (A, &cs_nonzero, NULL)) ; /* keep all nonzero entries */ +} diff --git a/lib/csparse/cs_dupl.c b/lib/csparse/cs_dupl.c new file mode 100644 index 0000000..7a60d63 --- /dev/null +++ b/lib/csparse/cs_dupl.c @@ -0,0 +1,34 @@ +#include "cs.h" +/* remove duplicate entries from A */ +int cs_dupl (cs *A) +{ + int i, j, p, q, nz = 0, n, m, *Ap, *Ai, *w ; + double *Ax ; + if (!CS_CSC (A)) return (0) ; /* check inputs */ + m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ; + w = cs_malloc (m, sizeof (int)) ; /* get workspace */ + if (!w) return (0) ; /* out of memory */ + for (i = 0 ; i < m ; i++) w [i] = -1 ; /* row i not yet seen */ + for (j = 0 ; j < n ; j++) + { + q = nz ; /* column j will start at q */ + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + i = Ai [p] ; /* A(i,j) is nonzero */ + if (w [i] >= q) + { + Ax [w [i]] += Ax [p] ; /* A(i,j) is a duplicate */ + } + else + { + w [i] = nz ; /* record where row i occurs */ + Ai [nz] = i ; /* keep A(i,j) */ + Ax [nz++] = Ax [p] ; + } + } + Ap [j] = q ; /* record start of column j */ + } + Ap [n] = nz ; /* finalize A */ + cs_free (w) ; /* free workspace */ + return (cs_sprealloc (A, 0)) ; /* remove extra space from A */ +} diff --git a/lib/csparse/cs_entry.c b/lib/csparse/cs_entry.c new file mode 100644 index 0000000..98f6788 --- /dev/null +++ b/lib/csparse/cs_entry.c @@ -0,0 +1,13 @@ +#include "cs.h" +/* add an entry to a triplet matrix; return 1 if ok, 0 otherwise */ +int cs_entry (cs *T, int i, int j, double x) +{ + if (!CS_TRIPLET (T) || i < 0 || j < 0) return (0) ; /* check inputs */ + if (T->nz >= T->nzmax && !cs_sprealloc (T,2*(T->nzmax))) return (0) ; + if (T->x) T->x [T->nz] = x ; + T->i [T->nz] = i ; + T->p [T->nz++] = j ; + T->m = CS_MAX (T->m, i+1) ; + T->n = CS_MAX (T->n, j+1) ; + return (1) ; +} diff --git a/lib/csparse/cs_ereach.c b/lib/csparse/cs_ereach.c new file mode 100644 index 0000000..c54213b --- /dev/null +++ b/lib/csparse/cs_ereach.c @@ -0,0 +1,23 @@ +#include "cs.h" +/* find nonzero pattern of Cholesky L(k,1:k-1) using etree and triu(A(:,k)) */ +int cs_ereach (const cs *A, int k, const int *parent, int *s, int *w) +{ + int i, p, n, len, top, *Ap, *Ai ; + if (!CS_CSC (A) || !parent || !s || !w) return (-1) ; /* check inputs */ + top = n = A->n ; Ap = A->p ; Ai = A->i ; + CS_MARK (w, k) ; /* mark node k as visited */ + for (p = Ap [k] ; p < Ap [k+1] ; p++) + { + i = Ai [p] ; /* A(i,k) is nonzero */ + if (i > k) continue ; /* only use upper triangular part of A */ + for (len = 0 ; !CS_MARKED (w,i) ; i = parent [i]) /* traverse up etree*/ + { + s [len++] = i ; /* L(k,i) is nonzero */ + CS_MARK (w, i) ; /* mark i as visited */ + } + while (len > 0) s [--top] = s [--len] ; /* push path onto stack */ + } + for (p = top ; p < n ; p++) CS_MARK (w, s [p]) ; /* unmark all nodes */ + CS_MARK (w, k) ; /* unmark node k */ + return (top) ; /* s [top..n-1] contains pattern of L(k,:)*/ +} diff --git a/lib/csparse/cs_etree.c b/lib/csparse/cs_etree.c new file mode 100644 index 0000000..30712bc --- /dev/null +++ b/lib/csparse/cs_etree.c @@ -0,0 +1,30 @@ +#include "cs.h" +/* compute the etree of A (using triu(A), or A'A without forming A'A */ +int *cs_etree (const cs *A, int ata) +{ + int i, k, p, m, n, inext, *Ap, *Ai, *w, *parent, *ancestor, *prev ; + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; + parent = cs_malloc (n, sizeof (int)) ; /* allocate result */ + w = cs_malloc (n + (ata ? m : 0), sizeof (int)) ; /* get workspace */ + if (!w || !parent) return (cs_idone (parent, NULL, w, 0)) ; + ancestor = w ; prev = w + n ; + if (ata) for (i = 0 ; i < m ; i++) prev [i] = -1 ; + for (k = 0 ; k < n ; k++) + { + parent [k] = -1 ; /* node k has no parent yet */ + ancestor [k] = -1 ; /* nor does k have an ancestor */ + for (p = Ap [k] ; p < Ap [k+1] ; p++) + { + i = ata ? (prev [Ai [p]]) : (Ai [p]) ; + for ( ; i != -1 && i < k ; i = inext) /* traverse from i to k */ + { + inext = ancestor [i] ; /* inext = ancestor of i */ + ancestor [i] = k ; /* path compression */ + if (inext == -1) parent [i] = k ; /* no anc., parent is k */ + } + if (ata) prev [Ai [p]] = k ; + } + } + return (cs_idone (parent, NULL, w, 1)) ; +} diff --git a/lib/csparse/cs_fkeep.c b/lib/csparse/cs_fkeep.c new file mode 100644 index 0000000..0e32051 --- /dev/null +++ b/lib/csparse/cs_fkeep.c @@ -0,0 +1,25 @@ +#include "cs.h" +/* drop entries for which fkeep(A(i,j)) is false; return nz if OK, else -1 */ +int cs_fkeep (cs *A, int (*fkeep) (int, int, double, void *), void *other) +{ + int j, p, nz = 0, n, *Ap, *Ai ; + double *Ax ; + if (!CS_CSC (A) || !fkeep) return (-1) ; /* check inputs */ + n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ; + for (j = 0 ; j < n ; j++) + { + p = Ap [j] ; /* get current location of col j */ + Ap [j] = nz ; /* record new location of col j */ + for ( ; p < Ap [j+1] ; p++) + { + if (fkeep (Ai [p], j, Ax ? Ax [p] : 1, other)) + { + if (Ax) Ax [nz] = Ax [p] ; /* keep A(i,j) */ + Ai [nz++] = Ai [p] ; + } + } + } + Ap [n] = nz ; /* finalize A */ + cs_sprealloc (A, 0) ; /* remove extra space from A */ + return (nz) ; +} diff --git a/lib/csparse/cs_gaxpy.c b/lib/csparse/cs_gaxpy.c new file mode 100644 index 0000000..580f600 --- /dev/null +++ b/lib/csparse/cs_gaxpy.c @@ -0,0 +1,17 @@ +#include "cs.h" +/* y = A*x+y */ +int cs_gaxpy (const cs *A, const double *x, double *y) +{ + int p, j, n, *Ap, *Ai ; + double *Ax ; + if (!CS_CSC (A) || !x || !y) return (0) ; /* check inputs */ + n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ; + for (j = 0 ; j < n ; j++) + { + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + y [Ai [p]] += Ax [p] * x [j] ; + } + } + return (1) ; +} diff --git a/lib/csparse/cs_happly.c b/lib/csparse/cs_happly.c new file mode 100644 index 0000000..b2d4346 --- /dev/null +++ b/lib/csparse/cs_happly.c @@ -0,0 +1,19 @@ +#include "cs.h" +/* apply the ith Householder vector to x */ +int cs_happly (const cs *V, int i, double beta, double *x) +{ + int p, *Vp, *Vi ; + double *Vx, tau = 0 ; + if (!CS_CSC (V) || !x) return (0) ; /* check inputs */ + Vp = V->p ; Vi = V->i ; Vx = V->x ; + for (p = Vp [i] ; p < Vp [i+1] ; p++) /* tau = v'*x */ + { + tau += Vx [p] * x [Vi [p]] ; + } + tau *= beta ; /* tau = beta*(v'*x) */ + for (p = Vp [i] ; p < Vp [i+1] ; p++) /* x = x - v*tau */ + { + x [Vi [p]] -= Vx [p] * tau ; + } + return (1) ; +} diff --git a/lib/csparse/cs_house.c b/lib/csparse/cs_house.c new file mode 100644 index 0000000..5935f4c --- /dev/null +++ b/lib/csparse/cs_house.c @@ -0,0 +1,23 @@ +#include "cs.h" +/* create a Householder reflection [v,beta,s]=house(x), overwrite x with v, + * where (I-beta*v*v')*x = s*e1. See Algo 5.1.1, Golub & Van Loan, 3rd ed. */ +double cs_house (double *x, double *beta, int n) +{ + double s, sigma = 0 ; + int i ; + if (!x || !beta) return (-1) ; /* check inputs */ + for (i = 1 ; i < n ; i++) sigma += x [i] * x [i] ; + if (sigma == 0) + { + s = fabs (x [0]) ; /* s = |x(0)| */ + (*beta) = (x [0] <= 0) ? 2 : 0 ; + x [0] = 1 ; + } + else + { + s = sqrt (x [0] * x [0] + sigma) ; /* s = norm (x) */ + x [0] = (x [0] <= 0) ? (x [0] - s) : (-sigma / (x [0] + s)) ; + (*beta) = -1. / (s * x [0]) ; + } + return (s) ; +} diff --git a/lib/csparse/cs_ipvec.c b/lib/csparse/cs_ipvec.c new file mode 100644 index 0000000..cb1e9a9 --- /dev/null +++ b/lib/csparse/cs_ipvec.c @@ -0,0 +1,9 @@ +#include "cs.h" +/* x(p) = b, for dense vectors x and b; p=NULL denotes identity */ +int cs_ipvec (const int *p, const double *b, double *x, int n) +{ + int k ; + if (!x || !b) return (0) ; /* check inputs */ + for (k = 0 ; k < n ; k++) x [p ? p [k] : k] = b [k] ; + return (1) ; +} diff --git a/lib/csparse/cs_leaf.c b/lib/csparse/cs_leaf.c new file mode 100644 index 0000000..1601cb9 --- /dev/null +++ b/lib/csparse/cs_leaf.c @@ -0,0 +1,22 @@ +#include "cs.h" +/* consider A(i,j), node j in ith row subtree and return lca(jprev,j) */ +int cs_leaf (int i, int j, const int *first, int *maxfirst, int *prevleaf, + int *ancestor, int *jleaf) +{ + int q, s, sparent, jprev ; + if (!first || !maxfirst || !prevleaf || !ancestor || !jleaf) return (-1) ; + *jleaf = 0 ; + if (i <= j || first [j] <= maxfirst [i]) return (-1) ; /* j not a leaf */ + maxfirst [i] = first [j] ; /* update max first[j] seen so far */ + jprev = prevleaf [i] ; /* jprev = previous leaf of ith subtree */ + prevleaf [i] = j ; + *jleaf = (jprev == -1) ? 1: 2 ; /* j is first or subsequent leaf */ + if (*jleaf == 1) return (i) ; /* if 1st leaf, q = root of ith subtree */ + for (q = jprev ; q != ancestor [q] ; q = ancestor [q]) ; + for (s = jprev ; s != q ; s = sparent) + { + sparent = ancestor [s] ; /* path compression */ + ancestor [s] = q ; + } + return (q) ; /* q = least common ancester (jprev,j) */ +} diff --git a/lib/csparse/cs_load.c b/lib/csparse/cs_load.c new file mode 100644 index 0000000..e1e8719 --- /dev/null +++ b/lib/csparse/cs_load.c @@ -0,0 +1,15 @@ +#include "cs.h" +/* load a triplet matrix from a file */ +cs *cs_load (FILE *f) +{ + int i, j ; + double x ; + cs *T ; + if (!f) return (NULL) ; /* check inputs */ + T = cs_spalloc (0, 0, 1, 1, 1) ; /* allocate result */ + while (fscanf (f, "%d %d %lg\n", &i, &j, &x) == 3) + { + if (!cs_entry (T, i, j, x)) return (cs_spfree (T)) ; + } + return (T) ; +} diff --git a/lib/csparse/cs_lsolve.c b/lib/csparse/cs_lsolve.c new file mode 100644 index 0000000..69376ac --- /dev/null +++ b/lib/csparse/cs_lsolve.c @@ -0,0 +1,18 @@ +#include "cs.h" +/* solve Lx=b where x and b are dense. x=b on input, solution on output. */ +int cs_lsolve (const cs *L, double *x) +{ + int p, j, n, *Lp, *Li ; + double *Lx ; + if (!CS_CSC (L) || !x) return (0) ; /* check inputs */ + n = L->n ; Lp = L->p ; Li = L->i ; Lx = L->x ; + for (j = 0 ; j < n ; j++) + { + x [j] /= Lx [Lp [j]] ; + for (p = Lp [j]+1 ; p < Lp [j+1] ; p++) + { + x [Li [p]] -= Lx [p] * x [j] ; + } + } + return (1) ; +} diff --git a/lib/csparse/cs_ltsolve.c b/lib/csparse/cs_ltsolve.c new file mode 100644 index 0000000..fc2ee4f --- /dev/null +++ b/lib/csparse/cs_ltsolve.c @@ -0,0 +1,18 @@ +#include "cs.h" +/* solve L'x=b where x and b are dense. x=b on input, solution on output. */ +int cs_ltsolve (const cs *L, double *x) +{ + int p, j, n, *Lp, *Li ; + double *Lx ; + if (!CS_CSC (L) || !x) return (0) ; /* check inputs */ + n = L->n ; Lp = L->p ; Li = L->i ; Lx = L->x ; + for (j = n-1 ; j >= 0 ; j--) + { + for (p = Lp [j]+1 ; p < Lp [j+1] ; p++) + { + x [j] -= Lx [p] * x [Li [p]] ; + } + x [j] /= Lx [Lp [j]] ; + } + return (1) ; +} diff --git a/lib/csparse/cs_lu.c b/lib/csparse/cs_lu.c new file mode 100644 index 0000000..c8883c6 --- /dev/null +++ b/lib/csparse/cs_lu.c @@ -0,0 +1,86 @@ +#include "cs.h" +/* [L,U,pinv]=lu(A, [q lnz unz]). lnz and unz can be guess */ +csn *cs_lu (const cs *A, const css *S, double tol) +{ + cs *L, *U ; + csn *N ; + double pivot, *Lx, *Ux, *x, a, t ; + int *Lp, *Li, *Up, *Ui, *pinv, *xi, *q, n, ipiv, k, top, p, i, col, lnz,unz; + if (!CS_CSC (A) || !S) return (NULL) ; /* check inputs */ + n = A->n ; + q = S->q ; lnz = S->lnz ; unz = S->unz ; + x = cs_malloc (n, sizeof (double)) ; /* get double workspace */ + xi = cs_malloc (2*n, sizeof (int)) ; /* get int workspace */ + N = cs_calloc (1, sizeof (csn)) ; /* allocate result */ + if (!x || !xi || !N) return (cs_ndone (N, NULL, xi, x, 0)) ; + N->L = L = cs_spalloc (n, n, lnz, 1, 0) ; /* allocate result L */ + N->U = U = cs_spalloc (n, n, unz, 1, 0) ; /* allocate result U */ + N->pinv = pinv = cs_malloc (n, sizeof (int)) ; /* allocate result pinv */ + if (!L || !U || !pinv) return (cs_ndone (N, NULL, xi, x, 0)) ; + Lp = L->p ; Up = U->p ; + for (i = 0 ; i < n ; i++) x [i] = 0 ; /* clear workspace */ + for (i = 0 ; i < n ; i++) pinv [i] = -1 ; /* no rows pivotal yet */ + for (k = 0 ; k <= n ; k++) Lp [k] = 0 ; /* no cols of L yet */ + lnz = unz = 0 ; + for (k = 0 ; k < n ; k++) /* compute L(:,k) and U(:,k) */ + { + /* --- Triangular solve --------------------------------------------- */ + Lp [k] = lnz ; /* L(:,k) starts here */ + Up [k] = unz ; /* U(:,k) starts here */ + if ((lnz + n > L->nzmax && !cs_sprealloc (L, 2*L->nzmax + n)) || + (unz + n > U->nzmax && !cs_sprealloc (U, 2*U->nzmax + n))) + { + return (cs_ndone (N, NULL, xi, x, 0)) ; + } + Li = L->i ; Lx = L->x ; Ui = U->i ; Ux = U->x ; + col = q ? (q [k]) : k ; + top = cs_spsolve (L, A, col, xi, x, pinv, 1) ; /* x = L\A(:,col) */ + /* --- Find pivot --------------------------------------------------- */ + ipiv = -1 ; + a = -1 ; + for (p = top ; p < n ; p++) + { + i = xi [p] ; /* x(i) is nonzero */ + if (pinv [i] < 0) /* row i is not yet pivotal */ + { + if ((t = fabs (x [i])) > a) + { + a = t ; /* largest pivot candidate so far */ + ipiv = i ; + } + } + else /* x(i) is the entry U(pinv[i],k) */ + { + Ui [unz] = pinv [i] ; + Ux [unz++] = x [i] ; + } + } + if (ipiv == -1 || a <= 0) return (cs_ndone (N, NULL, xi, x, 0)) ; + if (pinv [col] < 0 && fabs (x [col]) >= a*tol) ipiv = col ; + /* --- Divide by pivot ---------------------------------------------- */ + pivot = x [ipiv] ; /* the chosen pivot */ + Ui [unz] = k ; /* last entry in U(:,k) is U(k,k) */ + Ux [unz++] = pivot ; + pinv [ipiv] = k ; /* ipiv is the kth pivot row */ + Li [lnz] = ipiv ; /* first entry in L(:,k) is L(k,k) = 1 */ + Lx [lnz++] = 1 ; + for (p = top ; p < n ; p++) /* L(k+1:n,k) = x / pivot */ + { + i = xi [p] ; + if (pinv [i] < 0) /* x(i) is an entry in L(:,k) */ + { + Li [lnz] = i ; /* save unpermuted row in L */ + Lx [lnz++] = x [i] / pivot ; /* scale pivot column */ + } + x [i] = 0 ; /* x [0..n-1] = 0 for next k */ + } + } + /* --- Finalize L and U ------------------------------------------------- */ + Lp [n] = lnz ; + Up [n] = unz ; + Li = L->i ; /* fix row indices of L for final pinv */ + for (p = 0 ; p < lnz ; p++) Li [p] = pinv [Li [p]] ; + cs_sprealloc (L, 0) ; /* remove extra space from L and U */ + cs_sprealloc (U, 0) ; + return (cs_ndone (N, NULL, xi, x, 1)) ; /* success */ +} diff --git a/lib/csparse/cs_lusol.c b/lib/csparse/cs_lusol.c new file mode 100644 index 0000000..e7d79a1 --- /dev/null +++ b/lib/csparse/cs_lusol.c @@ -0,0 +1,26 @@ +#include "cs.h" +/* x=A\b where A is unsymmetric; b overwritten with solution */ +int cs_lusol (int order, const cs *A, double *b, double tol) +{ + double *x ; + css *S ; + csn *N ; + int n, ok ; + if (!CS_CSC (A) || !b) return (0) ; /* check inputs */ + n = A->n ; + S = cs_sqr (order, A, 0) ; /* ordering and symbolic analysis */ + N = cs_lu (A, S, tol) ; /* numeric LU factorization */ + x = cs_malloc (n, sizeof (double)) ; /* get workspace */ + ok = (S && N && x) ; + if (ok) + { + cs_ipvec (N->pinv, b, x, n) ; /* x = b(p) */ + cs_lsolve (N->L, x) ; /* x = L\x */ + cs_usolve (N->U, x) ; /* x = U\x */ + cs_ipvec (S->q, x, b, n) ; /* b(q) = x */ + } + cs_free (x) ; + cs_sfree (S) ; + cs_nfree (N) ; + return (ok) ; +} diff --git a/lib/csparse/cs_malloc.c b/lib/csparse/cs_malloc.c new file mode 100644 index 0000000..4d8e2e1 --- /dev/null +++ b/lib/csparse/cs_malloc.c @@ -0,0 +1,35 @@ +#include "cs.h" +#ifdef MATLAB_MEX_FILE +#define malloc mxMalloc +#define free mxFree +#define realloc mxRealloc +#define calloc mxCalloc +#endif + +/* wrapper for malloc */ +void *cs_malloc (int n, size_t size) +{ + return (malloc (CS_MAX (n,1) * size)) ; +} + +/* wrapper for calloc */ +void *cs_calloc (int n, size_t size) +{ + return (calloc (CS_MAX (n,1), size)) ; +} + +/* wrapper for free */ +void *cs_free (void *p) +{ + if (p) free (p) ; /* free p if it is not already NULL */ + return (NULL) ; /* return NULL to simplify the use of cs_free */ +} + +/* wrapper for realloc */ +void *cs_realloc (void *p, int n, size_t size, int *ok) +{ + void *pnew ; + pnew = realloc (p, CS_MAX (n,1) * size) ; /* realloc the block */ + *ok = (pnew != NULL) ; /* realloc fails if pnew is NULL */ + return ((*ok) ? pnew : p) ; /* return original p if failure */ +} diff --git a/lib/csparse/cs_maxtrans.c b/lib/csparse/cs_maxtrans.c new file mode 100644 index 0000000..7c8ea65 --- /dev/null +++ b/lib/csparse/cs_maxtrans.c @@ -0,0 +1,92 @@ +#include "cs.h" +/* find an augmenting path starting at column k and extend the match if found */ +static void cs_augment (int k, const cs *A, int *jmatch, int *cheap, int *w, + int *js, int *is, int *ps) +{ + int found = 0, p, i = -1, *Ap = A->p, *Ai = A->i, head = 0, j ; + js [0] = k ; /* start with just node k in jstack */ + while (head >= 0) + { + /* --- Start (or continue) depth-first-search at node j ------------- */ + j = js [head] ; /* get j from top of jstack */ + if (w [j] != k) /* 1st time j visited for kth path */ + { + w [j] = k ; /* mark j as visited for kth path */ + for (p = cheap [j] ; p < Ap [j+1] && !found ; p++) + { + i = Ai [p] ; /* try a cheap assignment (i,j) */ + found = (jmatch [i] == -1) ; + } + cheap [j] = p ; /* start here next time j is traversed*/ + if (found) + { + is [head] = i ; /* column j matched with row i */ + break ; /* end of augmenting path */ + } + ps [head] = Ap [j] ; /* no cheap match: start dfs for j */ + } + /* --- Depth-first-search of neighbors of j ------------------------- */ + for (p = ps [head] ; p < Ap [j+1] ; p++) + { + i = Ai [p] ; /* consider row i */ + if (w [jmatch [i]] == k) continue ; /* skip jmatch [i] if marked */ + ps [head] = p + 1 ; /* pause dfs of node j */ + is [head] = i ; /* i will be matched with j if found */ + js [++head] = jmatch [i] ; /* start dfs at column jmatch [i] */ + break ; + } + if (p == Ap [j+1]) head-- ; /* node j is done; pop from stack */ + } /* augment the match if path found: */ + if (found) for (p = head ; p >= 0 ; p--) jmatch [is [p]] = js [p] ; +} + +/* find a maximum transveral */ +int *cs_maxtrans (const cs *A, int seed) /*[jmatch [0..m-1]; imatch [0..n-1]]*/ +{ + int i, j, k, n, m, p, n2 = 0, m2 = 0, *Ap, *jimatch, *w, *cheap, *js, *is, + *ps, *Ai, *Cp, *jmatch, *imatch, *q ; + cs *C ; + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + n = A->n ; m = A->m ; Ap = A->p ; Ai = A->i ; + w = jimatch = cs_calloc (m+n, sizeof (int)) ; /* allocate result */ + if (!jimatch) return (NULL) ; + for (k = 0, j = 0 ; j < n ; j++) /* count nonempty rows and columns */ + { + n2 += (Ap [j] < Ap [j+1]) ; + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + w [Ai [p]] = 1 ; + k += (j == Ai [p]) ; /* count entries already on diagonal */ + } + } + if (k == CS_MIN (m,n)) /* quick return if diagonal zero-free */ + { + jmatch = jimatch ; imatch = jimatch + m ; + for (i = 0 ; i < k ; i++) jmatch [i] = i ; + for ( ; i < m ; i++) jmatch [i] = -1 ; + for (j = 0 ; j < k ; j++) imatch [j] = j ; + for ( ; j < n ; j++) imatch [j] = -1 ; + return (cs_idone (jimatch, NULL, NULL, 1)) ; + } + for (i = 0 ; i < m ; i++) m2 += w [i] ; + C = (m2 < n2) ? cs_transpose (A,0) : ((cs *) A) ; /* transpose if needed */ + if (!C) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, NULL, 0)) ; + n = C->n ; m = C->m ; Cp = C->p ; + jmatch = (m2 < n2) ? jimatch + n : jimatch ; + imatch = (m2 < n2) ? jimatch : jimatch + m ; + w = cs_malloc (5*n, sizeof (int)) ; /* get workspace */ + if (!w) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 0)) ; + cheap = w + n ; js = w + 2*n ; is = w + 3*n ; ps = w + 4*n ; + for (j = 0 ; j < n ; j++) cheap [j] = Cp [j] ; /* for cheap assignment */ + for (j = 0 ; j < n ; j++) w [j] = -1 ; /* all columns unflagged */ + for (i = 0 ; i < m ; i++) jmatch [i] = -1 ; /* nothing matched yet */ + q = cs_randperm (n, seed) ; /* q = random permutation */ + for (k = 0 ; k < n ; k++) /* augment, starting at column q[k] */ + { + cs_augment (q ? q [k]: k, C, jmatch, cheap, w, js, is, ps) ; + } + cs_free (q) ; + for (j = 0 ; j < n ; j++) imatch [j] = -1 ; /* find row match */ + for (i = 0 ; i < m ; i++) if (jmatch [i] >= 0) imatch [jmatch [i]] = i ; + return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 1)) ; +} diff --git a/lib/csparse/cs_multiply.c b/lib/csparse/cs_multiply.c new file mode 100644 index 0000000..c52ebea --- /dev/null +++ b/lib/csparse/cs_multiply.c @@ -0,0 +1,35 @@ +#include "cs.h" +/* C = A*B */ +cs *cs_multiply (const cs *A, const cs *B) +{ + int p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values, *Bi ; + double *x, *Bx, *Cx ; + cs *C ; + if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ; /* check inputs */ + if (A->n != B->m) return (NULL) ; + m = A->m ; anz = A->p [A->n] ; + n = B->n ; Bp = B->p ; Bi = B->i ; Bx = B->x ; bnz = Bp [n] ; + w = cs_calloc (m, sizeof (int)) ; /* get workspace */ + values = (A->x != NULL) && (Bx != NULL) ; + x = values ? cs_malloc (m, sizeof (double)) : NULL ; /* get workspace */ + C = cs_spalloc (m, n, anz + bnz, values, 0) ; /* allocate result */ + if (!C || !w || (values && !x)) return (cs_done (C, w, x, 0)) ; + Cp = C->p ; + for (j = 0 ; j < n ; j++) + { + if (nz + m > C->nzmax && !cs_sprealloc (C, 2*(C->nzmax)+m)) + { + return (cs_done (C, w, x, 0)) ; /* out of memory */ + } + Ci = C->i ; Cx = C->x ; /* C->i and C->x may be reallocated */ + Cp [j] = nz ; /* column j of C starts here */ + for (p = Bp [j] ; p < Bp [j+1] ; p++) + { + nz = cs_scatter (A, Bi [p], Bx ? Bx [p] : 1, w, x, j+1, C, nz) ; + } + if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Ci [p]] ; + } + Cp [n] = nz ; /* finalize the last column of C */ + cs_sprealloc (C, 0) ; /* remove extra space from C */ + return (cs_done (C, w, x, 1)) ; /* success; free workspace, return C */ +} diff --git a/lib/csparse/cs_norm.c b/lib/csparse/cs_norm.c new file mode 100644 index 0000000..24f8a2e --- /dev/null +++ b/lib/csparse/cs_norm.c @@ -0,0 +1,15 @@ +#include "cs.h" +/* 1-norm of a sparse matrix = max (sum (abs (A))), largest column sum */ +double cs_norm (const cs *A) +{ + int p, j, n, *Ap ; + double *Ax, norm = 0, s ; + if (!CS_CSC (A) || !A->x) return (-1) ; /* check inputs */ + n = A->n ; Ap = A->p ; Ax = A->x ; + for (j = 0 ; j < n ; j++) + { + for (s = 0, p = Ap [j] ; p < Ap [j+1] ; p++) s += fabs (Ax [p]) ; + norm = CS_MAX (norm, s) ; + } + return (norm) ; +} diff --git a/lib/csparse/cs_permute.c b/lib/csparse/cs_permute.c new file mode 100644 index 0000000..f08e4d1 --- /dev/null +++ b/lib/csparse/cs_permute.c @@ -0,0 +1,25 @@ +#include "cs.h" +/* C = A(p,q) where p and q are permutations of 0..m-1 and 0..n-1. */ +cs *cs_permute (const cs *A, const int *pinv, const int *q, int values) +{ + int t, j, k, nz = 0, m, n, *Ap, *Ai, *Cp, *Ci ; + double *Cx, *Ax ; + cs *C ; + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ; + C = cs_spalloc (m, n, Ap [n], values && Ax != NULL, 0) ; /* alloc result */ + if (!C) return (cs_done (C, NULL, NULL, 0)) ; /* out of memory */ + Cp = C->p ; Ci = C->i ; Cx = C->x ; + for (k = 0 ; k < n ; k++) + { + Cp [k] = nz ; /* column k of C is column q[k] of A */ + j = q ? (q [k]) : k ; + for (t = Ap [j] ; t < Ap [j+1] ; t++) + { + if (Cx) Cx [nz] = Ax [t] ; /* row i of A is row pinv[i] of C */ + Ci [nz++] = pinv ? (pinv [Ai [t]]) : Ai [t] ; + } + } + Cp [n] = nz ; /* finalize the last column of C */ + return (cs_done (C, NULL, NULL, 1)) ; +} diff --git a/lib/csparse/cs_pinv.c b/lib/csparse/cs_pinv.c new file mode 100644 index 0000000..0583e77 --- /dev/null +++ b/lib/csparse/cs_pinv.c @@ -0,0 +1,11 @@ +#include "cs.h" +/* pinv = p', or p = pinv' */ +int *cs_pinv (int const *p, int n) +{ + int k, *pinv ; + if (!p) return (NULL) ; /* p = NULL denotes identity */ + pinv = cs_malloc (n, sizeof (int)) ; /* allocate result */ + if (!pinv) return (NULL) ; /* out of memory */ + for (k = 0 ; k < n ; k++) pinv [p [k]] = k ;/* invert the permutation */ + return (pinv) ; /* return result */ +} diff --git a/lib/csparse/cs_post.c b/lib/csparse/cs_post.c new file mode 100644 index 0000000..80e66e6 --- /dev/null +++ b/lib/csparse/cs_post.c @@ -0,0 +1,24 @@ +#include "cs.h" +/* post order a forest */ +int *cs_post (const int *parent, int n) +{ + int j, k = 0, *post, *w, *head, *next, *stack ; + if (!parent) return (NULL) ; /* check inputs */ + post = cs_malloc (n, sizeof (int)) ; /* allocate result */ + w = cs_malloc (3*n, sizeof (int)) ; /* get workspace */ + if (!w || !post) return (cs_idone (post, NULL, w, 0)) ; + head = w ; next = w + n ; stack = w + 2*n ; + for (j = 0 ; j < n ; j++) head [j] = -1 ; /* empty linked lists */ + for (j = n-1 ; j >= 0 ; j--) /* traverse nodes in reverse order*/ + { + if (parent [j] == -1) continue ; /* j is a root */ + next [j] = head [parent [j]] ; /* add j to list of its parent */ + head [parent [j]] = j ; + } + for (j = 0 ; j < n ; j++) + { + if (parent [j] != -1) continue ; /* skip j if it is not a root */ + k = cs_tdfs (j, k, head, next, post, stack) ; + } + return (cs_idone (post, NULL, w, 1)) ; /* success; free w, return post */ +} diff --git a/lib/csparse/cs_print.c b/lib/csparse/cs_print.c new file mode 100644 index 0000000..8e3f33e --- /dev/null +++ b/lib/csparse/cs_print.c @@ -0,0 +1,36 @@ +#include "cs.h" +/* print a sparse matrix */ +int cs_print (const cs *A, int brief) +{ + int p, j, m, n, nzmax, nz, *Ap, *Ai ; + double *Ax ; + if (!A) { printf ("(null)\n") ; return (0) ; } + m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ; + nzmax = A->nzmax ; nz = A->nz ; + printf ("CSparse Version %d.%d.%d, %s. %s\n", CS_VER, CS_SUBVER, + CS_SUBSUB, CS_DATE, CS_COPYRIGHT) ; + if (nz < 0) + { + printf ("%d-by-%d, nzmax: %d nnz: %d, 1-norm: %g\n", m, n, nzmax, + Ap [n], cs_norm (A)) ; + for (j = 0 ; j < n ; j++) + { + printf (" col %d : locations %d to %d\n", j, Ap [j], Ap [j+1]-1); + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + printf (" %d : %g\n", Ai [p], Ax ? Ax [p] : 1) ; + if (brief && p > 20) { printf (" ...\n") ; return (1) ; } + } + } + } + else + { + printf ("triplet: %d-by-%d, nzmax: %d nnz: %d\n", m, n, nzmax, nz) ; + for (p = 0 ; p < nz ; p++) + { + printf (" %d %d : %g\n", Ai [p], Ap [p], Ax ? Ax [p] : 1) ; + if (brief && p > 20) { printf (" ...\n") ; return (1) ; } + } + } + return (1) ; +} diff --git a/lib/csparse/cs_pvec.c b/lib/csparse/cs_pvec.c new file mode 100644 index 0000000..03c3cc6 --- /dev/null +++ b/lib/csparse/cs_pvec.c @@ -0,0 +1,9 @@ +#include "cs.h" +/* x = b(p), for dense vectors x and b; p=NULL denotes identity */ +int cs_pvec (const int *p, const double *b, double *x, int n) +{ + int k ; + if (!x || !b) return (0) ; /* check inputs */ + for (k = 0 ; k < n ; k++) x [k] = b [p ? p [k] : k] ; + return (1) ; +} diff --git a/lib/csparse/cs_qr.c b/lib/csparse/cs_qr.c new file mode 100644 index 0000000..eb1774d --- /dev/null +++ b/lib/csparse/cs_qr.c @@ -0,0 +1,73 @@ +#include "cs.h" +/* sparse QR factorization [V,beta,pinv,R] = qr (A) */ +csn *cs_qr (const cs *A, const css *S) +{ + double *Rx, *Vx, *Ax, *x, *Beta ; + int i, k, p, m, n, vnz, p1, top, m2, len, col, rnz, *s, *leftmost, *Ap, *Ai, + *parent, *Rp, *Ri, *Vp, *Vi, *w, *pinv, *q ; + cs *R, *V ; + csn *N ; + if (!CS_CSC (A) || !S) return (NULL) ; + m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ; + q = S->q ; parent = S->parent ; pinv = S->pinv ; m2 = S->m2 ; + vnz = S->lnz ; rnz = S->unz ; leftmost = S->leftmost ; + w = cs_malloc (m2+n, sizeof (int)) ; /* get int workspace */ + x = cs_malloc (m2, sizeof (double)) ; /* get double workspace */ + N = cs_calloc (1, sizeof (csn)) ; /* allocate result */ + if (!w || !x || !N) return (cs_ndone (N, NULL, w, x, 0)) ; + s = w + m2 ; /* s is size n */ + for (k = 0 ; k < m2 ; k++) x [k] = 0 ; /* clear workspace x */ + N->L = V = cs_spalloc (m2, n, vnz, 1, 0) ; /* allocate result V */ + N->U = R = cs_spalloc (m2, n, rnz, 1, 0) ; /* allocate result R */ + N->B = Beta = cs_malloc (n, sizeof (double)) ; /* allocate result Beta */ + if (!R || !V || !Beta) return (cs_ndone (N, NULL, w, x, 0)) ; + Rp = R->p ; Ri = R->i ; Rx = R->x ; + Vp = V->p ; Vi = V->i ; Vx = V->x ; + for (i = 0 ; i < m2 ; i++) w [i] = -1 ; /* clear w, to mark nodes */ + rnz = 0 ; vnz = 0 ; + for (k = 0 ; k < n ; k++) /* compute V and R */ + { + Rp [k] = rnz ; /* R(:,k) starts here */ + Vp [k] = p1 = vnz ; /* V(:,k) starts here */ + w [k] = k ; /* add V(k,k) to pattern of V */ + Vi [vnz++] = k ; + top = n ; + col = q ? q [k] : k ; + for (p = Ap [col] ; p < Ap [col+1] ; p++) /* find R(:,k) pattern */ + { + i = leftmost [Ai [p]] ; /* i = min(find(A(i,q))) */ + for (len = 0 ; w [i] != k ; i = parent [i]) /* traverse up to k */ + { + s [len++] = i ; + w [i] = k ; + } + while (len > 0) s [--top] = s [--len] ; /* push path on stack */ + i = pinv [Ai [p]] ; /* i = permuted row of A(:,col) */ + x [i] = Ax [p] ; /* x (i) = A(:,col) */ + if (i > k && w [i] < k) /* pattern of V(:,k) = x (k+1:m) */ + { + Vi [vnz++] = i ; /* add i to pattern of V(:,k) */ + w [i] = k ; + } + } + for (p = top ; p < n ; p++) /* for each i in pattern of R(:,k) */ + { + i = s [p] ; /* R(i,k) is nonzero */ + cs_happly (V, i, Beta [i], x) ; /* apply (V(i),Beta(i)) to x */ + Ri [rnz] = i ; /* R(i,k) = x(i) */ + Rx [rnz++] = x [i] ; + x [i] = 0 ; + if (parent [i] == k) vnz = cs_scatter (V, i, 0, w, NULL, k, V, vnz); + } + for (p = p1 ; p < vnz ; p++) /* gather V(:,k) = x */ + { + Vx [p] = x [Vi [p]] ; + x [Vi [p]] = 0 ; + } + Ri [rnz] = k ; /* R(k,k) = norm (x) */ + Rx [rnz++] = cs_house (Vx+p1, Beta+k, vnz-p1) ; /* [v,beta]=house(x) */ + } + Rp [n] = rnz ; /* finalize R */ + Vp [n] = vnz ; /* finalize V */ + return (cs_ndone (N, NULL, w, x, 1)) ; /* success */ +} diff --git a/lib/csparse/cs_qrsol.c b/lib/csparse/cs_qrsol.c new file mode 100644 index 0000000..7d07966 --- /dev/null +++ b/lib/csparse/cs_qrsol.c @@ -0,0 +1,53 @@ +#include "cs.h" +/* x=A\b where A can be rectangular; b overwritten with solution */ +int cs_qrsol (int order, const cs *A, double *b) +{ + double *x ; + css *S ; + csn *N ; + cs *AT = NULL ; + int k, m, n, ok ; + if (!CS_CSC (A) || !b) return (0) ; /* check inputs */ + n = A->n ; + m = A->m ; + if (m >= n) + { + S = cs_sqr (order, A, 1) ; /* ordering and symbolic analysis */ + N = cs_qr (A, S) ; /* numeric QR factorization */ + x = cs_calloc (S ? S->m2 : 1, sizeof (double)) ; /* get workspace */ + ok = (S && N && x) ; + if (ok) + { + cs_ipvec (S->pinv, b, x, m) ; /* x(0:m-1) = b(p(0:m-1) */ + for (k = 0 ; k < n ; k++) /* apply Householder refl. to x */ + { + cs_happly (N->L, k, N->B [k], x) ; + } + cs_usolve (N->U, x) ; /* x = R\x */ + cs_ipvec (S->q, x, b, n) ; /* b(q(0:n-1)) = x(0:n-1) */ + } + } + else + { + AT = cs_transpose (A, 1) ; /* Ax=b is underdetermined */ + S = cs_sqr (order, AT, 1) ; /* ordering and symbolic analysis */ + N = cs_qr (AT, S) ; /* numeric QR factorization of A' */ + x = cs_calloc (S ? S->m2 : 1, sizeof (double)) ; /* get workspace */ + ok = (AT && S && N && x) ; + if (ok) + { + cs_pvec (S->q, b, x, m) ; /* x(q(0:m-1)) = b(0:m-1) */ + cs_utsolve (N->U, x) ; /* x = R'\x */ + for (k = m-1 ; k >= 0 ; k--) /* apply Householder refl. to x */ + { + cs_happly (N->L, k, N->B [k], x) ; + } + cs_pvec (S->pinv, x, b, n) ; /* b(0:n-1) = x(p(0:n-1)) */ + } + } + cs_free (x) ; + cs_sfree (S) ; + cs_nfree (N) ; + cs_spfree (AT) ; + return (ok) ; +} diff --git a/lib/csparse/cs_randperm.c b/lib/csparse/cs_randperm.c new file mode 100644 index 0000000..de6b007 --- /dev/null +++ b/lib/csparse/cs_randperm.c @@ -0,0 +1,22 @@ +#include "cs.h" +/* return a random permutation vector, the identity perm, or p = n-1:-1:0. + * seed = -1 means p = n-1:-1:0. seed = 0 means p = identity. otherwise + * p = random permutation. */ +int *cs_randperm (int n, int seed) +{ + int *p, k, j, t ; + if (seed == 0) return (NULL) ; /* return p = NULL (identity) */ + p = cs_malloc (n, sizeof (int)) ; /* allocate result */ + if (!p) return (NULL) ; /* out of memory */ + for (k = 0 ; k < n ; k++) p [k] = n-k-1 ; + if (seed == -1) return (p) ; /* return reverse permutation */ + srand (seed) ; /* get new random number seed */ + for (k = 0 ; k < n ; k++) + { + j = k + (rand ( ) % (n-k)) ; /* j = rand int in range k to n-1 */ + t = p [j] ; /* swap p[k] and p[j] */ + p [j] = p [k] ; + p [k] = t ; + } + return (p) ; +} diff --git a/lib/csparse/cs_reach.c b/lib/csparse/cs_reach.c new file mode 100644 index 0000000..29e8ca5 --- /dev/null +++ b/lib/csparse/cs_reach.c @@ -0,0 +1,19 @@ +#include "cs.h" +/* xi [top...n-1] = nodes reachable from graph of G*P' via nodes in B(:,k). + * xi [n...2n-1] used as workspace */ +int cs_reach (cs *G, const cs *B, int k, int *xi, const int *pinv) +{ + int p, n, top, *Bp, *Bi, *Gp ; + if (!CS_CSC (G) || !CS_CSC (B) || !xi) return (-1) ; /* check inputs */ + n = G->n ; Bp = B->p ; Bi = B->i ; Gp = G->p ; + top = n ; + for (p = Bp [k] ; p < Bp [k+1] ; p++) + { + if (!CS_MARKED (Gp, Bi [p])) /* start a dfs at unmarked node i */ + { + top = cs_dfs (Bi [p], G, top, xi, xi+n, pinv) ; + } + } + for (p = top ; p < n ; p++) CS_MARK (Gp, xi [p]) ; /* restore G */ + return (top) ; +} diff --git a/lib/csparse/cs_scatter.c b/lib/csparse/cs_scatter.c new file mode 100644 index 0000000..0077826 --- /dev/null +++ b/lib/csparse/cs_scatter.c @@ -0,0 +1,22 @@ +#include "cs.h" +/* x = x + beta * A(:,j), where x is a dense vector and A(:,j) is sparse */ +int cs_scatter (const cs *A, int j, double beta, int *w, double *x, int mark, + cs *C, int nz) +{ + int i, p, *Ap, *Ai, *Ci ; + double *Ax ; + if (!CS_CSC (A) || !w || !CS_CSC (C)) return (-1) ; /* check inputs */ + Ap = A->p ; Ai = A->i ; Ax = A->x ; Ci = C->i ; + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + i = Ai [p] ; /* A(i,j) is nonzero */ + if (w [i] < mark) + { + w [i] = mark ; /* i is new entry in column j */ + Ci [nz++] = i ; /* add i to pattern of C(:,j) */ + if (x) x [i] = beta * Ax [p] ; /* x(i) = beta*A(i,j) */ + } + else if (x) x [i] += beta * Ax [p] ; /* i exists in C(:,j) already */ + } + return (nz) ; +} diff --git a/lib/csparse/cs_scc.c b/lib/csparse/cs_scc.c new file mode 100644 index 0000000..e56d086 --- /dev/null +++ b/lib/csparse/cs_scc.c @@ -0,0 +1,41 @@ +#include "cs.h" +/* find the strongly connected components of a square matrix */ +csd *cs_scc (cs *A) /* matrix A temporarily modified, then restored */ +{ + int n, i, k, b, nb = 0, top, *xi, *pstack, *p, *r, *Ap, *ATp, *rcopy, *Blk ; + cs *AT ; + csd *D ; + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + n = A->n ; Ap = A->p ; + D = cs_dalloc (n, 0) ; /* allocate result */ + AT = cs_transpose (A, 0) ; /* AT = A' */ + xi = cs_malloc (2*n+1, sizeof (int)) ; /* get workspace */ + if (!D || !AT || !xi) return (cs_ddone (D, AT, xi, 0)) ; + Blk = xi ; rcopy = pstack = xi + n ; + p = D->p ; r = D->r ; ATp = AT->p ; + top = n ; + for (i = 0 ; i < n ; i++) /* first dfs(A) to find finish times (xi) */ + { + if (!CS_MARKED (Ap, i)) top = cs_dfs (i, A, top, xi, pstack, NULL) ; + } + for (i = 0 ; i < n ; i++) CS_MARK (Ap, i) ; /* restore A; unmark all nodes*/ + top = n ; + nb = n ; + for (k = 0 ; k < n ; k++) /* dfs(A') to find strongly connnected comp */ + { + i = xi [k] ; /* get i in reverse order of finish times */ + if (CS_MARKED (ATp, i)) continue ; /* skip node i if already ordered */ + r [nb--] = top ; /* node i is the start of a component in p */ + top = cs_dfs (i, AT, top, p, pstack, NULL) ; + } + r [nb] = 0 ; /* first block starts at zero; shift r up */ + for (k = nb ; k <= n ; k++) r [k-nb] = r [k] ; + D->nb = nb = n-nb ; /* nb = # of strongly connected components */ + for (b = 0 ; b < nb ; b++) /* sort each block in natural order */ + { + for (k = r [b] ; k < r [b+1] ; k++) Blk [p [k]] = b ; + } + for (b = 0 ; b <= nb ; b++) rcopy [b] = r [b] ; + for (i = 0 ; i < n ; i++) p [rcopy [Blk [i]]++] = i ; + return (cs_ddone (D, AT, xi, 1)) ; +} diff --git a/lib/csparse/cs_schol.c b/lib/csparse/cs_schol.c new file mode 100644 index 0000000..32c6734 --- /dev/null +++ b/lib/csparse/cs_schol.c @@ -0,0 +1,26 @@ +#include "cs.h" +/* ordering and symbolic analysis for a Cholesky factorization */ +css *cs_schol (int order, const cs *A) +{ + int n, *c, *post, *P ; + cs *C ; + css *S ; + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + n = A->n ; + S = cs_calloc (1, sizeof (css)) ; /* allocate result S */ + if (!S) return (NULL) ; /* out of memory */ + P = cs_amd (order, A) ; /* P = amd(A+A'), or natural */ + S->pinv = cs_pinv (P, n) ; /* find inverse permutation */ + cs_free (P) ; + if (order && !S->pinv) return (cs_sfree (S)) ; + C = cs_symperm (A, S->pinv, 0) ; /* C = spones(triu(A(P,P))) */ + S->parent = cs_etree (C, 0) ; /* find etree of C */ + post = cs_post (S->parent, n) ; /* postorder the etree */ + c = cs_counts (C, S->parent, post, 0) ; /* find column counts of chol(C) */ + cs_free (post) ; + cs_spfree (C) ; + S->cp = cs_malloc (n+1, sizeof (int)) ; /* allocate result S->cp */ + S->unz = S->lnz = cs_cumsum (S->cp, c, n) ; /* find column pointers for L */ + cs_free (c) ; + return ((S->lnz >= 0) ? S : cs_sfree (S)) ; +} diff --git a/lib/csparse/cs_spsolve.c b/lib/csparse/cs_spsolve.c new file mode 100644 index 0000000..d004f9c --- /dev/null +++ b/lib/csparse/cs_spsolve.c @@ -0,0 +1,28 @@ +#include "cs.h" +/* solve Gx=b(:,k), where G is either upper (lo=0) or lower (lo=1) triangular */ +int cs_spsolve (cs *G, const cs *B, int k, int *xi, double *x, const int *pinv, + int lo) +{ + int j, J, p, q, px, top, n, *Gp, *Gi, *Bp, *Bi ; + double *Gx, *Bx ; + if (!CS_CSC (G) || !CS_CSC (B) || !xi || !x) return (-1) ; + Gp = G->p ; Gi = G->i ; Gx = G->x ; n = G->n ; + Bp = B->p ; Bi = B->i ; Bx = B->x ; + top = cs_reach (G, B, k, xi, pinv) ; /* xi[top..n-1]=Reach(B(:,k)) */ + for (p = top ; p < n ; p++) x [xi [p]] = 0 ; /* clear x */ + for (p = Bp [k] ; p < Bp [k+1] ; p++) x [Bi [p]] = Bx [p] ; /* scatter B */ + for (px = top ; px < n ; px++) + { + j = xi [px] ; /* x(j) is nonzero */ + J = pinv ? (pinv [j]) : j ; /* j maps to col J of G */ + if (J < 0) continue ; /* column J is empty */ + x [j] /= Gx [lo ? (Gp [J]) : (Gp [J+1]-1)] ;/* x(j) /= G(j,j) */ + p = lo ? (Gp [J]+1) : (Gp [J]) ; /* lo: L(j,j) 1st entry */ + q = lo ? (Gp [J+1]) : (Gp [J+1]-1) ; /* up: U(j,j) last entry */ + for ( ; p < q ; p++) + { + x [Gi [p]] -= Gx [p] * x [j] ; /* x(i) -= G(i,j) * x(j) */ + } + } + return (top) ; /* return top of stack */ +} diff --git a/lib/csparse/cs_sqr.c b/lib/csparse/cs_sqr.c new file mode 100644 index 0000000..ad39bd2 --- /dev/null +++ b/lib/csparse/cs_sqr.c @@ -0,0 +1,88 @@ +#include "cs.h" +/* compute nnz(V) = S->lnz, S->pinv, S->leftmost, S->m2 from A and S->parent */ +static int cs_vcount (const cs *A, css *S) +{ + int i, k, p, pa, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i, *next, *head, + *tail, *nque, *pinv, *leftmost, *w, *parent = S->parent ; + S->pinv = pinv = cs_malloc (m+n, sizeof (int)) ; /* allocate pinv, */ + S->leftmost = leftmost = cs_malloc (m, sizeof (int)) ; /* and leftmost */ + w = cs_malloc (m+3*n, sizeof (int)) ; /* get workspace */ + if (!pinv || !w || !leftmost) + { + cs_free (w) ; /* pinv and leftmost freed later */ + return (0) ; /* out of memory */ + } + next = w ; head = w + m ; tail = w + m + n ; nque = w + m + 2*n ; + for (k = 0 ; k < n ; k++) head [k] = -1 ; /* queue k is empty */ + for (k = 0 ; k < n ; k++) tail [k] = -1 ; + for (k = 0 ; k < n ; k++) nque [k] = 0 ; + for (i = 0 ; i < m ; i++) leftmost [i] = -1 ; + for (k = n-1 ; k >= 0 ; k--) + { + for (p = Ap [k] ; p < Ap [k+1] ; p++) + { + leftmost [Ai [p]] = k ; /* leftmost[i] = min(find(A(i,:)))*/ + } + } + for (i = m-1 ; i >= 0 ; i--) /* scan rows in reverse order */ + { + pinv [i] = -1 ; /* row i is not yet ordered */ + k = leftmost [i] ; + if (k == -1) continue ; /* row i is empty */ + if (nque [k]++ == 0) tail [k] = i ; /* first row in queue k */ + next [i] = head [k] ; /* put i at head of queue k */ + head [k] = i ; + } + S->lnz = 0 ; + S->m2 = m ; + for (k = 0 ; k < n ; k++) /* find row permutation and nnz(V)*/ + { + i = head [k] ; /* remove row i from queue k */ + S->lnz++ ; /* count V(k,k) as nonzero */ + if (i < 0) i = S->m2++ ; /* add a fictitious row */ + pinv [i] = k ; /* associate row i with V(:,k) */ + if (--nque [k] <= 0) continue ; /* skip if V(k+1:m,k) is empty */ + S->lnz += nque [k] ; /* nque [k] is nnz (V(k+1:m,k)) */ + if ((pa = parent [k]) != -1) /* move all rows to parent of k */ + { + if (nque [pa] == 0) tail [pa] = tail [k] ; + next [tail [k]] = head [pa] ; + head [pa] = next [i] ; + nque [pa] += nque [k] ; + } + } + for (i = 0 ; i < m ; i++) if (pinv [i] < 0) pinv [i] = k++ ; + cs_free (w) ; + return (1) ; +} + +/* symbolic ordering and analysis for QR or LU */ +css *cs_sqr (int order, const cs *A, int qr) +{ + int n, k, ok = 1, *post ; + css *S ; + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + n = A->n ; + S = cs_calloc (1, sizeof (css)) ; /* allocate result S */ + if (!S) return (NULL) ; /* out of memory */ + S->q = cs_amd (order, A) ; /* fill-reducing ordering */ + if (order && !S->q) return (cs_sfree (S)) ; + if (qr) /* QR symbolic analysis */ + { + cs *C = order ? cs_permute (A, NULL, S->q, 0) : ((cs *) A) ; + S->parent = cs_etree (C, 1) ; /* etree of C'*C, where C=A(:,q) */ + post = cs_post (S->parent, n) ; + S->cp = cs_counts (C, S->parent, post, 1) ; /* col counts chol(C'*C) */ + cs_free (post) ; + ok = C && S->parent && S->cp && cs_vcount (C, S) ; + if (ok) for (S->unz = 0, k = 0 ; k < n ; k++) S->unz += S->cp [k] ; + ok = ok && S->lnz >= 0 && S->unz >= 0 ; /* int overflow guard */ + if (order) cs_spfree (C) ; + } + else + { + S->unz = 4*(A->p [n]) + n ; /* for LU factorization only, */ + S->lnz = S->unz ; /* guess nnz(L) and nnz(U) */ + } + return (ok ? S : cs_sfree (S)) ; /* return result S */ +} diff --git a/lib/csparse/cs_symperm.c b/lib/csparse/cs_symperm.c new file mode 100644 index 0000000..668526d --- /dev/null +++ b/lib/csparse/cs_symperm.c @@ -0,0 +1,39 @@ +#include "cs.h" +/* C = A(p,p) where A and C are symmetric the upper part stored; pinv not p */ +cs *cs_symperm (const cs *A, const int *pinv, int values) +{ + int i, j, p, q, i2, j2, n, *Ap, *Ai, *Cp, *Ci, *w ; + double *Cx, *Ax ; + cs *C ; + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ; + C = cs_spalloc (n, n, Ap [n], values && (Ax != NULL), 0) ; /* alloc result*/ + w = cs_calloc (n, sizeof (int)) ; /* get workspace */ + if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */ + Cp = C->p ; Ci = C->i ; Cx = C->x ; + for (j = 0 ; j < n ; j++) /* count entries in each column of C */ + { + j2 = pinv ? pinv [j] : j ; /* column j of A is column j2 of C */ + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + i = Ai [p] ; + if (i > j) continue ; /* skip lower triangular part of A */ + i2 = pinv ? pinv [i] : i ; /* row i of A is row i2 of C */ + w [CS_MAX (i2, j2)]++ ; /* column count of C */ + } + } + cs_cumsum (Cp, w, n) ; /* compute column pointers of C */ + for (j = 0 ; j < n ; j++) + { + j2 = pinv ? pinv [j] : j ; /* column j of A is column j2 of C */ + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + i = Ai [p] ; + if (i > j) continue ; /* skip lower triangular part of A*/ + i2 = pinv ? pinv [i] : i ; /* row i of A is row i2 of C */ + Ci [q = w [CS_MAX (i2, j2)]++] = CS_MIN (i2, j2) ; + if (Cx) Cx [q] = Ax [p] ; + } + } + return (cs_done (C, w, NULL, 1)) ; /* success; free workspace, return C */ +} diff --git a/lib/csparse/cs_tdfs.c b/lib/csparse/cs_tdfs.c new file mode 100644 index 0000000..ad19890 --- /dev/null +++ b/lib/csparse/cs_tdfs.c @@ -0,0 +1,24 @@ +#include "cs.h" +/* depth-first search and postorder of a tree rooted at node j */ +int cs_tdfs (int j, int k, int *head, const int *next, int *post, int *stack) +{ + int i, p, top = 0 ; + if (!head || !next || !post || !stack) return (-1) ; /* check inputs */ + stack [0] = j ; /* place j on the stack */ + while (top >= 0) /* while (stack is not empty) */ + { + p = stack [top] ; /* p = top of stack */ + i = head [p] ; /* i = youngest child of p */ + if (i == -1) + { + top-- ; /* p has no unordered children left */ + post [k++] = p ; /* node p is the kth postordered node */ + } + else + { + head [p] = next [i] ; /* remove i from children of p */ + stack [++top] = i ; /* start dfs on child node i */ + } + } + return (k) ; +} diff --git a/lib/csparse/cs_transpose.c b/lib/csparse/cs_transpose.c new file mode 100644 index 0000000..9daa84d --- /dev/null +++ b/lib/csparse/cs_transpose.c @@ -0,0 +1,25 @@ +#include "cs.h" +/* C = A' */ +cs *cs_transpose (const cs *A, int values) +{ + int p, q, j, *Cp, *Ci, n, m, *Ap, *Ai, *w ; + double *Cx, *Ax ; + cs *C ; + if (!CS_CSC (A)) return (NULL) ; /* check inputs */ + m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ; + C = cs_spalloc (n, m, Ap [n], values && Ax, 0) ; /* allocate result */ + w = cs_calloc (m, sizeof (int)) ; /* get workspace */ + if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */ + Cp = C->p ; Ci = C->i ; Cx = C->x ; + for (p = 0 ; p < Ap [n] ; p++) w [Ai [p]]++ ; /* row counts */ + cs_cumsum (Cp, w, m) ; /* row pointers */ + for (j = 0 ; j < n ; j++) + { + for (p = Ap [j] ; p < Ap [j+1] ; p++) + { + Ci [q = w [Ai [p]]++] = j ; /* place A(i,j) as entry C(j,i) */ + if (Cx) Cx [q] = Ax [p] ; + } + } + return (cs_done (C, w, NULL, 1)) ; /* success; free w and return C */ +} diff --git a/lib/csparse/cs_updown.c b/lib/csparse/cs_updown.c new file mode 100644 index 0000000..90dc4cc --- /dev/null +++ b/lib/csparse/cs_updown.c @@ -0,0 +1,37 @@ +#include "cs.h" +/* sparse Cholesky update/downdate, L*L' + sigma*w*w' (sigma = +1 or -1) */ +int cs_updown (cs *L, int sigma, const cs *C, const int *parent) +{ + int n, p, f, j, *Lp, *Li, *Cp, *Ci ; + double *Lx, *Cx, alpha, beta = 1, delta, gamma, w1, w2, *w, beta2 = 1 ; + if (!CS_CSC (L) || !CS_CSC (C) || !parent) return (0) ; /* check inputs */ + Lp = L->p ; Li = L->i ; Lx = L->x ; n = L->n ; + Cp = C->p ; Ci = C->i ; Cx = C->x ; + if ((p = Cp [0]) >= Cp [1]) return (1) ; /* return if C empty */ + w = cs_malloc (n, sizeof (double)) ; /* get workspace */ + if (!w) return (0) ; /* out of memory */ + f = Ci [p] ; + for ( ; p < Cp [1] ; p++) f = CS_MIN (f, Ci [p]) ; /* f = min (find (C)) */ + for (j = f ; j != -1 ; j = parent [j]) w [j] = 0 ; /* clear workspace w */ + for (p = Cp [0] ; p < Cp [1] ; p++) w [Ci [p]] = Cx [p] ; /* w = C */ + for (j = f ; j != -1 ; j = parent [j]) /* walk path f up to root */ + { + p = Lp [j] ; + alpha = w [j] / Lx [p] ; /* alpha = w(j) / L(j,j) */ + beta2 = beta*beta + sigma*alpha*alpha ; + if (beta2 <= 0) break ; /* not positive definite */ + beta2 = sqrt (beta2) ; + delta = (sigma > 0) ? (beta / beta2) : (beta2 / beta) ; + gamma = sigma * alpha / (beta2 * beta) ; + Lx [p] = delta * Lx [p] + ((sigma > 0) ? (gamma * w [j]) : 0) ; + beta = beta2 ; + for (p++ ; p < Lp [j+1] ; p++) + { + w1 = w [Li [p]] ; + w [Li [p]] = w2 = w1 - alpha * Lx [p] ; + Lx [p] = delta * Lx [p] + gamma * ((sigma > 0) ? w1 : w2) ; + } + } + cs_free (w) ; + return (beta2 > 0) ; +} diff --git a/lib/csparse/cs_usolve.c b/lib/csparse/cs_usolve.c new file mode 100644 index 0000000..40e0892 --- /dev/null +++ b/lib/csparse/cs_usolve.c @@ -0,0 +1,18 @@ +#include "cs.h" +/* solve Ux=b where x and b are dense. x=b on input, solution on output. */ +int cs_usolve (const cs *U, double *x) +{ + int p, j, n, *Up, *Ui ; + double *Ux ; + if (!CS_CSC (U) || !x) return (0) ; /* check inputs */ + n = U->n ; Up = U->p ; Ui = U->i ; Ux = U->x ; + for (j = n-1 ; j >= 0 ; j--) + { + x [j] /= Ux [Up [j+1]-1] ; + for (p = Up [j] ; p < Up [j+1]-1 ; p++) + { + x [Ui [p]] -= Ux [p] * x [j] ; + } + } + return (1) ; +} diff --git a/lib/csparse/cs_util.c b/lib/csparse/cs_util.c new file mode 100644 index 0000000..5684cd2 --- /dev/null +++ b/lib/csparse/cs_util.c @@ -0,0 +1,119 @@ +#include "cs.h" +/* allocate a sparse matrix (triplet form or compressed-column form) */ +cs *cs_spalloc (int m, int n, int nzmax, int values, int triplet) +{ + cs *A = cs_calloc (1, sizeof (cs)) ; /* allocate the cs struct */ + if (!A) return (NULL) ; /* out of memory */ + A->m = m ; /* define dimensions and nzmax */ + A->n = n ; + A->nzmax = nzmax = CS_MAX (nzmax, 1) ; + A->nz = triplet ? 0 : -1 ; /* allocate triplet or comp.col */ + A->p = cs_malloc (triplet ? nzmax : n+1, sizeof (int)) ; + A->i = cs_malloc (nzmax, sizeof (int)) ; + A->x = values ? cs_malloc (nzmax, sizeof (double)) : NULL ; + return ((!A->p || !A->i || (values && !A->x)) ? cs_spfree (A) : A) ; +} + +/* change the max # of entries sparse matrix */ +int cs_sprealloc (cs *A, int nzmax) +{ + int ok, oki, okj = 1, okx = 1 ; + if (!A) return (0) ; + if (nzmax <= 0) nzmax = (CS_CSC (A)) ? (A->p [A->n]) : A->nz ; + A->i = cs_realloc (A->i, nzmax, sizeof (int), &oki) ; + if (CS_TRIPLET (A)) A->p = cs_realloc (A->p, nzmax, sizeof (int), &okj) ; + if (A->x) A->x = cs_realloc (A->x, nzmax, sizeof (double), &okx) ; + ok = (oki && okj && okx) ; + if (ok) A->nzmax = nzmax ; + return (ok) ; +} + +/* free a sparse matrix */ +cs *cs_spfree (cs *A) +{ + if (!A) return (NULL) ; /* do nothing if A already NULL */ + cs_free (A->p) ; + cs_free (A->i) ; + cs_free (A->x) ; + return (cs_free (A)) ; /* free the cs struct and return NULL */ +} + +/* free a numeric factorization */ +csn *cs_nfree (csn *N) +{ + if (!N) return (NULL) ; /* do nothing if N already NULL */ + cs_spfree (N->L) ; + cs_spfree (N->U) ; + cs_free (N->pinv) ; + cs_free (N->B) ; + return (cs_free (N)) ; /* free the csn struct and return NULL */ +} + +/* free a symbolic factorization */ +css *cs_sfree (css *S) +{ + if (!S) return (NULL) ; /* do nothing if S already NULL */ + cs_free (S->pinv) ; + cs_free (S->q) ; + cs_free (S->parent) ; + cs_free (S->cp) ; + cs_free (S->leftmost) ; + return (cs_free (S)) ; /* free the css struct and return NULL */ +} + +/* allocate a cs_dmperm or cs_scc result */ +csd *cs_dalloc (int m, int n) +{ + csd *D ; + D = cs_calloc (1, sizeof (csd)) ; + if (!D) return (NULL) ; + D->p = cs_malloc (m, sizeof (int)) ; + D->r = cs_malloc (m+6, sizeof (int)) ; + D->q = cs_malloc (n, sizeof (int)) ; + D->s = cs_malloc (n+6, sizeof (int)) ; + return ((!D->p || !D->r || !D->q || !D->s) ? cs_dfree (D) : D) ; +} + +/* free a cs_dmperm or cs_scc result */ +csd *cs_dfree (csd *D) +{ + if (!D) return (NULL) ; /* do nothing if D already NULL */ + cs_free (D->p) ; + cs_free (D->q) ; + cs_free (D->r) ; + cs_free (D->s) ; + return (cs_free (D)) ; +} + +/* free workspace and return a sparse matrix result */ +cs *cs_done (cs *C, void *w, void *x, int ok) +{ + cs_free (w) ; /* free workspace */ + cs_free (x) ; + return (ok ? C : cs_spfree (C)) ; /* return result if OK, else free it */ +} + +/* free workspace and return int array result */ +int *cs_idone (int *p, cs *C, void *w, int ok) +{ + cs_spfree (C) ; /* free temporary matrix */ + cs_free (w) ; /* free workspace */ + return (ok ? p : cs_free (p)) ; /* return result if OK, else free it */ +} + +/* free workspace and return a numeric factorization (Cholesky, LU, or QR) */ +csn *cs_ndone (csn *N, cs *C, void *w, void *x, int ok) +{ + cs_spfree (C) ; /* free temporary matrix */ + cs_free (w) ; /* free workspace */ + cs_free (x) ; + return (ok ? N : cs_nfree (N)) ; /* return result if OK, else free it */ +} + +/* free workspace and return a csd result */ +csd *cs_ddone (csd *D, cs *C, void *w, int ok) +{ + cs_spfree (C) ; /* free temporary matrix */ + cs_free (w) ; /* free workspace */ + return (ok ? D : cs_dfree (D)) ; /* return result if OK, else free it */ +} diff --git a/lib/csparse/cs_utsolve.c b/lib/csparse/cs_utsolve.c new file mode 100644 index 0000000..9d6bf63 --- /dev/null +++ b/lib/csparse/cs_utsolve.c @@ -0,0 +1,18 @@ +#include "cs.h" +/* solve U'x=b where x and b are dense. x=b on input, solution on output. */ +int cs_utsolve (const cs *U, double *x) +{ + int p, j, n, *Up, *Ui ; + double *Ux ; + if (!CS_CSC (U) || !x) return (0) ; /* check inputs */ + n = U->n ; Up = U->p ; Ui = U->i ; Ux = U->x ; + for (j = 0 ; j < n ; j++) + { + for (p = Up [j] ; p < Up [j+1]-1 ; p++) + { + x [j] -= Ux [p] * x [Ui [p]] ; + } + x [j] /= Ux [Up [j+1]-1] ; + } + return (1) ; +}