Compare commits

...

68 commits

Author SHA1 Message Date
f8380fa246 add hurler action 2024-07-04 22:43:19 -05:00
1a377855ce implement full simulation with epoch save/restore 2024-07-04 18:40:34 -05:00
005a4cf24c implement all senses/actions 2024-07-04 07:36:49 -05:00
17422cad2a make web version 2024-07-04 07:03:08 -05:00
12392b5d91 create and destroy organisms 2024-06-05 09:42:07 -05:00
caefc84110 add genome copying 2024-06-04 16:41:40 -05:00
143aa32082 update lily-test 2024-06-04 16:35:08 -05:00
ccec744e78 implement sensing and acting 2024-06-04 12:40:05 -05:00
d56f099118 implement the jump instructions 2024-06-04 12:27:15 -05:00
eb660e1b10 implement loading instructions 2024-06-03 23:54:09 -05:00
00c193a3a9 implement storage instructions 2024-06-03 23:48:16 -05:00
9d397ce883 redesign instructions and implement bitwise operators 2024-06-03 23:02:31 -05:00
a84002b203 implement basic arithmetic 2024-06-03 19:40:57 -05:00
86fcf0abee copy genome to memory correctly 2024-06-03 17:22:47 -05:00
f6957d8a77 create/destroy minds 2024-06-03 17:15:35 -05:00
b16ae0c5ef switch to assembly-based minds 2024-06-03 11:25:41 -05:00
6288cdf903 fix memory leaks 2024-03-10 14:38:25 -05:00
defd656e41 fix bounding size bug on nr_mind_create 2024-03-08 20:57:25 -06:00
37d100ec64 fix OOR errors on adding/removing neurons 2024-03-08 20:52:28 -06:00
af54d42411 fix OOR error on shift/zero weights 2024-03-08 20:47:15 -06:00
5be6591a46 implement nr_mind_compute 2024-03-08 20:37:14 -06:00
c5be66eaa7 switch to csparse 3.1.3 and begin implementing nr_mind_t functions 2024-03-08 18:25:59 -06:00
ef90da1f75 add duplication and cs_sparse conversion functions for nr_sparse_t 2024-03-07 14:39:01 -06:00
7adc0874ea refactor: nr_isparse_* -> nr_sparse_* 2024-02-18 18:27:50 -06:00
7043ca9f62 implement nr_isparse_remove_neuron 2024-02-18 17:11:32 -06:00
41bf6e3ea7 implement nr_isparse_add_neuron 2024-02-17 12:31:45 -06:00
f96b7c1902 implement DOK nr_isparse_zero_weight 2024-02-17 03:20:15 -06:00
8a845e9bea begin switch to DOK isparse matrices 2024-02-17 01:42:27 -06:00
1a7bde7333 implement nr_isparse_zero_weight 2024-02-16 21:53:29 -06:00
724335883b implement isparse_shift_weight 2024-02-16 20:42:55 -06:00
fc7ad069a6 add csparse library 2024-02-15 15:22:10 -06:00
0317fb149c implement xorshift prng 2024-02-15 13:44:28 -06:00
af6553a000 remove haskell files 2024-02-15 11:01:09 -06:00
6167fb4a23 implement wrapping queries 2024-01-06 02:19:14 -06:00
cecfa5402f implement insertion and querying 2024-01-06 02:06:51 -06:00
d61721be17 remove Space.adjacent 2024-01-05 12:32:51 -06:00
4ede59d49b implement spatial distance 2024-01-05 11:42:23 -06:00
7f2ce9e6e0 implement spatial adjacency 2024-01-05 01:42:03 -06:00
c6fd06f123 clean up project structure 2024-01-04 22:35:09 -06:00
303bf3090e begin writing agent tests 2024-01-04 22:24:51 -06:00
b2a8a61b28 prevent applying invalid proposals 2023-12-08 14:21:19 -06:00
15a0b3d465 begin implementing applyLatticeProposals 2023-12-05 12:50:56 -06:00
899ad0ed13 refactor: Lattice -> World.Lattice 2023-12-04 20:10:18 -06:00
749df30696 implement Combinable for Proposal 2023-12-04 11:55:55 -06:00
1a3814b5ea implement Combinable for AgentPropList 2023-12-04 11:44:09 -06:00
135ce23bd1 implement Combinable for LatticePropList 2023-12-03 23:44:33 -06:00
450ac00012 begin implementing lattice proposal list merging 2023-11-30 17:34:58 -06:00
671c632d5c add Compatible typeclass 2023-11-30 14:19:20 -06:00
bf40a269da begin adding world types 2023-11-30 13:54:44 -06:00
587cbadb3e implement updateLattice 2023-11-29 20:48:46 -06:00
4facdd8a1a implement parseGenome 2023-11-29 17:40:38 -06:00
98eabf2989 add mutateGenome 2023-11-29 17:11:07 -06:00
d28e8e243d implement mutateGenomeRemoveGene 2023-11-29 16:45:14 -06:00
4c63f3d4db implement mutateGenomeAddGene 2023-11-29 16:25:38 -06:00
3047e46f34 imlement mutateGenomeRemoveNeuron 2023-11-29 15:45:44 -06:00
62a04f27ec implement mutateGenomeAddInternal 2023-11-29 12:41:09 -06:00
fa28b97de1 rewrite genome tests to use approximate float comparison 2023-11-29 12:12:46 -06:00
a506d45e16 implement individual gene mutation functions 2023-11-29 02:59:51 -06:00
e8d7a5237c show error messages on bad input/state 2023-11-22 10:39:33 -06:00
8fb358e847 fix warnings and bug in connectNeurons 2023-11-21 23:48:32 -06:00
435c52c733 implement self-connected computation 2023-11-21 23:41:38 -06:00
26be9f9dc7 implement first-pass internal neurons 2023-11-21 21:23:56 -06:00
a005bcdf7b implement output computation 2023-11-21 21:06:20 -06:00
fae22e6828 add failing compute tests 2023-11-21 14:16:38 -06:00
23e8acfc1e add connectNeurons 2023-11-21 11:34:48 -06:00
6ab486c73d add createEmptyNetwork 2023-11-21 10:35:07 -06:00
d128c42a96 initial stack configuration 2023-11-21 10:11:19 -06:00
f9fc4d26ec clear out js files 2023-11-16 14:50:00 -06:00
141 changed files with 6852 additions and 6015 deletions

View file

8
.gitignore vendored
View file

@ -1,5 +1,5 @@
yarn-error.log
node_modules/
.*
*.swp
*~
*.swp
.stack-work
*.swo
build

3
.gitmodules vendored Normal file
View file

@ -0,0 +1,3 @@
[submodule "deps/lily-test"]
path = deps/lily-test
url = https://sanine.net/git/lily-test

47
CMakeLists.txt Normal file
View file

@ -0,0 +1,47 @@
cmake_minimum_required(VERSION 3.8)
project(nerine)
set(DEPS_DIR ${CMAKE_CURRENT_SOURCE_DIR}/deps)
set(NERINE_SRC ${CMAKE_CURRENT_SOURCE_DIR}/src)
set(TEST_SRC ${CMAKE_CURRENT_SOURCE_DIR}/test)
# add_subdirectory(${DEPS_DIR}/csparse)
add_library(nerine
${NERINE_SRC}/nerine/rand.c
${NERINE_SRC}/nerine/genome.c
${NERINE_SRC}/nerine/mind.c
${NERINE_SRC}/nerine/organism.c
)
target_include_directories(nerine PUBLIC
${NERINE_SRC}
# ${DEPS_DIR}/csparse/Include
)
# target_link_libraries(nerine PUBLIC
# csparse
# )
add_executable(test-nerine
${TEST_SRC}/main.c
${TEST_SRC}/nerine/genome.c
${TEST_SRC}/nerine/mind.c
${TEST_SRC}/nerine/organism.c
)
target_include_directories(test-nerine PUBLIC
${TEST_SRC}
${NERINE_SRC}
${DEPS_DIR}/lily-test/
)
# clean up test filenames c:
string(LENGTH "${CMAKE_SOURCE_DIR}/" SOURCE_PATH_SIZE)
target_compile_definitions(test-nerine PUBLIC
SOURCE_PATH_SIZE=${SOURCE_PATH_SIZE}
)
set_target_properties(test-nerine PROPERTIES
C_STANDARD 99
CMAKE_C_FLAGS "-Wall -Wextra -Werror -Wfatal-errors -Wpedantic -pedantic-errors")
target_link_libraries(test-nerine PUBLIC nerine)

1
README.md Normal file
View file

@ -0,0 +1 @@
# nerine

61
deps/csparse/CMakeLists.txt vendored Normal file
View file

@ -0,0 +1,61 @@
cmake_minimum_required(VERSION 3.8)
project(CSPARSE)
set(SRC ${CMAKE_CURRENT_SOURCE_DIR}/Source)
add_library(csparse SHARED
${SRC}/cs_add.c
${SRC}/cs_dupl.c
${SRC}/cs_lsolve.c
${SRC}/cs_print.c
${SRC}/cs_symperm.c
${SRC}/cs_amd.c
${SRC}/cs_entry.c
${SRC}/cs_ltsolve.c
${SRC}/cs_pvec.c
${SRC}/cs_tdfs.c
${SRC}/cs_chol.c
${SRC}/cs_ereach.c
${SRC}/cs_lu.c
${SRC}/cs_qr.c
${SRC}/cs_transpose.c
${SRC}/cs_cholsol.c
${SRC}/cs_etree.c
${SRC}/cs_lusol.c
${SRC}/cs_qrsol.c
${SRC}/cs_updown.c
${SRC}/cs_compress.c
${SRC}/cs_fkeep.c
${SRC}/cs_malloc.c
${SRC}/cs_randperm.c
${SRC}/cs_usolve.c
${SRC}/cs_counts.c
${SRC}/cs_gaxpy.c
${SRC}/cs_maxtrans.c
${SRC}/cs_reach.c
${SRC}/cs_util.c
${SRC}/cs_cumsum.c
${SRC}/cs_happly.c
${SRC}/cs_multiply.c
${SRC}/cs_scatter.c
${SRC}/cs_utsolve.c
${SRC}/cs_dfs.c
${SRC}/cs_house.c
${SRC}/cs_norm.c
${SRC}/cs_scc.c
${SRC}/cs_dmperm.c
${SRC}/cs_ipvec.c
${SRC}/cs_permute.c
${SRC}/cs_schol.c
${SRC}/cs_droptol.c
${SRC}/cs_leaf.c
${SRC}/cs_pinv.c
${SRC}/cs_spsolve.c
${SRC}/cs_dropzeros.c
${SRC}/cs_load.c
${SRC}/cs_post.c
${SRC}/cs_sqr.c
)
target_include_directories(csparse PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/Include)
target_link_libraries(csparse PUBLIC m)

37
deps/csparse/Demo/Makefile vendored Normal file
View file

@ -0,0 +1,37 @@
CF = $(CFLAGS) $(CPPFLAGS) $(TARGET_ARCH) -O
I = -I../Include
CS = ../Lib/libcsparse.a
all: lib cs_demo1 cs_demo2 cs_demo3
- ./cs_demo1 < ../Matrix/t1
- ./cs_demo2 < ../Matrix/t1
- ./cs_demo2 < ../Matrix/ash219
- ./cs_demo2 < ../Matrix/bcsstk01
- ./cs_demo2 < ../Matrix/fs_183_1
- ./cs_demo2 < ../Matrix/mbeacxc
- ./cs_demo2 < ../Matrix/west0067
- ./cs_demo2 < ../Matrix/lp_afiro
- ./cs_demo2 < ../Matrix/bcsstk16
- ./cs_demo3 < ../Matrix/bcsstk01
- ./cs_demo3 < ../Matrix/bcsstk16
lib:
( cd ../Lib ; $(MAKE) )
cs_demo1: lib cs_demo1.c Makefile
$(CC) $(CF) $(I) -o cs_demo1 cs_demo1.c $(CS) -lm
cs_demo2: lib cs_demo2.c cs_demo.c cs_demo.h Makefile
$(CC) $(CF) $(I) -o cs_demo2 cs_demo2.c cs_demo.c $(CS) -lm
cs_demo3: lib cs_demo3.c cs_demo.c cs_demo.h Makefile
$(CC) $(CF) $(I) -o cs_demo3 cs_demo3.c cs_demo.c $(CS) -lm
clean:
- $(RM) *.o
purge: distclean
distclean: clean
- $(RM) -r cs_demo1 cs_demo2 cs_demo3 *.a *.dSYM *.obj *.dll

6
deps/csparse/Demo/README.txt vendored Normal file
View file

@ -0,0 +1,6 @@
CSparse/Demo: to compile a run the demos, just type "make" in this directory.
The printed residuals should all be small, except for the mbeacxc matrix
(which is numerically and structurally singular), and ash219 (which is a
least-squares problem). See cs_demo.out for the proper output of "make".
Timothy A. Davis, http://www.suitesparse.com

289
deps/csparse/Demo/cs_demo.c vendored Normal file
View file

@ -0,0 +1,289 @@
#include "cs_demo.h"
#include <time.h>
/* 1 if A is square & upper tri., -1 if square & lower tri., 0 otherwise */
static csi is_sym (cs *A)
{
csi is_upper, is_lower, j, p, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i ;
if (m != n) return (0) ;
is_upper = 1 ;
is_lower = 1 ;
for (j = 0 ; j < n ; j++)
{
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
if (Ai [p] > j) is_upper = 0 ;
if (Ai [p] < j) is_lower = 0 ;
}
}
return (is_upper ? 1 : (is_lower ? -1 : 0)) ;
}
/* true for off-diagonal entries */
static csi dropdiag (csi i, csi j, double aij, void *other) { return (i != j) ;}
/* C = A + triu(A,1)' */
static cs *make_sym (cs *A)
{
cs *AT, *C ;
AT = cs_transpose (A, 1) ; /* AT = A' */
cs_fkeep (AT, &dropdiag, NULL) ; /* drop diagonal entries from AT */
C = cs_add (A, AT, 1, 1) ; /* C = A+AT */
cs_spfree (AT) ;
return (C) ;
}
/* create a right-hand side */
static void rhs (double *x, double *b, csi m)
{
csi i ;
for (i = 0 ; i < m ; i++) b [i] = 1 + ((double) i) / m ;
for (i = 0 ; i < m ; i++) x [i] = b [i] ;
}
/* infinity-norm of x */
static double norm (double *x, csi n)
{
csi i ;
double normx = 0 ;
for (i = 0 ; i < n ; i++) normx = CS_MAX (normx, fabs (x [i])) ;
return (normx) ;
}
/* compute residual, norm(A*x-b,inf) / (norm(A,1)*norm(x,inf) + norm(b,inf)) */
static void print_resid (csi ok, cs *A, double *x, double *b, double *resid)
{
csi i, m, n ;
if (!ok) { printf (" (failed)\n") ; return ; }
m = A->m ; n = A->n ;
for (i = 0 ; i < m ; i++) resid [i] = -b [i] ; /* resid = -b */
cs_gaxpy (A, x, resid) ; /* resid = resid + A*x */
printf ("resid: %8.2e\n", norm (resid,m) / ((n == 0) ? 1 :
(cs_norm (A) * norm (x,n) + norm (b,m)))) ;
}
static double tic (void) { return (clock () / (double) CLOCKS_PER_SEC) ; }
static double toc (double t) { double s = tic () ; return (CS_MAX (0, s-t)) ; }
static void print_order (csi order)
{
switch (order)
{
case 0: printf ("natural ") ; break ;
case 1: printf ("amd(A+A') ") ; break ;
case 2: printf ("amd(S'*S) ") ; break ;
case 3: printf ("amd(A'*A) ") ; break ;
}
}
/* read a problem from a file; use %g for integers to avoid csi conflicts */
problem *get_problem (FILE *f, double tol)
{
cs *T, *A, *C ;
csi sym, m, n, mn, nz1, nz2 ;
problem *Prob ;
Prob = cs_calloc (1, sizeof (problem)) ;
if (!Prob) return (NULL) ;
T = cs_load (f) ; /* load triplet matrix T from a file */
Prob->A = A = cs_compress (T) ; /* A = compressed-column form of T */
cs_spfree (T) ; /* clear T */
if (!cs_dupl (A)) return (free_problem (Prob)) ; /* sum up duplicates */
Prob->sym = sym = is_sym (A) ; /* determine if A is symmetric */
m = A->m ; n = A->n ;
mn = CS_MAX (m,n) ;
nz1 = A->p [n] ;
cs_dropzeros (A) ; /* drop zero entries */
nz2 = A->p [n] ;
if (tol > 0) cs_droptol (A, tol) ; /* drop tiny entries (just to test) */
Prob->C = C = sym ? make_sym (A) : A ; /* C = A + triu(A,1)', or C=A */
if (!C) return (free_problem (Prob)) ;
printf ("\n--- Matrix: %g-by-%g, nnz: %g (sym: %g: nnz %g), norm: %8.2e\n",
(double) m, (double) n, (double) (A->p [n]), (double) sym,
(double) (sym ? C->p [n] : 0), cs_norm (C)) ;
if (nz1 != nz2) printf ("zero entries dropped: %g\n", (double) (nz1 - nz2));
if (nz2 != A->p [n]) printf ("tiny entries dropped: %g\n",
(double) (nz2 - A->p [n])) ;
Prob->b = cs_malloc (mn, sizeof (double)) ;
Prob->x = cs_malloc (mn, sizeof (double)) ;
Prob->resid = cs_malloc (mn, sizeof (double)) ;
return ((!Prob->b || !Prob->x || !Prob->resid) ? free_problem (Prob) : Prob) ;
}
/* free a problem */
problem *free_problem (problem *Prob)
{
if (!Prob) return (NULL) ;
cs_spfree (Prob->A) ;
if (Prob->sym) cs_spfree (Prob->C) ;
cs_free (Prob->b) ;
cs_free (Prob->x) ;
cs_free (Prob->resid) ;
return (cs_free (Prob)) ;
}
/* solve a linear system using Cholesky, LU, and QR, with various orderings */
csi demo2 (problem *Prob)
{
cs *A, *C ;
double *b, *x, *resid, t, tol ;
csi k, m, n, ok, order, nb, ns, *r, *s, *rr, sprank ;
csd *D ;
if (!Prob) return (0) ;
A = Prob->A ; C = Prob->C ; b = Prob->b ; x = Prob->x ; resid = Prob->resid;
m = A->m ; n = A->n ;
tol = Prob->sym ? 0.001 : 1 ; /* partial pivoting tolerance */
D = cs_dmperm (C, 1) ; /* randomized dmperm analysis */
if (!D) return (0) ;
nb = D->nb ; r = D->r ; s = D->s ; rr = D->rr ;
sprank = rr [3] ;
for (ns = 0, k = 0 ; k < nb ; k++)
{
ns += ((r [k+1] == r [k]+1) && (s [k+1] == s [k]+1)) ;
}
printf ("blocks: %g singletons: %g structural rank: %g\n",
(double) nb, (double) ns, (double) sprank) ;
cs_dfree (D) ;
for (order = 0 ; order <= 3 ; order += 3) /* natural and amd(A'*A) */
{
if (!order && m > 1000) continue ;
printf ("QR ") ;
print_order (order) ;
rhs (x, b, m) ; /* compute right-hand side */
t = tic () ;
ok = cs_qrsol (order, C, x) ; /* min norm(Ax-b) with QR */
printf ("time: %8.2f ", toc (t)) ;
print_resid (ok, C, x, b, resid) ; /* print residual */
}
if (m != n || sprank < n) return (1) ; /* return if rect. or singular*/
for (order = 0 ; order <= 3 ; order++) /* try all orderings */
{
if (!order && m > 1000) continue ;
printf ("LU ") ;
print_order (order) ;
rhs (x, b, m) ; /* compute right-hand side */
t = tic () ;
ok = cs_lusol (order, C, x, tol) ; /* solve Ax=b with LU */
printf ("time: %8.2f ", toc (t)) ;
print_resid (ok, C, x, b, resid) ; /* print residual */
}
if (!Prob->sym) return (1) ;
for (order = 0 ; order <= 1 ; order++) /* natural and amd(A+A') */
{
if (!order && m > 1000) continue ;
printf ("Chol ") ;
print_order (order) ;
rhs (x, b, m) ; /* compute right-hand side */
t = tic () ;
ok = cs_cholsol (order, C, x) ; /* solve Ax=b with Cholesky */
printf ("time: %8.2f ", toc (t)) ;
print_resid (ok, C, x, b, resid) ; /* print residual */
}
return (1) ;
}
/* free workspace for demo3 */
static csi done3 (csi ok, css *S, csn *N, double *y, cs *W, cs *E, csi *p)
{
cs_sfree (S) ;
cs_nfree (N) ;
cs_free (y) ;
cs_spfree (W) ;
cs_spfree (E) ;
cs_free (p) ;
return (ok) ;
}
/* Cholesky update/downdate */
csi demo3 (problem *Prob)
{
cs *A, *C, *W = NULL, *WW, *WT, *E = NULL, *W2 ;
csi n, k, *Li, *Lp, *Wi, *Wp, p1, p2, *p = NULL, ok ;
double *b, *x, *resid, *y = NULL, *Lx, *Wx, s, t, t1 ;
css *S = NULL ;
csn *N = NULL ;
if (!Prob || !Prob->sym || Prob->A->n == 0) return (0) ;
A = Prob->A ; C = Prob->C ; b = Prob->b ; x = Prob->x ; resid = Prob->resid;
n = A->n ;
if (!Prob->sym || n == 0) return (1) ;
rhs (x, b, n) ; /* compute right-hand side */
printf ("\nchol then update/downdate ") ;
print_order (1) ;
y = cs_malloc (n, sizeof (double)) ;
t = tic () ;
S = cs_schol (1, C) ; /* symbolic Chol, amd(A+A') */
printf ("\nsymbolic chol time %8.2f\n", toc (t)) ;
t = tic () ;
N = cs_chol (C, S) ; /* numeric Cholesky */
printf ("numeric chol time %8.2f\n", toc (t)) ;
if (!S || !N || !y) return (done3 (0, S, N, y, W, E, p)) ;
t = tic () ;
cs_ipvec (S->pinv, b, y, n) ; /* y = P*b */
cs_lsolve (N->L, y) ; /* y = L\y */
cs_ltsolve (N->L, y) ; /* y = L'\y */
cs_pvec (S->pinv, y, x, n) ; /* x = P'*y */
printf ("solve chol time %8.2f\n", toc (t)) ;
printf ("original: ") ;
print_resid (1, C, x, b, resid) ; /* print residual */
k = n/2 ; /* construct W */
W = cs_spalloc (n, 1, n, 1, 0) ;
if (!W) return (done3 (0, S, N, y, W, E, p)) ;
Lp = N->L->p ; Li = N->L->i ; Lx = N->L->x ;
Wp = W->p ; Wi = W->i ; Wx = W->x ;
Wp [0] = 0 ;
p1 = Lp [k] ;
Wp [1] = Lp [k+1] - p1 ;
s = Lx [p1] ;
srand (1) ;
for ( ; p1 < Lp [k+1] ; p1++)
{
p2 = p1 - Lp [k] ;
Wi [p2] = Li [p1] ;
Wx [p2] = s * rand () / ((double) RAND_MAX) ;
}
t = tic () ;
ok = cs_updown (N->L, +1, W, S->parent) ; /* update: L*L'+W*W' */
t1 = toc (t) ;
printf ("update: time: %8.2f\n", t1) ;
if (!ok) return (done3 (0, S, N, y, W, E, p)) ;
t = tic () ;
cs_ipvec (S->pinv, b, y, n) ; /* y = P*b */
cs_lsolve (N->L, y) ; /* y = L\y */
cs_ltsolve (N->L, y) ; /* y = L'\y */
cs_pvec (S->pinv, y, x, n) ; /* x = P'*y */
t = toc (t) ;
p = cs_pinv (S->pinv, n) ;
W2 = cs_permute (W, p, NULL, 1) ; /* E = C + (P'W)*(P'W)' */
WT = cs_transpose (W2,1) ;
WW = cs_multiply (W2, WT) ;
cs_spfree (WT) ;
cs_spfree (W2) ;
E = cs_add (C, WW, 1, 1) ;
cs_spfree (WW) ;
if (!E || !p) return (done3 (0, S, N, y, W, E, p)) ;
printf ("update: time: %8.2f (incl solve) ", t1+t) ;
print_resid (1, E, x, b, resid) ; /* print residual */
cs_nfree (N) ; /* clear N */
t = tic () ;
N = cs_chol (E, S) ; /* numeric Cholesky */
if (!N) return (done3 (0, S, N, y, W, E, p)) ;
cs_ipvec (S->pinv, b, y, n) ; /* y = P*b */
cs_lsolve (N->L, y) ; /* y = L\y */
cs_ltsolve (N->L, y) ; /* y = L'\y */
cs_pvec (S->pinv, y, x, n) ; /* x = P'*y */
t = toc (t) ;
printf ("rechol: time: %8.2f (incl solve) ", t) ;
print_resid (1, E, x, b, resid) ; /* print residual */
t = tic () ;
ok = cs_updown (N->L, -1, W, S->parent) ; /* downdate: L*L'-W*W' */
t1 = toc (t) ;
if (!ok) return (done3 (0, S, N, y, W, E, p)) ;
printf ("downdate: time: %8.2f\n", t1) ;
t = tic () ;
cs_ipvec (S->pinv, b, y, n) ; /* y = P*b */
cs_lsolve (N->L, y) ; /* y = L\y */
cs_ltsolve (N->L, y) ; /* y = L'\y */
cs_pvec (S->pinv, y, x, n) ; /* x = P'*y */
t = toc (t) ;
printf ("downdate: time: %8.2f (incl solve) ", t1+t) ;
print_resid (1, C, x, b, resid) ; /* print residual */
return (done3 (1, S, N, y, W, E, p)) ;
}

15
deps/csparse/Demo/cs_demo.h vendored Normal file
View file

@ -0,0 +1,15 @@
#include "cs.h"
typedef struct problem_struct
{
cs *A ;
cs *C ;
csi sym ;
double *x ;
double *b ;
double *resid ;
} problem ;
problem *get_problem (FILE *f, double tol) ;
csi demo2 (problem *Prob) ;
csi demo3 (problem *Prob) ;
problem *free_problem (problem *Prob) ;

177
deps/csparse/Demo/cs_demo.out vendored Normal file
View file

@ -0,0 +1,177 @@
( cd ../Lib ; make )
make[1]: Entering directory `/cise/research/sparse/SuiteSparse/CSparse/Lib'
make[1]: Nothing to be done for `all'.
make[1]: Leaving directory `/cise/research/sparse/SuiteSparse/CSparse/Lib'
cc -O -I../Include -o cs_demo1 cs_demo1.c ../Lib/libcsparse.a -lm
cc -O -I../Include -o cs_demo2 cs_demo2.c cs_demo.c ../Lib/libcsparse.a -lm
cc -O -I../Include -o cs_demo3 cs_demo3.c cs_demo.c ../Lib/libcsparse.a -lm
./cs_demo1 < ../Matrix/t1
T:
CSparse Version 3.1.3, Mar 26, 2014. Copyright (c) Timothy A. Davis, 2006-2014
triplet: 4-by-4, nzmax: 16 nnz: 10
2 2 : 3
1 0 : 3.1
3 3 : 1
0 2 : 3.2
1 1 : 2.9
3 0 : 3.5
3 1 : 0.4
1 3 : 0.9
0 0 : 4.5
2 1 : 1.7
A:
CSparse Version 3.1.3, Mar 26, 2014. Copyright (c) Timothy A. Davis, 2006-2014
4-by-4, nzmax: 10 nnz: 10, 1-norm: 11.1
col 0 : locations 0 to 2
1 : 3.1
3 : 3.5
0 : 4.5
col 1 : locations 3 to 5
1 : 2.9
3 : 0.4
2 : 1.7
col 2 : locations 6 to 7
2 : 3
0 : 3.2
col 3 : locations 8 to 9
3 : 1
1 : 0.9
AT:
CSparse Version 3.1.3, Mar 26, 2014. Copyright (c) Timothy A. Davis, 2006-2014
4-by-4, nzmax: 10 nnz: 10, 1-norm: 7.7
col 0 : locations 0 to 1
0 : 4.5
2 : 3.2
col 1 : locations 2 to 4
0 : 3.1
1 : 2.9
3 : 0.9
col 2 : locations 5 to 6
1 : 1.7
2 : 3
col 3 : locations 7 to 9
0 : 3.5
1 : 0.4
3 : 1
D:
CSparse Version 3.1.3, Mar 26, 2014. Copyright (c) Timothy A. Davis, 2006-2014
4-by-4, nzmax: 16 nnz: 16, 1-norm: 139.58
col 0 : locations 0 to 3
1 : 13.95
3 : 15.75
0 : 100.28
2 : 9.6
col 1 : locations 4 to 7
1 : 88.62
3 : 12.91
0 : 13.95
2 : 4.93
col 2 : locations 8 to 11
1 : 4.93
3 : 0.68
2 : 81.68
0 : 9.6
col 3 : locations 12 to 15
1 : 12.91
3 : 83.2
0 : 15.75
2 : 0.68
./cs_demo2 < ../Matrix/t1
--- Matrix: 4-by-4, nnz: 10 (sym: 0: nnz 0), norm: 1.11e+01
blocks: 1 singletons: 0 structural rank: 4
QR natural time: 0.00 resid: 3.06e-17
QR amd(A'*A) time: 0.00 resid: 3.06e-17
LU natural time: 0.00 resid: 1.53e-17
LU amd(A+A') time: 0.00 resid: 1.53e-17
LU amd(S'*S) time: 0.00 resid: 0.00e+00
LU amd(A'*A) time: 0.00 resid: 1.53e-17
./cs_demo2 < ../Matrix/ash219
--- Matrix: 219-by-85, nnz: 438 (sym: 0: nnz 0), norm: 9.00e+00
blocks: 1 singletons: 0 structural rank: 85
QR natural time: 0.00 resid: 1.61e-02
QR amd(A'*A) time: 0.00 resid: 1.61e-02
./cs_demo2 < ../Matrix/bcsstk01
--- Matrix: 48-by-48, nnz: 224 (sym: -1: nnz 400), norm: 3.57e+09
blocks: 1 singletons: 0 structural rank: 48
QR natural time: 0.00 resid: 2.62e-19
QR amd(A'*A) time: 0.00 resid: 5.27e-19
LU natural time: 0.00 resid: 2.17e-19
LU amd(A+A') time: 0.00 resid: 1.87e-19
LU amd(S'*S) time: 0.00 resid: 2.38e-19
LU amd(A'*A) time: 0.00 resid: 2.38e-19
Chol natural time: 0.00 resid: 2.64e-19
Chol amd(A+A') time: 0.00 resid: 2.55e-19
./cs_demo2 < ../Matrix/fs_183_1
--- Matrix: 183-by-183, nnz: 988 (sym: 0: nnz 0), norm: 1.70e+09
zero entries dropped: 71
tiny entries dropped: 10
blocks: 38 singletons: 37 structural rank: 183
QR natural time: 0.00 resid: 6.84e-28
QR amd(A'*A) time: 0.00 resid: 9.38e-28
LU natural time: 0.01 resid: 6.20e-28
LU amd(A+A') time: 0.00 resid: 1.55e-27
LU amd(S'*S) time: 0.00 resid: 6.98e-28
LU amd(A'*A) time: 0.00 resid: 6.98e-28
./cs_demo2 < ../Matrix/mbeacxc
--- Matrix: 492-by-490, nnz: 49920 (sym: 0: nnz 0), norm: 9.29e-01
blocks: 10 singletons: 8 structural rank: 448
QR natural time: 0.05 resid: nan
QR amd(A'*A) time: 0.06 resid: nan
./cs_demo2 < ../Matrix/west0067
--- Matrix: 67-by-67, nnz: 294 (sym: 0: nnz 0), norm: 6.14e+00
blocks: 2 singletons: 1 structural rank: 67
QR natural time: 0.00 resid: 7.14e-17
QR amd(A'*A) time: 0.00 resid: 3.10e-17
LU natural time: 0.00 resid: 3.89e-17
LU amd(A+A') time: 0.00 resid: 2.27e-17
LU amd(S'*S) time: 0.00 resid: 1.95e-17
LU amd(A'*A) time: 0.00 resid: 2.60e-17
./cs_demo2 < ../Matrix/lp_afiro
--- Matrix: 27-by-51, nnz: 102 (sym: 0: nnz 0), norm: 3.43e+00
blocks: 1 singletons: 0 structural rank: 27
QR natural time: 0.00 resid: 3.96e-16
QR amd(A'*A) time: 0.00 resid: 2.25e-16
./cs_demo2 < ../Matrix/bcsstk16
--- Matrix: 4884-by-4884, nnz: 147631 (sym: -1: nnz 290378), norm: 7.01e+09
blocks: 75 singletons: 74 structural rank: 4884
QR amd(A'*A) time: 0.61 resid: 1.39e-22
LU amd(A+A') time: 0.30 resid: 1.10e-22
LU amd(S'*S) time: 0.31 resid: 1.28e-22
LU amd(A'*A) time: 0.33 resid: 1.78e-22
Chol amd(A+A') time: 0.11 resid: 1.19e-22
./cs_demo3 < ../Matrix/bcsstk01
--- Matrix: 48-by-48, nnz: 224 (sym: -1: nnz 400), norm: 3.57e+09
chol then update/downdate amd(A+A')
symbolic chol time 0.00
numeric chol time 0.00
solve chol time 0.00
original: resid: 2.55e-19
update: time: 0.00
update: time: 0.00 (incl solve) resid: 9.66e-19
rechol: time: 0.00 (incl solve) resid: 1.55e-18
downdate: time: 0.00
downdate: time: 0.00 (incl solve) resid: 3.74e-17
./cs_demo3 < ../Matrix/bcsstk16
--- Matrix: 4884-by-4884, nnz: 147631 (sym: -1: nnz 290378), norm: 7.01e+09
chol then update/downdate amd(A+A')
symbolic chol time 0.00
numeric chol time 0.10
solve chol time 0.01
original: resid: 1.19e-22
update: time: 0.00
update: time: 0.00 (incl solve) resid: 1.12e-23
rechol: time: 0.11 (incl solve) resid: 1.17e-23
downdate: time: 0.00
downdate: time: 0.00 (incl solve) resid: 4.09e-22

27
deps/csparse/Demo/cs_demo1.c vendored Normal file
View file

@ -0,0 +1,27 @@
#include "cs.h"
int main (void)
{
cs *T, *A, *Eye, *AT, *C, *D ;
csi i, m ;
T = cs_load (stdin) ; /* load triplet matrix T from stdin */
printf ("T:\n") ; cs_print (T, 0) ; /* print T */
A = cs_compress (T) ; /* A = compressed-column form of T */
printf ("A:\n") ; cs_print (A, 0) ; /* print A */
cs_spfree (T) ; /* clear T */
AT = cs_transpose (A, 1) ; /* AT = A' */
printf ("AT:\n") ; cs_print (AT, 0) ; /* print AT */
m = A ? A->m : 0 ; /* m = # of rows of A */
T = cs_spalloc (m, m, m, 1, 1) ; /* create triplet identity matrix */
for (i = 0 ; i < m ; i++) cs_entry (T, i, i, 1) ;
Eye = cs_compress (T) ; /* Eye = speye (m) */
cs_spfree (T) ;
C = cs_multiply (A, AT) ; /* C = A*A' */
D = cs_add (C, Eye, 1, cs_norm (C)) ; /* D = C + Eye*norm (C,1) */
printf ("D:\n") ; cs_print (D, 0) ; /* print D */
cs_spfree (A) ; /* clear A AT C D Eye */
cs_spfree (AT) ;
cs_spfree (C) ;
cs_spfree (D) ;
cs_spfree (Eye) ;
return (0) ;
}

9
deps/csparse/Demo/cs_demo2.c vendored Normal file
View file

@ -0,0 +1,9 @@
#include "cs_demo.h"
/* cs_demo2: read a matrix and solve a linear system */
int main (void)
{
problem *Prob = get_problem (stdin, 1e-14) ;
demo2 (Prob) ;
free_problem (Prob) ;
return (0) ;
}

9
deps/csparse/Demo/cs_demo3.c vendored Normal file
View file

@ -0,0 +1,9 @@
#include "cs_demo.h"
/* cs_demo3: read a matrix and test Cholesky update/downdate */
int main (void)
{
problem *Prob = get_problem (stdin, 0) ;
demo3 (Prob) ;
free_problem (Prob) ;
return (0) ;
}

311
deps/csparse/Doc/ChangeLog vendored Normal file
View file

@ -0,0 +1,311 @@
Mar 26, 2014: version 3.1.3
* minor update to UFget
Apr 16, 2013: verison 3.1.2
* removed a useless line from cs_sqr.c (no integer overflow can
occur in S->lnz, S->unz at that point)
Jun 20, 2012: verison 3.1.1
* minor update for Windows (removed filesep)
Jun 1, 2012: version 3.1.0
* minor changes to enable the creation of CXSparse 3.1.0
Dec 7, 2011: version 3.0.2
* fixed the Makefile to better align with CFLAGS and other standards
* minor fix to MATLAB cs_install.m
Jan 18, 2010: version 3.0.0
* changed "int" to "csi", which is normally "ptrdiff_t". This way,
CSparse can be used on a 64-bit platform (and in 64-bit MATLAB
in particular), without the need to use CXSparse. If you wish
to use "int", simply edit CSparse/Include/cs.h and change the
definition of csi, or compile with -Dcsi=int.
Jan 25, 2011: version 2.2.5
* minor fix to "make install"
* minor change to cs_util.c, typecast return (cs_free (...))
* minor fixes to UFget: help file, UFsettings.txt
Nov 30, 2009: version 2.2.4
* added "make install" and "make uninstall"
* minor change to cs_make.m
Jan 20, 2009, v2.2.3
* all tabs expanded to 8 spaces (except for Makefiles)
* example corrected in cs_gaxpy.m
* minor change to cs_sparse to allow for i,j,x to be either
row or column vectors (or any mixture)
Sept 8, 2008, v2.2.2
* minor change to cs_make.m, to change "/" to filesep ("/" on Unix, "\"
on Windows), to handle limitations in some Windows compilers
Nov 1, 2007, v2.2.1
* very minor change to Include/cs.h: Changed name of 2nd argument of
cs_permute, to match the code. This has no affect on the code itself,
since the type ("int *") is unchanged. It's just a documentation issue.
* minor lint cleanup in mexFunctions
Mar 31, 2007, v2.2.0
* few changes to primary Source/ files. Changes mostly affect MATLAB
interface.
* Source/cs_house.c: correction to comment
* Souce/cs_updown.c: whitespace changed to reflect change in CXSparse,
no impact at all on CSparse itself.
* Doc/, Lib/ and Include/ directories created.
* Source/cs.h moved to Include/cs.h, version number changed to 2.2.0.
* modification to Makefiles, cs_make.m
* correction to help comments in cs_dmperm.m, cs_qr.m, cs_scc2.m, cs_scc.m
* if complex matrix passed to CSparse in MATLAB, error message
suggests using CXSparse instead
* minor performance fix to cs_sparse_mex.c
* cs_randperm added to MATLAB/Makefile; already appeared in cs_make.m
* minor improvements to MATLAB/ demos.
Mar 1, 2007, v2.1.0
* Source/cs_add.c: added test for matrix dimensions
* Source/cs_multiply.c: added test for matrix dimensions
* correction to MATLAB demo3 (no error in C code version of the demo)
* minor corrections to MATLAB m-file help comments.
Dec 12, 2006, v2.0.7
* minor MATLAB cleanup
Dec 6, 2006, v2.0.6
* Update to UFget. Now relies on the MATLAB urlwrite function instead of
my own Java code.
Nov 2006, v2.0.5
* Added UFgrep to UFget toolbox.
* No changes to C Source code, except for version number and date.
* Added two test matrices: ibm32a and ibm32b. ibm32a is the Harwell/
Boeing matrix ibm32, but with the last column removed. ibm32b
is the transpose of ibm32a. With optimization enabled (-O),
2 lines in cs_maxtrans.c are not tested; these matrices correct
that problem.
* Fixed UFget. Earlier version could not download matrices larger than
32MB.
* Modified UFget/UFweb, to reflect changes in the UF Sparse Matrix
Collection.
* Added ccspy.m and cs_scc2.m MATLAB functions
* Added examples to help info in each *.m MATLAB file
* modified cs_dmspy to speed up the plotting of large matrices with many
singletons
* minor change to cspy: now draws a box around the matrix.
* minor changes to MATLAB demos and tests.
Oct 13, 2006, v2.0.4
* minor modification to cs_updown.c. "n" was incorrectly declared "double".
It should be "int". This was safe, just a little confusing (n was only
used in an argument to cs_malloc, and is thus typecast).
Sept 28, 2006, v2.0.3
* minor modifications to MATLAB interface, to allow CSparse to be used
in MATLAB 6.5.
* added examples to m-files, other minor m-file cleanup.
* bug fix to cspy, to handle NaN's properly.
Aug 23, 2006: v2.0.2
* change to cs_updown mexFunction, to handle multiple rank updates
* symbolic links removed from Tcov/ directory in the distribution.
They are now created by Tcov/Makefile as needed. This makes the
zip files cleaner. Tcov/*test.c test files renamed.
July 19, 2006:
* minor fix to cs_must_compile.m and cs_make.m, to allow CSparse to be
compiled in MATLAB 6.5 with cs_make.
* minor fix to cspy for complex matrices (imaginary part was ignored).
* no change to version number or date, since there are no changes that
affect the appearance of CSparse in the book ("Direct Methods for
Sparse Linear Systems", SIAM, 2006).
June 24, 2006:
* minor typos in comments corrected. No change in code itself,
and no change in version number or date.
May 27, 2006, v2.0.1: (this version is printed in the book)
* minor bug fix. cs_util.c modified, so that cs_sprealloc (T,0) works
properly for a triplet matrix T (setting T->nzmax equal to T->nz).
The line in v2.0.0:
nzmax = (nzmax <= 0) ? (A->p [A->n]) : nzmax ;
changes to the following in v2.0.1:
if (nzmax <= 0) nzmax = (CS_CSC (A)) ? (A->p [A->n]) : A->nz ;
* minor typographical changes arising from the editting of the book.
Apr 12, 2006, v2.0.0:
* random permutation option added to cs_maxtrans and cs_dmperm, to help
avoid rare cases where the O(|A|n) time complexity is reached in
practice (GHS_indef/boyd2 in the UF sparse matrix collection, for
example). New cs_randperm function added.
Apr 10, 2006:
* stylistic changes for the book (except for the bug fix):
* "int order" parameter of cs_amd, cs_lusol, cs_cholsol, cs_qrsol, cs_sqr,
cs_schol changed. Now 0 means no ordering, 1 is A+A', 2 is S*S', and
3 is A*A'. In v1.2 and earlier, "order" took on a value ranging from
-1 to 2. "int n" parameter rearranged for cs_ipvec, cs_pvec, cs_post
(int n moved to the end). cs_triplet renamed cs_compress.
To ensure that these changes are propagated into user code that calls
CSparse, the "order" parameter has been placed as the first parameter
in all these routines. Your compiler will complain (gcc will halt) if
you upgrade from v1.2 to v2.0 without changing your code. This is much
better than a silent error in which you get the wrong ordering by
mistake (with a huge penalty in run-time performance and no compiler
warnings).
New syntax (v2.0 and later):
----------------------------
order = 0: natural ordering
order = 1: amd (A+A')
order = 2: amd (S'*S), where S=A except dense rows dropped
order = 3: amd (A'*A)
int cs_cholsol (int order, const cs *A, double *b) ;
int cs_lusol (int order, const cs *A, double *b, double tol) ;
int cs_qrsol (int order, const cs *A, double *b) ;
int *cs_amd (int order, const cs *A) ;
css *cs_schol (int order, const cs *A) ;
css *cs_sqr (int order, const cs *A, int qr) ;
int *cs_post (const int *parent, int n) ;
int cs_ipvec (const int *p, const double *b, double *x, int n) ;
int cs_pvec (const int *p, const double *b, double *x, int n) ;
cs *cs_compress (const cs *T) ;
Old syntax (v1.2 and earlier):
------------------------------
order = -1: natural ordering
order = 0: amd (A+A')
order = 1: amd (S'*S), where S=A except dense rows dropped
order = 2: amd (A'*A)
int cs_cholsol (const cs *A, double *b, int order) ;
int cs_lusol (const cs *A, double *b, int order, double tol) ;
int cs_qrsol (const cs *A, double *b, int order) ;
int *cs_amd (const cs *A, int order) ;
css *cs_schol (const cs *A, int order) ;
css *cs_sqr (const cs *A, int order, int qr) ;
int *cs_post (int n, const int *parent) ;
int cs_ipvec (int n, const int *p, const double *b, double *x) ;
int cs_pvec (int n, const int *p, const double *b, double *x) ;
cs *cs_triplet (const cs *T) ;
* CS_OVERFLOW macro removed (removed from cs_malloc.c; not needed).
* S->leftmost added to css (it was tacked onto S->pinv before).
* P,Q,R,S components of csd struct changed to p,q,r,s.
* Pinv and Q components of css struct changed to pinv and q.
* CS_CSC and CS_TRIPLET macros added, to clarify which CSparse functions
accept cs matrices in compressed column form, triplet form, or both.
* check for negative row/column indices added to cs_entry.
* cs_ereach and cs_leaf functions added.
* call to cs_sprealloc added to cs_fkeep.
* bug fixes in cs_counts and cs_amd (memory leak under rare out-of-memory
conditions).
Mar 15, 2006:
* cs_scc modified so that the row and columns of each component are put in
their natural order. cs_dmperm modified so that columns of each block
(instead of rows) are placed in their natural order.
* cs_splsolve renamed cs_spsolve, generalized to handle both Lx=b and
Ux=b, and non-unit diagonal, and placed in its own file (cs_spsolve.c;
it was a static function in cs_lu.c). cs_lsolve_mex.c and
cs_splsolve_mex.c merged (the latter was removed). cs_spsolve.c file
added
* cs_dmspy changed, so that block borders line up better with the matrix.
Mar 6, 2006:
* Makefile modified so that the Tcov tests (which may not be portable)
are not compiled and run by the default "make" in the CSparse directory.
To compile everything, including the Tcov tests, use "make all".
Trivial change to cs.h.
Feb 27, 2006:
* cs_reach, cs_dfs, cs_splsolve, cs_lu, and cs_scc changed to remove O(n)
initialized workspace.
* cs_reach and cs_splsolve now user-callable (were static in cs_lu.c).
Feb 20, 2006:
* various changes to simplify the construction of CXSparse
Feb 7, 2006:
* changed prototypes, adding "const" where appropriate.

19
deps/csparse/Doc/License.txt vendored Normal file
View file

@ -0,0 +1,19 @@
CSparse: a Concise Sparse matrix package.
Copyright (c) 2006, Timothy A. Davis.
http://www.suitesparse.com
--------------------------------------------------------------------------------
CSparse is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
CSparse is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this Module; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

504
deps/csparse/Doc/lesser.txt vendored Normal file
View file

@ -0,0 +1,504 @@
GNU LESSER GENERAL PUBLIC LICENSE
Version 2.1, February 1999
Copyright (C) 1991, 1999 Free Software Foundation, Inc.
51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Everyone is permitted to copy and distribute verbatim copies
of this license document, but changing it is not allowed.
[This is the first released version of the Lesser GPL. It also counts
as the successor of the GNU Library Public License, version 2, hence
the version number 2.1.]
Preamble
The licenses for most software are designed to take away your
freedom to share and change it. By contrast, the GNU General Public
Licenses are intended to guarantee your freedom to share and change
free software--to make sure the software is free for all its users.
This license, the Lesser General Public License, applies to some
specially designated software packages--typically libraries--of the
Free Software Foundation and other authors who decide to use it. You
can use it too, but we suggest you first think carefully about whether
this license or the ordinary General Public License is the better
strategy to use in any particular case, based on the explanations below.
When we speak of free software, we are referring to freedom of use,
not price. Our General Public Licenses are designed to make sure that
you have the freedom to distribute copies of free software (and charge
for this service if you wish); that you receive source code or can get
it if you want it; that you can change the software and use pieces of
it in new free programs; and that you are informed that you can do
these things.
To protect your rights, we need to make restrictions that forbid
distributors to deny you these rights or to ask you to surrender these
rights. These restrictions translate to certain responsibilities for
you if you distribute copies of the library or if you modify it.
For example, if you distribute copies of the library, whether gratis
or for a fee, you must give the recipients all the rights that we gave
you. You must make sure that they, too, receive or can get the source
code. If you link other code with the library, you must provide
complete object files to the recipients, so that they can relink them
with the library after making changes to the library and recompiling
it. And you must show them these terms so they know their rights.
We protect your rights with a two-step method: (1) we copyright the
library, and (2) we offer you this license, which gives you legal
permission to copy, distribute and/or modify the library.
To protect each distributor, we want to make it very clear that
there is no warranty for the free library. Also, if the library is
modified by someone else and passed on, the recipients should know
that what they have is not the original version, so that the original
author's reputation will not be affected by problems that might be
introduced by others.
Finally, software patents pose a constant threat to the existence of
any free program. We wish to make sure that a company cannot
effectively restrict the users of a free program by obtaining a
restrictive license from a patent holder. Therefore, we insist that
any patent license obtained for a version of the library must be
consistent with the full freedom of use specified in this license.
Most GNU software, including some libraries, is covered by the
ordinary GNU General Public License. This license, the GNU Lesser
General Public License, applies to certain designated libraries, and
is quite different from the ordinary General Public License. We use
this license for certain libraries in order to permit linking those
libraries into non-free programs.
When a program is linked with a library, whether statically or using
a shared library, the combination of the two is legally speaking a
combined work, a derivative of the original library. The ordinary
General Public License therefore permits such linking only if the
entire combination fits its criteria of freedom. The Lesser General
Public License permits more lax criteria for linking other code with
the library.
We call this license the "Lesser" General Public License because it
does Less to protect the user's freedom than the ordinary General
Public License. It also provides other free software developers Less
of an advantage over competing non-free programs. These disadvantages
are the reason we use the ordinary General Public License for many
libraries. However, the Lesser license provides advantages in certain
special circumstances.
For example, on rare occasions, there may be a special need to
encourage the widest possible use of a certain library, so that it becomes
a de-facto standard. To achieve this, non-free programs must be
allowed to use the library. A more frequent case is that a free
library does the same job as widely used non-free libraries. In this
case, there is little to gain by limiting the free library to free
software only, so we use the Lesser General Public License.
In other cases, permission to use a particular library in non-free
programs enables a greater number of people to use a large body of
free software. For example, permission to use the GNU C Library in
non-free programs enables many more people to use the whole GNU
operating system, as well as its variant, the GNU/Linux operating
system.
Although the Lesser General Public License is Less protective of the
users' freedom, it does ensure that the user of a program that is
linked with the Library has the freedom and the wherewithal to run
that program using a modified version of the Library.
The precise terms and conditions for copying, distribution and
modification follow. Pay close attention to the difference between a
"work based on the library" and a "work that uses the library". The
former contains code derived from the library, whereas the latter must
be combined with the library in order to run.
GNU LESSER GENERAL PUBLIC LICENSE
TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
0. This License Agreement applies to any software library or other
program which contains a notice placed by the copyright holder or
other authorized party saying it may be distributed under the terms of
this Lesser General Public License (also called "this License").
Each licensee is addressed as "you".
A "library" means a collection of software functions and/or data
prepared so as to be conveniently linked with application programs
(which use some of those functions and data) to form executables.
The "Library", below, refers to any such software library or work
which has been distributed under these terms. A "work based on the
Library" means either the Library or any derivative work under
copyright law: that is to say, a work containing the Library or a
portion of it, either verbatim or with modifications and/or translated
straightforwardly into another language. (Hereinafter, translation is
included without limitation in the term "modification".)
"Source code" for a work means the preferred form of the work for
making modifications to it. For a library, complete source code means
all the source code for all modules it contains, plus any associated
interface definition files, plus the scripts used to control compilation
and installation of the library.
Activities other than copying, distribution and modification are not
covered by this License; they are outside its scope. The act of
running a program using the Library is not restricted, and output from
such a program is covered only if its contents constitute a work based
on the Library (independent of the use of the Library in a tool for
writing it). Whether that is true depends on what the Library does
and what the program that uses the Library does.
1. You may copy and distribute verbatim copies of the Library's
complete source code as you receive it, in any medium, provided that
you conspicuously and appropriately publish on each copy an
appropriate copyright notice and disclaimer of warranty; keep intact
all the notices that refer to this License and to the absence of any
warranty; and distribute a copy of this License along with the
Library.
You may charge a fee for the physical act of transferring a copy,
and you may at your option offer warranty protection in exchange for a
fee.
2. You may modify your copy or copies of the Library or any portion
of it, thus forming a work based on the Library, and copy and
distribute such modifications or work under the terms of Section 1
above, provided that you also meet all of these conditions:
a) The modified work must itself be a software library.
b) You must cause the files modified to carry prominent notices
stating that you changed the files and the date of any change.
c) You must cause the whole of the work to be licensed at no
charge to all third parties under the terms of this License.
d) If a facility in the modified Library refers to a function or a
table of data to be supplied by an application program that uses
the facility, other than as an argument passed when the facility
is invoked, then you must make a good faith effort to ensure that,
in the event an application does not supply such function or
table, the facility still operates, and performs whatever part of
its purpose remains meaningful.
(For example, a function in a library to compute square roots has
a purpose that is entirely well-defined independent of the
application. Therefore, Subsection 2d requires that any
application-supplied function or table used by this function must
be optional: if the application does not supply it, the square
root function must still compute square roots.)
These requirements apply to the modified work as a whole. If
identifiable sections of that work are not derived from the Library,
and can be reasonably considered independent and separate works in
themselves, then this License, and its terms, do not apply to those
sections when you distribute them as separate works. But when you
distribute the same sections as part of a whole which is a work based
on the Library, the distribution of the whole must be on the terms of
this License, whose permissions for other licensees extend to the
entire whole, and thus to each and every part regardless of who wrote
it.
Thus, it is not the intent of this section to claim rights or contest
your rights to work written entirely by you; rather, the intent is to
exercise the right to control the distribution of derivative or
collective works based on the Library.
In addition, mere aggregation of another work not based on the Library
with the Library (or with a work based on the Library) on a volume of
a storage or distribution medium does not bring the other work under
the scope of this License.
3. You may opt to apply the terms of the ordinary GNU General Public
License instead of this License to a given copy of the Library. To do
this, you must alter all the notices that refer to this License, so
that they refer to the ordinary GNU General Public License, version 2,
instead of to this License. (If a newer version than version 2 of the
ordinary GNU General Public License has appeared, then you can specify
that version instead if you wish.) Do not make any other change in
these notices.
Once this change is made in a given copy, it is irreversible for
that copy, so the ordinary GNU General Public License applies to all
subsequent copies and derivative works made from that copy.
This option is useful when you wish to copy part of the code of
the Library into a program that is not a library.
4. You may copy and distribute the Library (or a portion or
derivative of it, under Section 2) in object code or executable form
under the terms of Sections 1 and 2 above provided that you accompany
it with the complete corresponding machine-readable source code, which
must be distributed under the terms of Sections 1 and 2 above on a
medium customarily used for software interchange.
If distribution of object code is made by offering access to copy
from a designated place, then offering equivalent access to copy the
source code from the same place satisfies the requirement to
distribute the source code, even though third parties are not
compelled to copy the source along with the object code.
5. A program that contains no derivative of any portion of the
Library, but is designed to work with the Library by being compiled or
linked with it, is called a "work that uses the Library". Such a
work, in isolation, is not a derivative work of the Library, and
therefore falls outside the scope of this License.
However, linking a "work that uses the Library" with the Library
creates an executable that is a derivative of the Library (because it
contains portions of the Library), rather than a "work that uses the
library". The executable is therefore covered by this License.
Section 6 states terms for distribution of such executables.
When a "work that uses the Library" uses material from a header file
that is part of the Library, the object code for the work may be a
derivative work of the Library even though the source code is not.
Whether this is true is especially significant if the work can be
linked without the Library, or if the work is itself a library. The
threshold for this to be true is not precisely defined by law.
If such an object file uses only numerical parameters, data
structure layouts and accessors, and small macros and small inline
functions (ten lines or less in length), then the use of the object
file is unrestricted, regardless of whether it is legally a derivative
work. (Executables containing this object code plus portions of the
Library will still fall under Section 6.)
Otherwise, if the work is a derivative of the Library, you may
distribute the object code for the work under the terms of Section 6.
Any executables containing that work also fall under Section 6,
whether or not they are linked directly with the Library itself.
6. As an exception to the Sections above, you may also combine or
link a "work that uses the Library" with the Library to produce a
work containing portions of the Library, and distribute that work
under terms of your choice, provided that the terms permit
modification of the work for the customer's own use and reverse
engineering for debugging such modifications.
You must give prominent notice with each copy of the work that the
Library is used in it and that the Library and its use are covered by
this License. You must supply a copy of this License. If the work
during execution displays copyright notices, you must include the
copyright notice for the Library among them, as well as a reference
directing the user to the copy of this License. Also, you must do one
of these things:
a) Accompany the work with the complete corresponding
machine-readable source code for the Library including whatever
changes were used in the work (which must be distributed under
Sections 1 and 2 above); and, if the work is an executable linked
with the Library, with the complete machine-readable "work that
uses the Library", as object code and/or source code, so that the
user can modify the Library and then relink to produce a modified
executable containing the modified Library. (It is understood
that the user who changes the contents of definitions files in the
Library will not necessarily be able to recompile the application
to use the modified definitions.)
b) Use a suitable shared library mechanism for linking with the
Library. A suitable mechanism is one that (1) uses at run time a
copy of the library already present on the user's computer system,
rather than copying library functions into the executable, and (2)
will operate properly with a modified version of the library, if
the user installs one, as long as the modified version is
interface-compatible with the version that the work was made with.
c) Accompany the work with a written offer, valid for at
least three years, to give the same user the materials
specified in Subsection 6a, above, for a charge no more
than the cost of performing this distribution.
d) If distribution of the work is made by offering access to copy
from a designated place, offer equivalent access to copy the above
specified materials from the same place.
e) Verify that the user has already received a copy of these
materials or that you have already sent this user a copy.
For an executable, the required form of the "work that uses the
Library" must include any data and utility programs needed for
reproducing the executable from it. However, as a special exception,
the materials to be distributed need not include anything that is
normally distributed (in either source or binary form) with the major
components (compiler, kernel, and so on) of the operating system on
which the executable runs, unless that component itself accompanies
the executable.
It may happen that this requirement contradicts the license
restrictions of other proprietary libraries that do not normally
accompany the operating system. Such a contradiction means you cannot
use both them and the Library together in an executable that you
distribute.
7. You may place library facilities that are a work based on the
Library side-by-side in a single library together with other library
facilities not covered by this License, and distribute such a combined
library, provided that the separate distribution of the work based on
the Library and of the other library facilities is otherwise
permitted, and provided that you do these two things:
a) Accompany the combined library with a copy of the same work
based on the Library, uncombined with any other library
facilities. This must be distributed under the terms of the
Sections above.
b) Give prominent notice with the combined library of the fact
that part of it is a work based on the Library, and explaining
where to find the accompanying uncombined form of the same work.
8. You may not copy, modify, sublicense, link with, or distribute
the Library except as expressly provided under this License. Any
attempt otherwise to copy, modify, sublicense, link with, or
distribute the Library is void, and will automatically terminate your
rights under this License. However, parties who have received copies,
or rights, from you under this License will not have their licenses
terminated so long as such parties remain in full compliance.
9. You are not required to accept this License, since you have not
signed it. However, nothing else grants you permission to modify or
distribute the Library or its derivative works. These actions are
prohibited by law if you do not accept this License. Therefore, by
modifying or distributing the Library (or any work based on the
Library), you indicate your acceptance of this License to do so, and
all its terms and conditions for copying, distributing or modifying
the Library or works based on it.
10. Each time you redistribute the Library (or any work based on the
Library), the recipient automatically receives a license from the
original licensor to copy, distribute, link with or modify the Library
subject to these terms and conditions. You may not impose any further
restrictions on the recipients' exercise of the rights granted herein.
You are not responsible for enforcing compliance by third parties with
this License.
11. If, as a consequence of a court judgment or allegation of patent
infringement or for any other reason (not limited to patent issues),
conditions are imposed on you (whether by court order, agreement or
otherwise) that contradict the conditions of this License, they do not
excuse you from the conditions of this License. If you cannot
distribute so as to satisfy simultaneously your obligations under this
License and any other pertinent obligations, then as a consequence you
may not distribute the Library at all. For example, if a patent
license would not permit royalty-free redistribution of the Library by
all those who receive copies directly or indirectly through you, then
the only way you could satisfy both it and this License would be to
refrain entirely from distribution of the Library.
If any portion of this section is held invalid or unenforceable under any
particular circumstance, the balance of the section is intended to apply,
and the section as a whole is intended to apply in other circumstances.
It is not the purpose of this section to induce you to infringe any
patents or other property right claims or to contest validity of any
such claims; this section has the sole purpose of protecting the
integrity of the free software distribution system which is
implemented by public license practices. Many people have made
generous contributions to the wide range of software distributed
through that system in reliance on consistent application of that
system; it is up to the author/donor to decide if he or she is willing
to distribute software through any other system and a licensee cannot
impose that choice.
This section is intended to make thoroughly clear what is believed to
be a consequence of the rest of this License.
12. If the distribution and/or use of the Library is restricted in
certain countries either by patents or by copyrighted interfaces, the
original copyright holder who places the Library under this License may add
an explicit geographical distribution limitation excluding those countries,
so that distribution is permitted only in or among countries not thus
excluded. In such case, this License incorporates the limitation as if
written in the body of this License.
13. The Free Software Foundation may publish revised and/or new
versions of the Lesser General Public License from time to time.
Such new versions will be similar in spirit to the present version,
but may differ in detail to address new problems or concerns.
Each version is given a distinguishing version number. If the Library
specifies a version number of this License which applies to it and
"any later version", you have the option of following the terms and
conditions either of that version or of any later version published by
the Free Software Foundation. If the Library does not specify a
license version number, you may choose any version ever published by
the Free Software Foundation.
14. If you wish to incorporate parts of the Library into other free
programs whose distribution conditions are incompatible with these,
write to the author to ask for permission. For software which is
copyrighted by the Free Software Foundation, write to the Free
Software Foundation; we sometimes make exceptions for this. Our
decision will be guided by the two goals of preserving the free status
of all derivatives of our free software and of promoting the sharing
and reuse of software generally.
NO WARRANTY
15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO
WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW.
EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR
OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY
KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE
LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME
THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN
WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY
AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU
FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR
CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE
LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING
RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A
FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF
SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH
DAMAGES.
END OF TERMS AND CONDITIONS
How to Apply These Terms to Your New Libraries
If you develop a new library, and you want it to be of the greatest
possible use to the public, we recommend making it free software that
everyone can redistribute and change. You can do so by permitting
redistribution under these terms (or, alternatively, under the terms of the
ordinary General Public License).
To apply these terms, attach the following notices to the library. It is
safest to attach them to the start of each source file to most effectively
convey the exclusion of warranty; and each file should have at least the
"copyright" line and a pointer to where the full notice is found.
<one line to give the library's name and a brief idea of what it does.>
Copyright (C) <year> <name of author>
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Also add information on how to contact you by electronic and paper mail.
You should also get your employer (if you work as a programmer) or your
school, if any, to sign a "copyright disclaimer" for the library, if
necessary. Here is a sample; alter the names:
Yoyodyne, Inc., hereby disclaims all copyright interest in the
library `Frob' (a library for tweaking knobs) written by James Random Hacker.
<signature of Ty Coon>, 1 April 1990
Ty Coon, President of Vice
That's all there is to it!

152
deps/csparse/Include/cs.h vendored Normal file
View file

@ -0,0 +1,152 @@
#ifndef _CS_H
#define _CS_H
#include <stdlib.h>
#include <limits.h>
#include <math.h>
#include <stdio.h>
#include <stddef.h>
#ifdef MATLAB_MEX_FILE
#include "mex.h"
#endif
#define CS_VER 3 /* CSparse Version */
#define CS_SUBVER 1
#define CS_SUBSUB 3
#define CS_DATE "Mar 26, 2014" /* CSparse release date */
#define CS_COPYRIGHT "Copyright (c) Timothy A. Davis, 2006-2014"
#ifdef MATLAB_MEX_FILE
#undef csi
#define csi mwSignedIndex
#endif
#ifndef csi
#define csi ptrdiff_t
#endif
/* --- primary CSparse routines and data structures ------------------------- */
typedef struct cs_sparse /* matrix in compressed-column or triplet form */
{
csi nzmax ; /* maximum number of entries */
csi m ; /* number of rows */
csi n ; /* number of columns */
csi *p ; /* column pointers (size n+1) or col indices (size nzmax) */
csi *i ; /* row indices, size nzmax */
double *x ; /* numerical values, size nzmax */
csi nz ; /* # of entries in triplet matrix, -1 for compressed-col */
} cs ;
cs *cs_add (const cs *A, const cs *B, double alpha, double beta) ;
csi cs_cholsol (csi order, const cs *A, double *b) ;
cs *cs_compress (const cs *T) ;
csi cs_dupl (cs *A) ;
csi cs_entry (cs *T, csi i, csi j, double x) ;
csi cs_gaxpy (const cs *A, const double *x, double *y) ;
cs *cs_load (FILE *f) ;
csi cs_lusol (csi order, const cs *A, double *b, double tol) ;
cs *cs_multiply (const cs *A, const cs *B) ;
double cs_norm (const cs *A) ;
csi cs_print (const cs *A, csi brief) ;
csi cs_qrsol (csi order, const cs *A, double *b) ;
cs *cs_transpose (const cs *A, csi values) ;
/* utilities */
void *cs_calloc (csi n, size_t size) ;
void *cs_free (void *p) ;
void *cs_realloc (void *p, csi n, size_t size, csi *ok) ;
cs *cs_spalloc (csi m, csi n, csi nzmax, csi values, csi triplet) ;
cs *cs_spfree (cs *A) ;
csi cs_sprealloc (cs *A, csi nzmax) ;
void *cs_malloc (csi n, size_t size) ;
/* --- secondary CSparse routines and data structures ----------------------- */
typedef struct cs_symbolic /* symbolic Cholesky, LU, or QR analysis */
{
csi *pinv ; /* inverse row perm. for QR, fill red. perm for Chol */
csi *q ; /* fill-reducing column permutation for LU and QR */
csi *parent ; /* elimination tree for Cholesky and QR */
csi *cp ; /* column pointers for Cholesky, row counts for QR */
csi *leftmost ; /* leftmost[i] = min(find(A(i,:))), for QR */
csi 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 */
csi *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 */
{
csi *p ; /* size m, row permutation */
csi *q ; /* size n, column permutation */
csi *r ; /* size nb+1, block k is rows r[k] to r[k+1]-1 in A(p,q) */
csi *s ; /* size nb+1, block k is cols s[k] to s[k+1]-1 in A(p,q) */
csi nb ; /* # of blocks in fine dmperm decomposition */
csi rr [5] ; /* coarse row decomposition */
csi cc [5] ; /* coarse column decomposition */
} csd ;
csi *cs_amd (csi order, const cs *A) ;
csn *cs_chol (const cs *A, const css *S) ;
csd *cs_dmperm (const cs *A, csi seed) ;
csi cs_droptol (cs *A, double tol) ;
csi cs_dropzeros (cs *A) ;
csi cs_happly (const cs *V, csi i, double beta, double *x) ;
csi cs_ipvec (const csi *p, const double *b, double *x, csi n) ;
csi cs_lsolve (const cs *L, double *x) ;
csi 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 csi *pinv, const csi *q, csi values) ;
csi *cs_pinv (const csi *p, csi n) ;
csi cs_pvec (const csi *p, const double *b, double *x, csi n) ;
csn *cs_qr (const cs *A, const css *S) ;
css *cs_schol (csi order, const cs *A) ;
css *cs_sqr (csi order, const cs *A, csi qr) ;
cs *cs_symperm (const cs *A, const csi *pinv, csi values) ;
csi cs_updown (cs *L, csi sigma, const cs *C, const csi *parent) ;
csi cs_usolve (const cs *U, double *x) ;
csi 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 -------------------------------------------- */
csi *cs_counts (const cs *A, const csi *parent, const csi *post, csi ata) ;
double cs_cumsum (csi *p, csi *c, csi n) ;
csi cs_dfs (csi j, cs *G, csi top, csi *xi, csi *pstack, const csi *pinv) ;
csi cs_ereach (const cs *A, csi k, const csi *parent, csi *s, csi *w) ;
csi *cs_etree (const cs *A, csi ata) ;
csi cs_fkeep (cs *A, csi (*fkeep) (csi, csi, double, void *), void *other) ;
double cs_house (double *x, double *beta, csi n) ;
csi cs_leaf (csi i, csi j, const csi *first, csi *maxfirst, csi *prevleaf,
csi *ancestor, csi *jleaf) ;
csi *cs_maxtrans (const cs *A, csi seed) ;
csi *cs_post (const csi *parent, csi n) ;
csi *cs_randperm (csi n, csi seed) ;
csi cs_reach (cs *G, const cs *B, csi k, csi *xi, const csi *pinv) ;
csi cs_scatter (const cs *A, csi j, double beta, csi *w, double *x, csi mark,
cs *C, csi nz) ;
csd *cs_scc (cs *A) ;
csi cs_spsolve (cs *G, const cs *B, csi k, csi *xi, double *x,
const csi *pinv, csi lo) ;
csi cs_tdfs (csi j, csi k, csi *head, const csi *next, csi *post,
csi *stack) ;
/* utilities */
csd *cs_dalloc (csi m, csi n) ;
csd *cs_ddone (csd *D, cs *C, void *w, csi ok) ;
cs *cs_done (cs *C, void *w, void *x, csi ok) ;
csi *cs_idone (csi *p, cs *C, void *w, csi ok) ;
csn *cs_ndone (csn *N, cs *C, void *w, void *x, csi 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

33
deps/csparse/Lib/Makefile vendored Normal file
View file

@ -0,0 +1,33 @@
CF = $(CFLAGS) $(CPPFLAGS) $(TARGET_ARCH) -O
I = -I../Include
RANLIB = ranlib
ARCHIVE = $(AR) $(ARFLAGS)
all: libcsparse.a
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
$(CS): ../Include/cs.h Makefile
%.o: ../Source/%.c ../Include/cs.h
$(CC) $(CF) $(I) -c $<
libcsparse.a: $(CS)
$(ARCHIVE) libcsparse.a $(CS)
- $(RANLIB) libcsparse.a
clean:
- $(RM) *.o
purge: distclean
distclean: clean
- $(RM) *.a *.obj *.dll

40
deps/csparse/Makefile vendored Normal file
View file

@ -0,0 +1,40 @@
#------------------------------------------------------------------------------
# CSparse Makefile
#------------------------------------------------------------------------------
VERSION = 3.1.3
C:
( cd Lib ; $(MAKE) )
( cd Demo ; $(MAKE) )
all: C cov
library:
( cd Lib ; $(MAKE) )
cov:
( cd Tcov ; $(MAKE) )
clean:
( cd Lib ; $(MAKE) clean )
( cd Demo ; $(MAKE) clean )
( cd Tcov ; $(MAKE) clean )
( cd MATLAB/CSparse ; $(RM) *.o )
( cd MATLAB/Test ; $(RM) *.o )
purge:
( cd Lib ; $(MAKE) purge )
( cd Demo ; $(MAKE) purge )
( cd Tcov ; $(MAKE) purge )
( cd MATLAB/CSparse ; $(RM) *.o *.mex* )
( cd MATLAB/Test ; $(RM) *.o *.mex* )
distclean: purge
# do not install CSparse; use CXSparse instead
install:
# uninstall CSparse: do nothing
uninstall:

434
deps/csparse/README.txt vendored Normal file
View file

@ -0,0 +1,434 @@
###
This is a copy of the CSParse library without the /MATLAB and Matrix/
folders.
###
CSparse: a Concise Sparse Matrix package.
VERSION 3.1.3, Copyright (c) 2006-2014, Timothy A. Davis, Mar 26, 2014
http://www.suitesparse.com
Refer to "Direct Methods for Sparse Linear Systems," Timothy A. Davis,
SIAM, Philadelphia, 2006. No detailed user guide is included in this
package; the user guide is the book itself.
The algorithms contained in CSparse have been chosen with five goals in mind:
(1) they must embody much of the theory behind sparse matrix algorithms,
(2) they must be either asymptotically optimal in their run time and memory
usage or be fast in practice,
(3) they must be concise so as to be easily understood and short enough to
print in the book,
(4) they must cover a wide spectrum of matrix operations, and
(5) they must be accurate and robust.
The focus is on direct methods; iterative methods and solvers for
eigenvalue problems are beyond the scope of this package.
No detailed user guide is included in this package; the user guide is the
book itself. Some indication of how to call the CSparse C routines is given
by the M-files in the MATLAB/CSparse directory.
Complex matrices are not supported, except for methods that operate only
on the nonzero pattern of a matrix. A complex version of CSparse appears
as a separate package, CXSparse ("Concise Extended Sparse matrix package").
The performance of the sparse factorization methods in CSparse will not be
competitive with UMFPACK or CHOLMOD, but the codes are much more concise and
easy to understand (see the above goals). Other methods are competitive.
Some of the MATLAB tests require the AMD package.
See http://www.suitesparse.com for CSparse and the AMD ordering
package. See the Doc/License.txt file for the license (GNU LGPL).
To compile the C-library (./Source) and C demo programs (./Demo) just type
"make" in this directory. This should work on any system with the "make"
command. To run the exhaustive tests, type "make" in the Tcov directory
(Linux is assumed). To compile the MATLAB mexFunctions type "make mex" in
this directory, or just "make" in the MATLAB directory. To remove all files
not in the original distribution, type "make distclean". I recommend that you
use a different level of optimization than "cc -O", which was chosen so that
the Makefile is portable. See Source/Makefile.
You can simply type "cs_install" while in the CSparse/MATLAB directory to
compile and install CSparse for use in MATLAB. This is especially useful for
a typical Microsoft Windows system, which does not include "make". For more
details, see CSparse/MATLAB/README.txt.
--------------------------------------------------------------------------------
Contents:
--------------------------------------------------------------------------------
Demo/ demo C programs that use CSparse
Doc/ license and change log
Makefile Makefile for the whole package
MATLAB/ MATLAB interface, demos, and tests for CSparse
Matrix/ sample matrices
README.txt this file
Source/ primary CSparse source files (C only, no MATLAB)
Tcov/ CSparse tests
--------------------------------------------------------------------------------
./Doc: license and change log
--------------------------------------------------------------------------------
ChangeLog changes in CSparse since first release
lesser.txt the GNU LGPL
License.txt license (GNU LGPL)
--------------------------------------------------------------------------------
./Source: Primary source code for CSparse
--------------------------------------------------------------------------------
cs_add.c add sparse matrices
cs_amd.c approximate minimum degree
cs_chol.c sparse Cholesky
cs_cholsol.c x=A\b using sparse Cholesky
cs_compress.c convert a triplet form to compressed-column form
cs_counts.c column counts for Cholesky and QR
cs_cumsum.c cumulative sum
cs_dfs.c depth-first-search
cs_dmperm.c Dulmage-Mendelsohn permutation
cs_droptol.c drop small entries from a sparse matrix
cs_dropzeros.c drop zeros from a sparse matrix
cs_dupl.c remove (and sum) duplicates
cs_entry.c add an entry to a triplet matrix
cs_ereach.c nonzero pattern of Cholesky L(k,:) from etree and triu(A(:,k))
cs_etree.c find elimination tree
cs_fkeep.c drop entries from a sparse matrix
cs_gaxpy.c sparse matrix times dense matrix
cs.h include file for CSparse
cs_happly.c apply Householder reflection
cs_house.c compute Householder reflection
cs_ipvec.c x(p)=b
cs_leaf.c determine if j is a leaf of the skeleton matrix and find lca
cs_load.c load a sparse matrix from a file
cs_lsolve.c x=L\b
cs_ltsolve.c x=L'\b
cs_lu.c sparse LU factorization
cs_lusol.c x=A\b using sparse LU factorization
cs_malloc.c memory manager
cs_maxtrans.c maximum transveral (permutation for zero-free diagonal)
cs_multiply.c sparse matrix multiply
cs_norm.c sparse matrix norm
cs_permute.c permute a sparse matrix
cs_pinv.c invert a permutation vector
cs_post.c postorder an elimination tree
cs_print.c print a sparse matrix
cs_pvec.c x=b(p)
cs_qr.c sparse QR
cs_qrsol.c solve a least-squares problem
cs_randperm.c random permutation
cs_reach.c find nonzero pattern of x=L\b for sparse L and b
cs_scatter.c scatter a sparse vector
cs_scc.c strongly-connected components
cs_schol.c symbolic Cholesky
cs_spsolve.c x=G\b where G, x, and b are sparse, and G upper/lower triangular
cs_sqr.c symbolic QR (also can be used for LU)
cs_symperm.c symmetric permutation of a sparse matrix
cs_tdfs.c depth-first-search of a tree
cs_transpose.c transpose a sparse matrix
cs_updown.c sparse rank-1 Cholesky update/downate
cs_usolve.c x=U\b
cs_util.c various utilities (allocate/free matrices, workspace, etc)
cs_utsolve.c x=U'\b
Makefile Makefile for CSparse
README.txt README file for CSparse
--------------------------------------------------------------------------------
./Demo: C program demos
--------------------------------------------------------------------------------
cs_demo1.c read a matrix from a file and perform basic matrix operations
cs_demo2.c read a matrix from a file and solve a linear system
cs_demo3.c read a matrix, solve a linear system, update/downdate
cs_demo.c support routines for cs_demo*.c
cs_demo.h include file for demo programs
cs_demo.out output of "make", which runs the demos on some matrices
Makefile Makefile for Demo programs
readhb.f read a Rutherford-Boeing matrix
README.txt Demo README file
--------------------------------------------------------------------------------
./MATLAB: MATLAB interface, demos, and tests
--------------------------------------------------------------------------------
cs_install.m MATLAB function for compiling and installing CSparse for MATLAB
CSparse/ MATLAB interface for CSparse
Demo/ MATLAB demos for CSparse
Makefile MATLAB interface Makefile
README.txt MATLAB README file
Test/ MATLAB test for CSparse, and "textbook" routines
UFget/ MATLAB interface to UF Sparse Matrix Collection
--------------------------------------------------------------------------------
./MATLAB/CSparse: MATLAB interface for CSparse
--------------------------------------------------------------------------------
Contents.m Contents of MATLAB interface to CSparse
cs_add.m add two sparse matrices
cs_add_mex.c
cs_amd.m approximate minimum degree
cs_amd_mex.c
cs_chol.m sparse Cholesky
cs_chol_mex.c
cs_cholsol.m x=A\b using a sparse Cholesky
cs_cholsol_mex.c
cs_counts.m column counts for Cholesky or QR (like "symbfact" in MATLAB)
cs_counts_mex.c
cs_dmperm.m Dulmage-Mendelsohn permutation
cs_dmperm_mex.c
cs_dmsol.m x=A\b using dmperm
cs_dmspy.m plot a picture of a dmperm-permuted matrix
cs_droptol.m drop small entries
cs_droptol_mex.c
cs_esep.m find edge separator
cs_etree.m compute elimination tree
cs_etree_mex.c
cs_gaxpy.m sparse matrix times dense vector
cs_gaxpy_mex.c
cs_lsolve.m x=L\b where L is lower triangular
cs_lsolve_mex.c
cs_ltsolve.m x=L'\b where L is lower triangular
cs_ltsolve_mex.c
cs_lu.m sparse LU factorization
cs_lu_mex.c
cs_lusol.m x=A\b using sparse LU factorization
cs_lusol_mex.c
cs_make.m compiles CSparse for use in MATLAB
cs_mex.c support routines for CSparse mexFunctions
cs_mex.h
cs_multiply.m sparse matrix multiply
cs_multiply_mex.c
cs_must_compile.m determine if a source file needs to be compiled with mex
cs_nd.m nested dissection
cs_nsep.m find node separator
cs_permute.m permute a sparse matrix
cs_permute_mex.c
cs_print.m print a sparse matrix
cs_print_mex.c
cs_qleft.m apply Householder vectors to the left
cs_qright.m apply Householder vectors to the right
cs_qr.m sparse QR factorization
cs_qr_mex.c
cs_qrsol.m solve a sparse least squares problem
cs_qrsol_mex.c
cs_randperm.m randdom permutation
cs_randperm_mex.c
cs_scc.m strongly-connected components
cs_scc_mex.c
cs_sep.m convert an edge separator into a node separator
cs_sparse.m convert a triplet form matrix to a compress-column form
cs_sparse_mex.c
cs_symperm.m symmetric permutation of a sparse matrix
cs_symperm_mex.c
cs_sqr.m symbolic QR ordering and analysis
cs_sqr_mex.c
cs_thumb_mex.c compute small "thumbnail" of a sparse matrix (for cspy).
cs_transpose.m transpose a sparse matrix
cs_transpose_mex.c
cs_updown.m sparse Cholesky update/downdate
cs_updown_mex.c
cs_usolve.m x=U\b where U is upper triangular
cs_usolve_mex.c
cs_utsolve.m x=U'\b where U is upper triangular
cs_utsolve_mex.c
cspy.m a color "spy"
Makefile Makefile for CSparse MATLAB interface
README.txt README file for CSparse MATLAB interface
--------------------------------------------------------------------------------
./MATLAB/Demo: MATLAB demos for CSparse
--------------------------------------------------------------------------------
Contents.m Contents of MATLAB demo for CSparse
cs_demo.m run all MATLAB demos for CSparse
cs_demo1.m MATLAB version of Demo/cs_demo1.c
cs_demo2.m MATLAB version of Demo/cs_demo2.c
cs_demo3.m MATLAB version of Demo/cs_demo3.c
private/ private functions for MATLAB demos
README.txt README file for CSparse MATLAB demo
--------------------------------------------------------------------------------
./MATLAB/Demo/private: private functions for MATLAB demos
--------------------------------------------------------------------------------
demo2.m demo 2
demo3.m demo 3
ex_1.m example 1
ex2.m example 2
ex3.m example 3
frand.m generate a random finite-element matrix
get_problem.m get a matrix
is_sym.m determine if a matrix is symmetric
mesh2d1.m construct a 2D mesh (method 1)
mesh2d2.m construct a 2D mesh (method 2)
mesh3d1.m construct a 3D mesh (method 1)
mesh3d2.m construct a 3D mesh (method 2)
print_order.m print the ordering method used
resid.m compute residual
rhs.m create right-hand-side
--------------------------------------------------------------------------------
./MATLAB/Test: Extensive test of CSparse, in MATLAB
--------------------------------------------------------------------------------
Makefile Makefile for MATLAB Test directory
README.txt README file for MATLAB/Test
Contents.m Contents of MATLAB/Test, "textbook" files only
chol_downdate.m downdate a Cholesky factorization.
chol_left.m left-looking Cholesky factorization.
chol_left2.m left-looking Cholesky factorization, more details.
chol_right.m right-looking Cholesky factorization.
chol_super.m left-looking "supernodal" Cholesky factorization.
chol_up.m up-looking Cholesky factorization.
chol_update.m update a Cholesky factorization.
chol_updown.m update or downdate a Cholesky factorization.
cond1est.m 1-norm condition estimate.
cs_fiedler.m the Fiedler vector of a connected graph.
givens2.m find a Givens rotation.
house.m find a Householder reflection.
lu_left.m left-looking LU factorization.
lu_right.m right-looking LU factorization.
lu_rightp.m right-looking LU factorization, with partial pivoting.
lu_rightpr.m recursive right-looking LU, with partial pivoting.
lu_rightr.m recursive right-looking LU.
norm1est.m 1-norm estimate.
qr_givens.m Givens-rotation QR factorization.
qr_givens_full.m Givens-rotation QR factorization, for full matrices.
qr_left.m left-looking Householder QR factorization.
qr_right.m right-looking Householder QR factorization.
cs_fiedler.m Fiedler vector
cs_frand.m generate a random finite-element matrix
cs_frand_mex.c
cs_ipvec.m x(p)=b
cs_ipvec_mex.c
cs_maxtransr.m recursive maximum matching algorithm
cs_maxtransr_mex.c
cs_pvec.m x=b(p)
cs_pvec_mex.c interface for cs_pvec
cs_reach.m non-recursive reach (interface to CSparse cs_reach)
cs_reach_mex.c non-recursive x=spones(L\sparse(b))
cs_reachr.m recursive reach (interface to CSparse cs_reachr)
cs_reachr_mex.c
cs_rowcnt.m row counts for sparse Cholesky
cs_rowcnt_mex.c row counts for sparse Cholesky
cs_sparse2.m same as cs_sparse, to test cs_entry function
cs_sparse2_mex.c like cs_sparse, but for testing cs_entry
cs_test_make.m compiles MATLAB tests
check_if_same.m check if two inputs are identical or not
choldn.m Cholesky downdate
cholup.m Cholesky update, using Given's rotations
cholupdown.m Cholesky update/downdate (Bischof, Pan, and Tang method)
cs_q1.m construct Q from Householder vectors
cs_test_make.m compiles the CSparse, Demo, and Test mexFunctions.
dmperm_test.m test cs_dmperm
chol_example.m simple Cholesky factorization example
etree_sample.m construct a sample etree and symbolic factorization
gqr3.m QR factorization, based on Givens rotations
happly.m apply Householder reflection to a vector
hmake1.m construct a Householder reflection
mynormest1.m estimate norm(A,1), using LU factorization (L*U = P*A*Q).
myqr.m QR factorization using Householder reflections
another_colormap.m try another color map
cspy_test.m test cspy and cs_dmspy
qr2.m QR factorization based on Householder reflections
sample_colormap.m try a colormap for use in cspy
signum.m compute and display the sign of a column vector x
sqr_example.m test cs_sqr
dmspy_test.m test cspy, cs_dmspy, and cs_dmperm
test_qr.m test various QR factorization methods
test_randperms.m test random permutations
testh.m test Householder reflections
test_qr1.m test QR factorizations
test_qrsol.m test cs_qrsol
test_sep.m test cs_sep, and compare with Gilbert's meshpart vtxsep
testall.m test all CSparse functions (run tests 1 to 28 below)
test1.m test cs_transpose
test2.m test cs_sparse
test3.m test cs_lsolve, cs_ltsolve, cs_usolve, cs_chol
test4.m test cs_multiply
test5.m test cs_add
test6.m test cs_reach, cs_reachr, cs_lsolve, cs_usolve
test7.m test cs_lu
test8.m test cs_cholsol, cs_lusol
test9.m test cs_qr
test10.m test cs_qr
test11.m test cs_rowcnt
test12.m test cs_qr and compare with svd
test13.m test cs_counts, cs_etree
test14.m test cs_droptol
test15.m test cs_amd
test16.m test cs_amd
test17.m test cs_qr, cs_qright, cs_q1, cs_qrleft, cs_qrsol
test18.m test iterative refinement after backslash
test19.m test cs_dmperm, cs_maxtransr, cs_dmspy, cs_scc
test20.m test cholupdown
test21.m test cs_updown
test22.m test cond1est
test23.m test cs_dmspy
test24.m test cs_fielder
test25.m test cs_nd
test26.m test cs_dmsol and cs_dmspy
test27.m test cs_qr, cs_utsolve, cs_qrsol
test28.m test cs_randperm, cs_dmperm
--------------------------------------------------------------------------------
./MATLAB/UFget: MATLAB interface for the UF Sparse Matrix Collection
--------------------------------------------------------------------------------
Contents.m Contents of UFget
mat/ default directory where downloaded matrices will be put
README.txt README file for UFget
UFget_defaults.m default parameter settings
UFget_example.m example of use
UFget_install.m installs UFget temporarily (for current session)
UFget_java.class read a url and load it in into MATLAB (compiled Java code)
UFget_java.java read a url and load it in into MATLAB (Java source code)
UFget_lookup.m look up a matrix in the index
UFget.m UFget itself (primary user interface)
UFweb.m open url for a matrix or collection
mat/UF_Index.mat index of matrices in UF Sparse Matrix Collection
--------------------------------------------------------------------------------
./Matrix: Sample matrices, most from Rutherford/Boeing collection
--------------------------------------------------------------------------------
ash219 overdetermined pattern of Holland survey. Ashkenazi, 1974.
bcsstk01 stiffness matrix for small generalized eigenvalue problem
bcsstk16 stiffness matrix, Corp of Engineers dam
fs_183_1 unsymmetric facsimile convergence matrix
lp_afiro NETLIB afiro linear programming problem
mbeacxc US economy, 1972. Dan Szyld, while at NYU
t1 small example used in Chapter 2
west0067 Cavett problem with 5 components (chemical eng., Westerberg)
--------------------------------------------------------------------------------
./Tcov: Exhaustive test coverage of CSparse
--------------------------------------------------------------------------------
covall same as covall.linux
covall.linux find coverage (Linux)
covall.sol find coverage (Solaris)
cov.awk coverage summary
cover print uncovered lines
covs print uncovered lines
cstcov_malloc_test.c malloc test
cstcov_malloc_test.h
cstcov_test.c main program for Tcov tests
gcovs run gcov (Linux)
Makefile Makefile for Tcov tests
nil an empty matrix
zero a 1-by-1 zero matrix
README.txt README file for Tcov directory

5
deps/csparse/Source/README.txt vendored Normal file
View file

@ -0,0 +1,5 @@
CSparse/Source directory: primary ANSI C source code files for CSparse.
All of these files are printed verbatim in the book. To compile the
libcsparse.a C-callable library, just type "make" in this directory.
Timothy A. Davis, http://www.suitesparse.com

28
deps/csparse/Source/cs_add.c vendored Normal file
View file

@ -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)
{
csi 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 (csi)) ; /* 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 */
}

364
deps/csparse/Source/cs_amd.c vendored Normal file
View file

@ -0,0 +1,364 @@
#include "cs.h"
/* clear w */
static csi cs_wclear (csi mark, csi lemax, csi *w, csi n)
{
csi 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 csi cs_diag (csi i, csi j, double aij, void *other) { return (i != j) ; }
/* p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */
csi *cs_amd (csi order, const cs *A) /* order 0:natural, 1:Chol, 2:LU, 3:QR */
{
cs *C, *A2, *AT ;
csi *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 ;
csi 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 (csi)) ; /* allocate result */
W = cs_malloc (8*(n+1), sizeof (csi)) ; /* 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 = ((h<0) ? (-h):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)) ;
}

58
deps/csparse/Source/cs_chol.c vendored Normal file
View file

@ -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 ;
csi 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 (csi)) ; /* get csi 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 */
}

26
deps/csparse/Source/cs_cholsol.c vendored Normal file
View file

@ -0,0 +1,26 @@
#include "cs.h"
/* x=A\b where A is symmetric positive definite; b overwritten with solution */
csi cs_cholsol (csi order, const cs *A, double *b)
{
double *x ;
css *S ;
csn *N ;
csi 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) ;
}

22
deps/csparse/Source/cs_compress.c vendored Normal file
View file

@ -0,0 +1,22 @@
#include "cs.h"
/* C = compressed-column form of a triplet matrix T */
cs *cs_compress (const cs *T)
{
csi 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 (csi)) ; /* 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 */
}

61
deps/csparse/Source/cs_counts.c vendored Normal file
View file

@ -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 csi *post, csi *w, csi **head, csi **next)
{
csi 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 ;
}
}
csi *cs_counts (const cs *A, const csi *parent, const csi *post, csi ata)
{
csi 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 (csi)) ; /* allocate result */
w = cs_malloc (s, sizeof (csi)) ; /* 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 */
}

17
deps/csparse/Source/cs_cumsum.c vendored Normal file
View file

@ -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 (csi *p, csi *c, csi n)
{
csi 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 csi 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]) */
}

36
deps/csparse/Source/cs_dfs.c vendored Normal file
View file

@ -0,0 +1,36 @@
#include "cs.h"
/* depth-first-search of the graph of a matrix, starting at node j */
csi cs_dfs (csi j, cs *G, csi top, csi *xi, csi *pstack, const csi *pinv)
{
csi 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) ;
}

144
deps/csparse/Source/cs_dmperm.c vendored Normal file
View file

@ -0,0 +1,144 @@
#include "cs.h"
/* breadth-first search for coarse decomposition (C0,C1,R1 or R0,R3,C3) */
static csi cs_bfs (const cs *A, csi n, csi *wi, csi *wj, csi *queue,
const csi *imatch, const csi *jmatch, csi mark)
{
csi *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 (csi n, const csi *wj, const csi *imatch, csi *p, csi *q,
csi *cc, csi *rr, csi set, csi mark)
{
csi kc = cc [set], j ;
csi 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 (csi m, const csi *wi, csi *p, csi *rr, csi set)
{
csi 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 csi cs_rprune (csi i, csi j, double aij, void *other)
{
csi *rr = (csi *) other ;
return (i >= rr [1] && i < rr [2]) ;
}
/* Given A, compute coarse and then fine dmperm */
csd *cs_dmperm (const cs *A, csi seed)
{
csi 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)) ;
}

9
deps/csparse/Source/cs_droptol.c vendored Normal file
View file

@ -0,0 +1,9 @@
#include "cs.h"
static csi cs_tol (csi i, csi j, double aij, void *tol)
{
return (fabs (aij) > *((double *) tol)) ;
}
csi cs_droptol (cs *A, double tol)
{
return (cs_fkeep (A, &cs_tol, &tol)) ; /* keep all large entries */
}

9
deps/csparse/Source/cs_dropzeros.c vendored Normal file
View file

@ -0,0 +1,9 @@
#include "cs.h"
static csi cs_nonzero (csi i, csi j, double aij, void *other)
{
return (aij != 0) ;
}
csi cs_dropzeros (cs *A)
{
return (cs_fkeep (A, &cs_nonzero, NULL)) ; /* keep all nonzero entries */
}

34
deps/csparse/Source/cs_dupl.c vendored Normal file
View file

@ -0,0 +1,34 @@
#include "cs.h"
/* remove duplicate entries from A */
csi cs_dupl (cs *A)
{
csi 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 (csi)) ; /* 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 */
}

13
deps/csparse/Source/cs_entry.c vendored Normal file
View file

@ -0,0 +1,13 @@
#include "cs.h"
/* add an entry to a triplet matrix; return 1 if ok, 0 otherwise */
csi cs_entry (cs *T, csi i, csi 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) ;
}

23
deps/csparse/Source/cs_ereach.c vendored Normal file
View file

@ -0,0 +1,23 @@
#include "cs.h"
/* find nonzero pattern of Cholesky L(k,1:k-1) using etree and triu(A(:,k)) */
csi cs_ereach (const cs *A, csi k, const csi *parent, csi *s, csi *w)
{
csi 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,:)*/
}

30
deps/csparse/Source/cs_etree.c vendored Normal file
View file

@ -0,0 +1,30 @@
#include "cs.h"
/* compute the etree of A (using triu(A), or A'A without forming A'A */
csi *cs_etree (const cs *A, csi ata)
{
csi 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 (csi)) ; /* allocate result */
w = cs_malloc (n + (ata ? m : 0), sizeof (csi)) ; /* 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)) ;
}

25
deps/csparse/Source/cs_fkeep.c vendored Normal file
View file

@ -0,0 +1,25 @@
#include "cs.h"
/* drop entries for which fkeep(A(i,j)) is false; return nz if OK, else -1 */
csi cs_fkeep (cs *A, csi (*fkeep) (csi, csi, double, void *), void *other)
{
csi 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) ;
}

17
deps/csparse/Source/cs_gaxpy.c vendored Normal file
View file

@ -0,0 +1,17 @@
#include "cs.h"
/* y = A*x+y */
csi cs_gaxpy (const cs *A, const double *x, double *y)
{
csi 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) ;
}

19
deps/csparse/Source/cs_happly.c vendored Normal file
View file

@ -0,0 +1,19 @@
#include "cs.h"
/* apply the ith Householder vector to x */
csi cs_happly (const cs *V, csi i, double beta, double *x)
{
csi 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) ;
}

23
deps/csparse/Source/cs_house.c vendored Normal file
View file

@ -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, csi n)
{
double s, sigma = 0 ;
csi 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) ;
}

9
deps/csparse/Source/cs_ipvec.c vendored Normal file
View file

@ -0,0 +1,9 @@
#include "cs.h"
/* x(p) = b, for dense vectors x and b; p=NULL denotes identity */
csi cs_ipvec (const csi *p, const double *b, double *x, csi n)
{
csi k ;
if (!x || !b) return (0) ; /* check inputs */
for (k = 0 ; k < n ; k++) x [p ? p [k] : k] = b [k] ;
return (1) ;
}

22
deps/csparse/Source/cs_leaf.c vendored Normal file
View file

@ -0,0 +1,22 @@
#include "cs.h"
/* consider A(i,j), node j in ith row subtree and return lca(jprev,j) */
csi cs_leaf (csi i, csi j, const csi *first, csi *maxfirst, csi *prevleaf,
csi *ancestor, csi *jleaf)
{
csi 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) */
}

15
deps/csparse/Source/cs_load.c vendored Normal file
View file

@ -0,0 +1,15 @@
#include "cs.h"
/* load a triplet matrix from a file */
cs *cs_load (FILE *f)
{
double i, j ; /* use double for integers to avoid csi conflicts */
double x ;
cs *T ;
if (!f) return (NULL) ; /* check inputs */
T = cs_spalloc (0, 0, 1, 1, 1) ; /* allocate result */
while (fscanf (f, "%lg %lg %lg\n", &i, &j, &x) == 3)
{
if (!cs_entry (T, (csi) i, (csi) j, x)) return (cs_spfree (T)) ;
}
return (T) ;
}

18
deps/csparse/Source/cs_lsolve.c vendored Normal file
View file

@ -0,0 +1,18 @@
#include "cs.h"
/* solve Lx=b where x and b are dense. x=b on input, solution on output. */
csi cs_lsolve (const cs *L, double *x)
{
csi 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) ;
}

18
deps/csparse/Source/cs_ltsolve.c vendored Normal file
View file

@ -0,0 +1,18 @@
#include "cs.h"
/* solve L'x=b where x and b are dense. x=b on input, solution on output. */
csi cs_ltsolve (const cs *L, double *x)
{
csi 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) ;
}

86
deps/csparse/Source/cs_lu.c vendored Normal file
View file

@ -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 ;
csi *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 (csi)) ; /* get csi 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 (csi)) ; /* 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 */
}

26
deps/csparse/Source/cs_lusol.c vendored Normal file
View file

@ -0,0 +1,26 @@
#include "cs.h"
/* x=A\b where A is unsymmetric; b overwritten with solution */
csi cs_lusol (csi order, const cs *A, double *b, double tol)
{
double *x ;
css *S ;
csn *N ;
csi 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) ;
}

35
deps/csparse/Source/cs_malloc.c vendored Normal file
View file

@ -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 (csi n, size_t size)
{
return (malloc (CS_MAX (n,1) * size)) ;
}
/* wrapper for calloc */
void *cs_calloc (csi 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, csi n, size_t size, csi *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 */
}

92
deps/csparse/Source/cs_maxtrans.c vendored Normal file
View file

@ -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 (csi k, const cs *A, csi *jmatch, csi *cheap, csi *w,
csi *js, csi *is, csi *ps)
{
csi 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 */
csi *cs_maxtrans (const cs *A, csi seed) /*[jmatch [0..m-1]; imatch [0..n-1]]*/
{
csi 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 (csi)) ; /* 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 (csi)) ; /* 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)) ;
}

35
deps/csparse/Source/cs_multiply.c vendored Normal file
View file

@ -0,0 +1,35 @@
#include "cs.h"
/* C = A*B */
cs *cs_multiply (const cs *A, const cs *B)
{
csi 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 (csi)) ; /* 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 */
}

15
deps/csparse/Source/cs_norm.c vendored Normal file
View file

@ -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)
{
csi 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) ;
}

25
deps/csparse/Source/cs_permute.c vendored Normal file
View file

@ -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 csi *pinv, const csi *q, csi values)
{
csi 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)) ;
}

11
deps/csparse/Source/cs_pinv.c vendored Normal file
View file

@ -0,0 +1,11 @@
#include "cs.h"
/* pinv = p', or p = pinv' */
csi *cs_pinv (csi const *p, csi n)
{
csi k, *pinv ;
if (!p) return (NULL) ; /* p = NULL denotes identity */
pinv = cs_malloc (n, sizeof (csi)) ; /* 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 */
}

24
deps/csparse/Source/cs_post.c vendored Normal file
View file

@ -0,0 +1,24 @@
#include "cs.h"
/* post order a forest */
csi *cs_post (const csi *parent, csi n)
{
csi j, k = 0, *post, *w, *head, *next, *stack ;
if (!parent) return (NULL) ; /* check inputs */
post = cs_malloc (n, sizeof (csi)) ; /* allocate result */
w = cs_malloc (3*n, sizeof (csi)) ; /* 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 */
}

39
deps/csparse/Source/cs_print.c vendored Normal file
View file

@ -0,0 +1,39 @@
#include "cs.h"
/* print a sparse matrix; use %g for integers to avoid differences with csi */
csi cs_print (const cs *A, csi brief)
{
csi 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 ("%g-by-%g, nzmax: %g nnz: %g, 1-norm: %g\n", (double) m,
(double) n, (double) nzmax, (double) (Ap [n]), cs_norm (A)) ;
for (j = 0 ; j < n ; j++)
{
printf (" col %g : locations %g to %g\n", (double) j,
(double) (Ap [j]), (double) (Ap [j+1]-1)) ;
for (p = Ap [j] ; p < Ap [j+1] ; p++)
{
printf (" %g : %g\n", (double) (Ai [p]), Ax ? Ax [p] : 1) ;
if (brief && p > 20) { printf (" ...\n") ; return (1) ; }
}
}
}
else
{
printf ("triplet: %g-by-%g, nzmax: %g nnz: %g\n", (double) m,
(double) n, (double) nzmax, (double) nz) ;
for (p = 0 ; p < nz ; p++)
{
printf (" %g %g : %g\n", (double) (Ai [p]), (double) (Ap [p]),
Ax ? Ax [p] : 1) ;
if (brief && p > 20) { printf (" ...\n") ; return (1) ; }
}
}
return (1) ;
}

9
deps/csparse/Source/cs_pvec.c vendored Normal file
View file

@ -0,0 +1,9 @@
#include "cs.h"
/* x = b(p), for dense vectors x and b; p=NULL denotes identity */
csi cs_pvec (const csi *p, const double *b, double *x, csi n)
{
csi k ;
if (!x || !b) return (0) ; /* check inputs */
for (k = 0 ; k < n ; k++) x [k] = b [p ? p [k] : k] ;
return (1) ;
}

73
deps/csparse/Source/cs_qr.c vendored Normal file
View file

@ -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 ;
csi 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 (csi)) ; /* get csi 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 */
}

53
deps/csparse/Source/cs_qrsol.c vendored Normal file
View file

@ -0,0 +1,53 @@
#include "cs.h"
/* x=A\b where A can be rectangular; b overwritten with solution */
csi cs_qrsol (csi order, const cs *A, double *b)
{
double *x ;
css *S ;
csn *N ;
cs *AT = NULL ;
csi 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) ;
}

22
deps/csparse/Source/cs_randperm.c vendored Normal file
View file

@ -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. */
csi *cs_randperm (csi n, csi seed)
{
csi *p, k, j, t ;
if (seed == 0) return (NULL) ; /* return p = NULL (identity) */
p = cs_malloc (n, sizeof (csi)) ; /* 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 integer in range k to n-1 */
t = p [j] ; /* swap p[k] and p[j] */
p [j] = p [k] ;
p [k] = t ;
}
return (p) ;
}

19
deps/csparse/Source/cs_reach.c vendored Normal file
View file

@ -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 */
csi cs_reach (cs *G, const cs *B, csi k, csi *xi, const csi *pinv)
{
csi 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) ;
}

22
deps/csparse/Source/cs_scatter.c vendored Normal file
View file

@ -0,0 +1,22 @@
#include "cs.h"
/* x = x + beta * A(:,j), where x is a dense vector and A(:,j) is sparse */
csi cs_scatter (const cs *A, csi j, double beta, csi *w, double *x, csi mark,
cs *C, csi nz)
{
csi 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) ;
}

41
deps/csparse/Source/cs_scc.c vendored Normal file
View file

@ -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 */
{
csi 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 (csi)) ; /* 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)) ;
}

26
deps/csparse/Source/cs_schol.c vendored Normal file
View file

@ -0,0 +1,26 @@
#include "cs.h"
/* ordering and symbolic analysis for a Cholesky factorization */
css *cs_schol (csi order, const cs *A)
{
csi 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 (csi)) ; /* 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)) ;
}

28
deps/csparse/Source/cs_spsolve.c vendored Normal file
View file

@ -0,0 +1,28 @@
#include "cs.h"
/* solve Gx=b(:,k), where G is either upper (lo=0) or lower (lo=1) triangular */
csi cs_spsolve (cs *G, const cs *B, csi k, csi *xi, double *x, const csi *pinv,
csi lo)
{
csi 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 */
}

87
deps/csparse/Source/cs_sqr.c vendored Normal file
View file

@ -0,0 +1,87 @@
#include "cs.h"
/* compute nnz(V) = S->lnz, S->pinv, S->leftmost, S->m2 from A and S->parent */
static csi cs_vcount (const cs *A, css *S)
{
csi 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 (csi)) ; /* allocate pinv, */
S->leftmost = leftmost = cs_malloc (m, sizeof (csi)) ; /* and leftmost */
w = cs_malloc (m+3*n, sizeof (csi)) ; /* 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 (csi order, const cs *A, csi qr)
{
csi 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] ;
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 */
}

39
deps/csparse/Source/cs_symperm.c vendored Normal file
View file

@ -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 csi *pinv, csi values)
{
csi 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 (csi)) ; /* 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 */
}

24
deps/csparse/Source/cs_tdfs.c vendored Normal file
View file

@ -0,0 +1,24 @@
#include "cs.h"
/* depth-first search and postorder of a tree rooted at node j */
csi cs_tdfs (csi j, csi k, csi *head, const csi *next, csi *post, csi *stack)
{
csi 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) ;
}

25
deps/csparse/Source/cs_transpose.c vendored Normal file
View file

@ -0,0 +1,25 @@
#include "cs.h"
/* C = A' */
cs *cs_transpose (const cs *A, csi values)
{
csi 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 (csi)) ; /* 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 */
}

37
deps/csparse/Source/cs_updown.c vendored Normal file
View file

@ -0,0 +1,37 @@
#include "cs.h"
/* sparse Cholesky update/downdate, L*L' + sigma*w*w' (sigma = +1 or -1) */
csi cs_updown (cs *L, csi sigma, const cs *C, const csi *parent)
{
csi 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) ;
}

18
deps/csparse/Source/cs_usolve.c vendored Normal file
View file

@ -0,0 +1,18 @@
#include "cs.h"
/* solve Ux=b where x and b are dense. x=b on input, solution on output. */
csi cs_usolve (const cs *U, double *x)
{
csi 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) ;
}

119
deps/csparse/Source/cs_util.c vendored Normal file
View file

@ -0,0 +1,119 @@
#include "cs.h"
/* allocate a sparse matrix (triplet form or compressed-column form) */
cs *cs_spalloc (csi m, csi n, csi nzmax, csi values, csi 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 (csi)) ;
A->i = cs_malloc (nzmax, sizeof (csi)) ;
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 */
csi cs_sprealloc (cs *A, csi nzmax)
{
csi 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 (csi), &oki) ;
if (CS_TRIPLET (A)) A->p = cs_realloc (A->p, nzmax, sizeof (csi), &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 *) 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 ((csn *) 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 ((css *) cs_free (S)) ; /* free the css struct and return NULL */
}
/* allocate a cs_dmperm or cs_scc result */
csd *cs_dalloc (csi m, csi n)
{
csd *D ;
D = cs_calloc (1, sizeof (csd)) ;
if (!D) return (NULL) ;
D->p = cs_malloc (m, sizeof (csi)) ;
D->r = cs_malloc (m+6, sizeof (csi)) ;
D->q = cs_malloc (n, sizeof (csi)) ;
D->s = cs_malloc (n+6, sizeof (csi)) ;
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 ((csd *) cs_free (D)) ; /* free the csd struct and return NULL */
}
/* free workspace and return a sparse matrix result */
cs *cs_done (cs *C, void *w, void *x, csi 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 csi array result */
csi *cs_idone (csi *p, cs *C, void *w, csi ok)
{
cs_spfree (C) ; /* free temporary matrix */
cs_free (w) ; /* free workspace */
return (ok ? p : (csi *) cs_free (p)) ; /* return result, or free it */
}
/* free workspace and return a numeric factorization (Cholesky, LU, or QR) */
csn *cs_ndone (csn *N, cs *C, void *w, void *x, csi 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, csi 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 */
}

18
deps/csparse/Source/cs_utsolve.c vendored Normal file
View file

@ -0,0 +1,18 @@
#include "cs.h"
/* solve U'x=b where x and b are dense. x=b on input, solution on output. */
csi cs_utsolve (const cs *U, double *x)
{
csi 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) ;
}

102
deps/csparse/Tcov/Makefile vendored Normal file
View file

@ -0,0 +1,102 @@
# To run with valgrind:
V =
# V = valgrind -q
# Linux test coverage
CC = gcc
CFLAGS = -O -g -fprofile-arcs -ftest-coverage \
-Wall -W -Wshadow -Wmissing-prototypes -Wstrict-prototypes \
-Wredundant-decls -Wnested-externs -Wdisabled-optimization -ansi \
-Wno-unused-parameter -I../Include -I../Demo
run: all run1 runbook run3 runtest
./covall
all: cs_demo1 cs_demo2 cs_demo3 cstcov_test
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_transpose.o cs_compress.o cs_usolve.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 \
cstcov_malloc_test.o cs_utsolve.o cs_reach.o cs_spsolve.o cs_ereach.o \
cs_leaf.o cs_randperm.o
$(CS): ../Include/cs.h cstcov_malloc_test.h
.PRECIOUS: cs_%.c cs_dem%.c
cs_dem%.c:
- ln -s ../Demo/$@
cs_%.c:
- ln -s ../Source/$@
cs_demo1: $(CS) cs_demo1.c
$(CC) $(CFLAGS) -o cs_demo1 cs_demo1.c $(CS) -lm
cs_demo2: $(CS) cs_demo2.c cs_demo.c
$(CC) $(CFLAGS) -o cs_demo2 cs_demo2.c cs_demo.c $(CS) -lm
cs_demo3: $(CS) cs_demo3.c cs_demo.c
$(CC) $(CFLAGS) -o cs_demo3 cs_demo3.c cs_demo.c $(CS) -lm
cstcov_test: $(CS) cstcov_test.c cs_demo.c
$(CC) $(CFLAGS) -o cstcov_test cstcov_test.c cs_demo.c $(CS) -lm
# tiny tests
run1: all
- $(V) ./cs_demo1 < nil
- $(V) ./cs_demo1 < zero
- $(V) ./cs_demo2 < nil
- $(V) ./cs_demo2 < zero
- $(V) ./cs_demo3 < nil
# test coverage for book:
runbook: all
- $(V) ./cs_demo1 < ../Matrix/t1
- $(V) ./cs_demo2 < ../Matrix/bcsstk01
- $(V) ./cs_demo2 < ../Matrix/fs_183_1
- $(V) ./cs_demo2 < ../Matrix/mbeacxc
- $(V) ./cs_demo2 < ../Matrix/west0067
- $(V) ./cs_demo2 < ../Matrix/lp_afiro
- $(V) ./cs_demo3 < ../Matrix/bcsstk16
# other tests
run3: all
- $(V) ./cs_demo2 < ../Matrix/t1
- $(V) ./cs_demo2 < ../Matrix/ash219
- $(V) ./cs_demo3 < ../Matrix/bcsstk01
- $(V) ./cs_demo2 < ../Matrix/bcsstk16
- $(V) ./cs_demo2 < ../Matrix/ibm32a
- $(V) ./cs_demo2 < ../Matrix/ibm32b
# exhaustive memory tests
runtest: all
- $(V) ./cstcov_test nil > test_nil.out
- $(V) ./cstcov_test zero > test_zero.out
- $(V) ./cstcov_test ../Matrix/t1 > test_t1.out
- $(V) ./cstcov_test ../Matrix/bcsstk01 > test_k1.out
- $(V) ./cstcov_test ../Matrix/fs_183_1 > test_fs.out
- $(V) ./cstcov_test ../Matrix/west0067 > test_we.out
- $(V) ./cstcov_test ../Matrix/ash219 > test_ash.out
- $(V) ./cstcov_test ../Matrix/lp_afiro > test_afiro.out
readhb: readhb.f
f77 -o readhb readhb.f
readhb.f:
- ln -s readhb.f
clean:
- $(RM) *.o *.bbg *.da *.gcov *.gcda *.gcno
purge: distclean
# remove everything for distribution, including all symbolic links
distclean: clean
- $(RM) cs_demo1 cs_demo2 readhb *.out *.a cs_demo3 cstcov_test cov.sort
- $(RM) -r cs_*.c *.dSYM

15
deps/csparse/Tcov/README.txt vendored Normal file
View file

@ -0,0 +1,15 @@
CSparse/Tcov: comprehensive test coverage for CSparse. Requires Linux.
Type "make" to compile, and then "make run" to run the tests.
The test coverage is in cover.out. The test output is
printed on stdout, except for cs_test (which prints its output in various
*.out files).
If the test is successful, the last line printed should be
"statements not yet tested: 0", and all printed residuals should be small.
Note that you will get warnings about unused parameters for some functions.
These warnings can be safely ignored. They are parameters for functions that
are passed to cs_fkeep, and all functions used in this manner must have the
same calling sequence, even if some of the parameters are not used.
Timothy A. Davis, http://www.suitesparse.com

17
deps/csparse/Tcov/cov.awk vendored Normal file
View file

@ -0,0 +1,17 @@
/cannot/
/function/ { f = $8 }
/file/ { f = $8 }
/lines/ {
k = match ($1, "%") ;
p = substr ($1, 1, k-1) ;
if ((p+0) != 100)
{
printf "%8s %s\n", p, f
}
}

7
deps/csparse/Tcov/covall vendored Executable file
View file

@ -0,0 +1,7 @@
#!/bin/csh
./gcovs cs*.c |& awk -f cov.awk | sort -n > cov.out
sort -n cov.out > cov.sort
./covs > covs.out
echo -n "statments not yet tested: "
grep "#####" *gcov | wc -l
./cover *v > cover.out

7
deps/csparse/Tcov/covall.linux vendored Executable file
View file

@ -0,0 +1,7 @@
#!/bin/csh
./gcovs cs*.c |& awk -f cov.awk | sort -n > cov.out
sort -n cov.out > cov.sort
./covs > covs.out
echo -n "statments not yet tested: "
grep "#####" *gcov | wc -l
./cover *v > cover.out

5
deps/csparse/Tcov/covall.sol vendored Executable file
View file

@ -0,0 +1,5 @@
#!/bin/csh
tcov -x cm.profile cs*.c >& /dev/null
echo -n "statments not yet tested: "
./covs > covs.out
grep "#####" *tcov | wc -l

9
deps/csparse/Tcov/cover vendored Executable file
View file

@ -0,0 +1,9 @@
#!/bin/csh
# usage: cover files
echo '================================================================='
foreach file ($argv[1-])
echo $file
echo '================================================================='
grep "#####" -A5 -B5 $file
echo '================================================================='
end

7
deps/csparse/Tcov/covs vendored Executable file
View file

@ -0,0 +1,7 @@
#!/bin/csh
echo '================================================================='
foreach file (*.?cov)
echo $file
grep "#####" $file
echo '================================================================='
end

34
deps/csparse/Tcov/cstcov_malloc_test.c vendored Normal file
View file

@ -0,0 +1,34 @@
#include "cstcov_malloc_test.h"
csi malloc_count = INT_MAX ;
/* wrapper for malloc */
void *cs_malloc (csi n, size_t size)
{
if (--malloc_count < 0) return (NULL) ; /* pretend to fail */
return (malloc (CS_MAX (n,1) * size)) ;
}
/* wrapper for calloc */
void *cs_calloc (csi n, size_t size)
{
if (--malloc_count < 0) return (NULL) ; /* pretend to fail */
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, csi n, size_t size, csi *ok)
{
void *pnew ;
*ok = 0 ;
if (--malloc_count < 0) return (p) ; /* pretend to fail */
pnew = realloc (p, CS_MAX (n,1) * size) ; /* realloc the block */
*ok = (pnew != NULL) ;
return ((*ok) ? pnew : p) ; /* return original p if failure */
}

View file

@ -0,0 +1,5 @@
#include "cs.h"
#ifndef EXTERN
#define EXTERN extern
#endif
EXTERN csi malloc_count ;

30
deps/csparse/Tcov/cstcov_test.c vendored Normal file
View file

@ -0,0 +1,30 @@
#include "cs_demo.h"
/* cs_test: read a matrix and run cs_demo2 and cs_demo3, using malloc tests. */
#include "cstcov_malloc_test.h"
int main (int argc, char **argv)
{
FILE *f ;
problem *Prob ;
int trials, ok, demo ;
if (argc < 2) return (-1) ;
printf ("cs_test, file: %s\n", argv [1]) ;
for (demo = 2 ; demo <= 3 ; demo++)
{
printf ("demo: %g\n", (double) demo) ;
for (trials = 0 ; trials < 4000 ; trials++)
{
malloc_count = trials ;
f = fopen (argv [1], "r") ;
if (!f) return (-1) ;
Prob = get_problem (f, (demo == 2) ? 1e-14 : 0) ;
fclose (f) ;
if (Prob) ok = (demo == 2) ? demo2 (Prob) : demo3 (Prob) ;
free_problem (Prob) ;
if (malloc_count > 0) break ;
}
printf ("demo %g # trials: %g\n", (double) demo, (double) trials) ;
}
return (0) ;
}

6
deps/csparse/Tcov/gcovs vendored Executable file
View file

@ -0,0 +1,6 @@
# usage: gcovs files
echo '================================================================='
foreach file ($argv[1-])
gcov -f $file
echo '================================================================='
end

0
deps/csparse/Tcov/nil vendored Normal file
View file

1
deps/csparse/Tcov/zero vendored Normal file
View file

@ -0,0 +1 @@
0 0 0

1
deps/lily-test vendored Submodule

@ -0,0 +1 @@
Subproject commit 6833afddd9fa9b76834cd18ac604340fcaafd47c

View file

@ -1,10 +0,0 @@
{
"devDependencies": {
"jest": "^29.5.0"
},
"scripts": {
"test": "NODE_OPTIONS=--experimental-vm-modules jest",
"wtest": "test.bat"
},
"type": "module"
}

17
public/index.html Normal file
View file

@ -0,0 +1,17 @@
<!doctype html>
<html>
<head>
<meta charset="utf-8">
<title>nerine</title>
<script src="main.js"></script>
</head>
<body>
<h1>nerine</h1>
<h2 id="status-header">epoch 0, step 0</h2>
<h3 id="stats-header"></h3>
<form>
<input type="file" onchange="load_start(this)"></input>
</form>
<canvas id="render-canvas" width="500" height="500"></canvas>
</body>
</html>

945
public/main.js Normal file
View file

@ -0,0 +1,945 @@
// ================================================================
//
// util
//
// ================================================================
let random_seed = 0x12345678;
function random() {
random_seed = (random_seed * 1664525 + 1013904223) | 0;
const value = (random_seed >>> 0) / 0x100000000;
return value;
}
// modulo value for unsigned integers
const UINT_MOD = 0x10000;
const UINT_MAX = 0xffff;
function rand_uint() {
return Math.floor(UINT_MOD * random());
}
// nonnegative modulo
function mod(n, m) {
const k = n % m;
if (k < 0) {
return k+m;
} else {
return k;
}
}
function pos_mod([x, y], [w, h]) {
return [ mod(x, w), mod(y, h) ];
}
function distance2(p1, p2) {
const [x1, y1] = p1;
const [x2, y2] = p2;
const dx = x2 - x1;
const dy = y2 - y1;
return (dx*dx) + (dy*dy);
}
function rand_idx(arr) {
return Math.floor(arr.length * random());
}
function rand_choice(arr) {
return arr[rand_idx(arr)];
}
function to_shuffled(arr) {
const copy = [...arr];
for (let i=0; i<arr.length-1; i++) {
const j = Math.floor((arr.length-i) * random()) + i;
const [x, y] = [copy[i], copy[j]];
copy[i] = y;
copy[j] = x;
}
return copy;
}
function get_step(direction) {
switch (direction) {
case 'north':
return [0, 1];
case 'south':
return [0, -1];
case 'east':
return [1, 0];
case 'west':
return [-1, 0];
default:
throw new Error(`unknown direction: ${direction}`);
}
}
function rotate([dx, dy], direction) {
switch (direction) {
case 'north':
return [dx, dy];
case 'south':
return [-dx, -dy];
case 'east':
return [dy, -dx];
case 'west':
return [-dy, dx];
default:
throw new Error(`unknown direction: ${direction}`);
}
}
// ================================================================
//
// mind
//
// ================================================================
function mind_get_next_instruction(mind) {
const { pc, memory } = mind;
const operand_idx = (pc + 1) % memory.length;
return [ memory[pc], memory[operand_idx] ];
}
function mind_execute_instruction(mind, sense, act, state) {
const { memory } = mind;
const [ instruction, operand ] = mind_get_next_instruction(mind);
const opcode = instruction & 0xf; // bottom 4 bits
let result = [];
switch (opcode) {
case 0: // add
mind.acc += memory[operand % memory.length];
mind.acc %= UINT_MOD;
mind.pc += 2;
break;
case 1: // bit shift left
mind.acc <<= memory[operand % memory.length];
mind.acc %= UINT_MOD;
mind.pc += 2;
break;
case 2: // bit shift right
mind.acc >>= memory[operand % memory.length];
mind.acc %= UINT_MOD;
mind.pc += 2;
break;
case 3: // bitwise AND
mind.acc &= memory[operand % memory.length];
mind.pc += 2;
break;
case 4: // bitwise XOR
mind.acc ^= memory[operand % memory.length];
mind.pc += 2;
break;
case 5: // store
memory[operand % memory.length] = mind.acc;
mind.pc += 2;
break;
case 6: // store pc
memory[operand % memory.length] = mind.pc;
mind.pc += 2;
break;
case 7: // copy to cell[acc]
memory[mind.acc % memory.length] = memory[operand % memory.length];
mind.pc += 2;
break;
case 8: // copy from cell[acc]
memory[operand % memory.length] = memory[mind.acc % memory.length];
mind.pc += 2;
break;
case 9: // load
mind.acc = memory[operand % memory.length];
mind.pc += 2;
break;
case 10: // loada
mind.acc = memory[mind.acc % memory.length];
mind.pc += 2;
break;
case 11: // jump
mind.pc = memory[operand % memory.length];
break;
case 12: // jnz
if (mind.acc !== 0) {
mind.pc = memory[operand % memory.length];
} else {
mind.pc += 2;
}
break;
case 13: // jez
if (mind.acc === 0) {
mind.pc = memory[operand % memory.length];
} else {
mind.pc += 2;
}
break;
case 14: // sense
mind.acc = sense(state, operand);
mind.pc += 2;
break;
case 15: // act
result = act(state, operand);
mind.pc += 2;
break;
default:
throw new Error(`invalid opcode: ${opcode}`);
}
mind.pc %= memory.length;
return result;
}
// ================================================================
//
// genome
//
// ================================================================
function genome_instantiate(genome) {
const memory = genome.pairs.reduce(
(acc, [addr, value]) => {
acc[addr % acc.length] = value;
return acc;
},
[...Array(genome.size)].map(() => rand_uint()),
);
return {
memory,
pc: 0,
acc: rand_uint(),
};
}
function genome_mut_append(genome) {
return {
...genome,
pairs: [...genome.pairs, [rand_uint(), rand_uint()]],
}
}
function genome_mut_dup(genome) {
return {
...genome,
pairs: [ ...genome.pairs, genome.pairs[rand_idx(genome.pairs)] ],
};
}
function genome_mut_delete(genome) {
return {
...genome,
pairs: genome.pairs.toSpliced(rand_idx(genome.pairs), 1),
};
}
function genome_mut_addr(genome) {
const idx = rand_idx(genome.pairs);
const [ addr, value ] = genome.pairs[idx];
return {
...genome,
pairs: genome.pairs.toSpliced(idx, 1, [ addr & rand_uint(), value ]),
}
}
function genome_mut_value(genome) {
const idx = rand_idx(genome.pairs);
const [ addr, value ] = genome.pairs[idx];
return {
...genome,
pairs: genome.pairs.toSpliced(idx, 1, [ addr, value & rand_uint() ]),
}
}
function genome_mut_size(genome) {
return {
...genome,
size: Math.max(genome.size + random() > 0.5 ? 1 : -1, 16),
};
}
function genome_mut(genome) {
const mutations = [
genome_mut_append,
genome_mut_dup,
genome_mut_delete,
genome_mut_addr,
genome_mut_value,
genome_mut_size,
];
const idx = rand_idx(mutations);
return mutations[idx](genome);
}
// ================================================================
//
// simulation
//
// ================================================================
function world_create(width, height) {
const world = {
width,
height,
next_id: 0,
creatures: {},
creature_pos: {},
pos_creature: {},
time: 0,
sound_frames: [],
current_sound_frame: [],
energy_sources: [...Array(4)].map(() => pos_mod([rand_uint(), rand_uint()], [ width, height ])),
};
world.act = make_act(world);
world.sense = make_sense(world);
return world;
}
// check if a position is currently occupied
function world_occupied(world, pos) {
return world.pos_creature[pos] !== undefined;
}
// add a new creature to the world
// return the id of the inserted creature on success, or null on failure
function world_creature_insert(world, creature, pos) {
pos = pos_mod(pos, [world.width, world.height]);
if (!world_occupied(world, pos)) {
const id = world.next_id;
world.next_id += 1;
world.creatures[id] = creature;
world.creature_pos[id] = pos;
world.pos_creature[pos] = id;
return id;
} else {
return null;
}
}
function world_creature_remove(world, id) {
const pos = world.creature_pos[id];
delete world.creatures[id];
delete world.creature_pos[id];
delete world.pos_creature[pos];
}
// attempt to move a creature
// do nothing if the target position is blocked
function world_creature_move(world, id, new_pos) {
new_pos = pos_mod(new_pos, [world.width, world.height]);
if (!world_occupied(world, new_pos)) {
const old_pos = world.creature_pos[id];
delete world.pos_creature[old_pos];
world.pos_creature[new_pos] = id;
world.creature_pos[id] = new_pos;
} else {
// blocked, do nothing
}
}
const MAX_SOUND_AGE = 64;
function world_step(world) {
// if (random() * world.width * world.height < world.creatures.length) {
// // delete a random creature if the population is getting too high
// world_creature_remove(world, rand_choice([...Object.keys(world.creatures)]));
// }
[...Object.values(world.creatures)].forEach(x => creature_step(x, world));
world.time += 1;
world.sound_frames = [
world.current_sound_frame,
...world.sound_frames,
].toSpliced(MAX_SOUND_AGE);
world.current_sound_frame = [];
}
// ================================================================
//
// creatures
//
// ================================================================
function creature_create(genome) {
const mind = genome_instantiate(genome);
return {
id: null,
genome,
mind,
frozen: 'unfrozen',
direction: rand_choice(['north', 'south', 'east', 'west']),
local_energy: 0,
foreign_energy: 0,
};
}
function creature_step(creature, world) {
if (creature.frozen !== 'frozen') {
const pos = world.creature_pos[creature.id];
const energy_level = world.energy_sources.reduce(
(acc, source_pos) => acc + (10/Math.sqrt(distance2(pos, source_pos))),
0,
);
if (random() < energy_level) {
creature.local_energy += 1;
}
}
mind_execute_instruction(creature.mind, world.sense, world.act, creature);
}
function creature_forward(creature, world) {
const [dx, dy] = get_step(creature.direction);
const [x0, y0] = world.creature_pos[creature.id];
return pos_mod([x0+dx, y0+dy], [world.width, world.height]);
}
function creature_move(creature, world) {
const target = creature_forward(creature, world);
world_creature_move(world, creature.id, target);
}
const ENERGY_CONV = 16;
// try to spend `amount` of (foreign) energy
// return true if successful and false if not enough energy was available
function creature_use_energy(creature, amount) {
const total_energy = Math.floor(creature.local_energy / ENERGY_CONV) + creature.foreign_energy;
if (total_energy < amount) {
// not enough energy :c
return false;
}
if (creature.foreign_energy >= amount) {
// foreign energy covers it!
creature.foreign_energy -= amount;
return true;
}
// split between foreign and local energy
const remaining = amount - creature.foreign_energy;
creature.foreign_energy = 0;
creature.local_energy -= ENERGY_CONV * remaining;
return true;
}
// --===== senses =====--
// detect the presence of another creature
function creature_sense_presence(creature, world, pos) {
creature.mind.acc = world_occupied(world, pos) ? 1 : 0;
}
// detect if another creature is frozen (or at least pretending to be)
function creature_sense_frozen(creature, world, pos) {
if (!world_occupied(world, pos)) {
// no creature, so not frozen
creature.mind.acc = 0;
} else {
const target_id = world.pos_creature[pos];
const target = world.creatures[target_id];
creature.mind.acc = target.frozen === 'frozen' || target.frozen === 'pretend' ? 1 : 0;
}
}
// read another creature's mind
function creature_read_mind(creature, world, pos, pc) {
if (!world_occupied(world, pos)) {
creature.mind.acc = 0;
} else {
const id = world.pos_creature[pos];
const target = world.creatures[id];
creature.mind.acc = target.mind.memory[pc];
}
}
// detect sounds as they were in the recent past
function creature_listen(creature, world, time) {
const cpos = world.creature_pos[creature.id];
const frame = world.sound_frames[time % MAX_SOUND_AGE] || [];
const value = frame.reduce(
(acc, source) => {
const { pos, value } = source;
return acc + Math.abs(Math.floor(value/distance2(cpos, pos)));
},
0
);
creature.mind.acc = value > UINT_MAX ? UINT_MAX : value;
}
function creature_check_time(creature, world) {
creature.mind.acc = world.time;
}
// check how much local energy we have
function creature_check_local(creature) {
creature.mind.acc = mod(creature.local_energy, UINT_MOD);
}
// check how much foreign energy we have
function creature_check_foreign(creature) {
creature.mind.acc = mod(creature.foreign_energy, UINT_MOD);
}
function make_sense(world) {
return (creature, value) => {
const choice = value & 0x7;
const rest = value >> 3;
const dx = (rest & 0xf) - 15;
const dy = ((rest >> 4) & 0xf) - 15;
const target_pos = rotate([dx, dy], creature.direction);
switch (choice) {
case 0:
case 1:
creature_sense_presence(creature, world, target_pos);
break;
case 2:
creature_sense_frozen(creature, world, target_pos);
break;
case 3:
creature_read_mind(creature, world, target_pos, (rest >> 8) & 0x1f);
break;
case 4:
creature_listen(creature, world, rest);
break;
case 5:
creature_check_time(creature, world);
break;
case 6:
creature_check_local(creature);
break;
case 7:
creature_check_foreign(creature);
break;
default:
throw new Error(`impossible sense choice: ${choice}`);
}
}
}
// --===== actions =====--
// create offspring
function creature_reproduce(creature, world) {
if (creature.frozen === 'frozen') { return; }
const pos = creature_forward(creature, world);
if (world_occupied(world, pos)) {
// blocked!
return;
}
const energy_cost = creature.genome.size + creature.genome.pairs.length;
if (!creature_use_energy(creature, energy_cost)) {
// not enough energy
return;
}
const max_mutations = Math.max(0.1 * energy_cost, 1);
const new_genome = [...Array(Math.floor(max_mutations * random()))]
.reduce(genome_mut, creature.genome)
const child = creature_create(new_genome);
child.id = world_creature_insert(world, child, pos);
}
// create a sound
function creature_speak(creature, world, value) {
world.current_sound_frame.push({ pos: world.creature_pos[creature.id], value });
}
// freeze a target
function creature_freeze(creature, world) {
if (creature.frozen === 'frozen') { return; }
const pos = creature_forward(creature, world);
if (!world_occupied(world, pos)) {
// no target...
return;
} else {
const target_id = world.pos_creature[pos];
const target = world.creatures[target_id];
target.frozen = 'frozen';
}
}
// unfreeze a target
function creature_unfreeze(creature, world) {
if (creature.frozen === 'frozen') { return; }
const pos = creature_forward(creature, world);
if (!world_occupied(world, pos)) {
// no target...
return;
} else {
const target_id = world.pos_creature[pos];
const target = world.creatures[target_id];
if (target.frozen === 'frozen') {
target.frozen === 'unfrozen';
}
}
}
// set pretend frozen
function creature_set_pretend(creature, world, value) {
if (creature.frozen === 'frozen') { return; }
if (value === 0) {
creature.frozen = 'unfrozen';
} else {
creature.frozen = 'pretend';
}
}
// give energy
function creature_give(creature, world, value) {
if (creature.frozen === 'frozen') { return; }
const amount = Math.min(creature.local_energy, value);
const pos = creature_forward(creature, world);
if (!world_occupied(world, pos)) {
// no target...
return;
} else {
const target_id = world.pos_creature[pos];
const target = world.creatures[target_id];
target.foreign_energy += amount;
creature.local_energy -= amount;
}
}
// steal energy from frozen creatures
function creature_steal(creature, world) {
if (creature.frozen === 'frozen') { return; }
const pos = creature_forward(creature, world);
if (!world_occupied(world, pos)) {
// no target...
return;
} else {
const target_id = world.pos_creature[pos];
const target = world.creatures[target_id];
creature.foreign_energy += Math.floor(0.5 * target.local_energy);
target.local_energy = 0;
}
}
// hurl a creature over yourself
function creature_hurl(creature, world) {
if (creature.frozen === 'frozen') { return; }
const pos = creature_forward(creature, world);
if (!world_occupied(world, pos)) {
// no target...
return;
}
const dup = { ...creature };
creature_turn_right(dup);
creature_turn_right(dup);
const behind = creature_forward(dup, world);
if (world_occupied(world, behind)) {
// blocked
return;
}
const target_id = world.pos_creature[pos];
world_creature_move(world, target_id, behind);
}
// move
function creature_move(creature, world) {
if (creature.frozen === 'frozen') { return; }
const pos = creature_forward(creature, world);
world_creature_move(world, creature.id, pos);
}
// turn left
function creature_turn_left(creature) {
if (creature.frozen === 'frozen') { return; }
switch (creature.direction) {
case 'north':
creature.direction = 'west';
break;
case 'west':
creature.direction = 'south';
break;
case 'south':
creature.direction = 'east';
break;
case 'east':
creature.direction = 'north';
break;
default:
throw new Error(`bad direction: ${creature.direction}`);
}
}
// turn right
function creature_turn_right(creature) {
if (creature.frozen === 'frozen') { return; }
switch (creature.direction) {
case 'north':
creature.direction = 'east';
break;
case 'west':
creature.direction = 'north';
break;
case 'south':
creature.direction = 'west';
break;
case 'east':
creature.direction = 'south';
break;
default:
throw new Error(`bad direction: ${creature.direction}`);
}
}
function make_act(world) {
return (creature, value) => {
const choice = value & 0xf;
const rest = value >> 4;
switch (choice) {
case 0:
case 1:
case 2:
creature_reproduce(creature, world);
break;
case 3:
creature_hurl(creature, world);
break;
case 4:
creature_speak(creature, world, rest);
break;
case 5:
creature_freeze(creature, world);
break;
case 6:
creature_unfreeze(creature, world);
break;
case 7:
creature_set_pretend(creature, world, rest);
break;
case 8:
creature_give(creature, world, rest);
break;
case 9:
creature_steal(creature, world);
break;
case 10:
case 11:
creature_move(creature, world);
break;
case 12:
case 13:
creature_turn_left(creature);
break;
case 14:
case 15:
creature_turn_right(creature);
break;
}
};
}
// ================================================================
//
// setup
//
// ================================================================
let world = world_create(256, 256);
function world_insert_all(world, creatures) {
creatures.forEach(
x => {
let id = world_creature_insert(world, x, [ rand_uint(), rand_uint() ]);
while (id === null) {
id = world_creature_insert(world, x, [ rand_uint(), rand_uint() ]);
}
x.id = id;
}
);
}
function random_creature() {
const genome = [...Array(200)]
.reduce(genome_mut_append, { size: 128, pairs: [] });
return creature_create(genome);
}
let epoch = null;
function random_start() {
random_seed = Math.floor(0x100000000 * Math.random());
world = world_create(256, 256);
// have no saved epochs
const repl_genome = {
size: 16,
pairs: [
[0, 15], [1, 12], // turn left
[2, 15], [3, 0], // replicate
],
};
const travel_genome = {
size: 128,
pairs: [
[0, 15], [1, 10], // forward
[2, 15], [3, 10], // forward
[4, 15], [5, 10], // forward
[6, 15], [7, 10], // forward
[8, 15], [9, 10], // forward
[10, 15], [11, 10], // forward
[12, 15], [13, 10], // forward
[14, 15], [15, 0], // reproduce
[16, 15], [17, 12], // turn left
[18, 11], [19, 15], // jump to start
],
};
const replicators = [...Array(48)].map(() => creature_create(repl_genome));
const travellers = [...Array(16)].map(() => creature_create(travel_genome));
const randoms = [...Array(512)].map(random_creature);
const creatures = [ ...replicators, ...travellers, ...randoms ];
world_insert_all(world, creatures);
epoch = 0;
}
function load_start(input) {
console.log('start');
const file = input.files[0];
const reader = new FileReader();
reader.readAsText(file);
reader.onload = () => {
console.log(reader.result);
const result = JSON.parse(reader.result);
random_seed = result.random_seed;
epoch = result.epoch;
world = world_create(256, 256);
world_insert_all(world, result.creatures);
}
reader.onerror = () => {
alert(reader.error);
}
}
window.onload = () => {
const header = document.getElementById("status-header");
const stats = document.getElementById("stats-header");
const canvas = document.getElementById("render-canvas");
canvas.width = 3 * world.width;
canvas.height = 3 * world.height;
const ctx = canvas.getContext('2d');
ctx.scale(3, 3);
function world_render(world) {
ctx.fillStyle = '#ffffff';
ctx.fillRect(0, 0, world.width, world.height);
[...Object.values(world.creatures)].forEach(creature => {
const m = creature.mind.memory;
ctx.fillStyle = `rgb(${m[0] % 256}, ${m[1] % 256}, ${m[2] % 256})`;
const [x, y] = world.creature_pos[creature.id];
ctx.fillRect(x, y, 1, 1);
});
}
// from https://stackoverflow.com/questions/13405129/create-and-save-a-file-with-javascript
function download(data, filename, type) {
var file = new Blob([data], {type: type});
if (window.navigator.msSaveOrOpenBlob) // IE10+
window.navigator.msSaveOrOpenBlob(file, filename);
else { // Others
var a = document.createElement("a"),
url = URL.createObjectURL(file);
a.href = url;
a.download = filename;
document.body.appendChild(a);
a.click();
setTimeout(function() {
document.body.removeChild(a);
window.URL.revokeObjectURL(url);
}, 0);
}
}
function epoch_step() {
const choices = [
...[...Object.values(world.creatures)]
.sort((a, b) => a.foreign_energy - b.foreign_energy)
.filter(creature => creature.frozen !== 'frozen')
.toSpliced(500),
...[...Array(12)].map(random_creature),
];
world.creatures = {};
world.creature_pos = {};
world.pos_creature = {};
world.time = 0;
world_insert_all(world, choices);
epoch += 1;
download(
JSON.stringify({ epoch, random_seed, creatures: choices }),
`epoch${epoch}.json`,
'application/json',
);
step();
}
function step() {
header.innerText = `epoch ${epoch}, step ${world.time}`;
world_step(world);
world_render(world);
if (world.time % 4 === 0) {
const creatures = [...Object.values(world.creatures)];
stats.innerText = `frozen: ${0.1*Math.floor(creatures.filter(x => x.frozen === 'frozen').length / creatures.length * 1000)}%`;
}
if (world.time < 8192) {
setTimeout(step, 0);
} else {
epoch_step();
}
};
random_start();
step();
}

View file

@ -1,12 +0,0 @@
src/genome
==========
Genomes represent the neural network that underlies a creature. They are an array of tuples of the form `[ source, sink, weight ]`, where:
* `source` is an integer in the range `[0, N - num_outputs)`, representing the index of the source neuron
* `sink` is an integer in the range `[num_inputs, N)`, representing the index of the sink neuron
* `weight` is a floating-point value in the range `[-4, 4]`, representing the weight of the connection
`num_input` and `num_output` are fixed by the environment (as they are the input senses and output actions of the creature, respectively).
`N` is not fixed, but is instead determined by the maximum index present in the genome. As long as the maximum index is greater than or
equal to `num_input + num_output`, the genome is considered valid.

View file

@ -1,218 +0,0 @@
'use strict';
import { random_choice } from '../util.js';
import { network } from '../mind/topology.js';
// check if a given genome is valid and compute its size
export function validate_genome(genome) {
const { n_input, n_internal, n_output } = genome;
console.log(n_input + n_internal);
return genome.genes.reduce(
(acc, [source, sink, weight]) => acc && (
(source < n_input+n_internal) &&
(sink >= n_input)
),
true
);
}
// parse a genome into a useable neural net
export function parse_genome(genome) {
const { n_input, n_internal, n_output } = genome;
const n = genome.genes.reduce(
(acc, [source, sink, weight]) => acc.connect(source, sink, weight),
network(n_input, n_internal, n_output)
);
return n;
}
// --===== mutations =====--
function clamp(value, min, max) {
if (value > max) { return max; }
if (value < min) { return min; }
return value;
}
// adjust the source input of a gene
export function mut_gene_source(n_input, n_internal, n_output, gene, r) {
const [source, sink, weight] = gene;
const new_source = r < 0.5 ? source-1 : source+1;
return [
clamp(new_source, 0, n_input+n_internal-1),
sink,
weight,
];
}
// adjust the sink target of a gene
export function mut_gene_sink(n_input, n_internal, n_output, gene, r) {
const [source, sink, weight] = gene;
const new_sink = r < 0.5 ? sink-1 : sink+1;
return [
source,
clamp(new_sink, n_input+n_internal, n_input+n_internal+n_output-1),
weight,
];
}
// modify a gene's weight
// only adjusts the weight by performing a weighted average, so as to
// more gently modify the generated net
export function mut_gene_weight(weight_max, gene, r) {
const [source, sink, weight] = gene;
const rr = (2*r)-1;
const move = weight_max * rr;
const new_weight = (2*weight + move)/3;
return [
source,
sink,
clamp(new_weight, -weight_max, weight_max),
];
}
// expand the size of the neural net encoded by the genome
// relabels internal indices so that there is one extra internal neuron
export function mut_genome_expand(genome, r) {
const expand_index = Math.floor(genome.n_internal * r) + genome.n_input;
const new_genes = genome.genes.map(([source, sink, weight]) => [
source >= expand_index ? source+1 : source,
sink >= expand_index ? sink+1 : sink,
weight,
]);
return {
...genome,
n_internal: genome.n_internal+1,
genes: new_genes,
};
}
// contract the size of the neural net encoded by the genome
// relabels internal indices so that there is one less internal neuron
export function mut_genome_contract(genome, r) {
const { n_input, n_internal, n_output } = genome;
const contract_idx = Math.floor(n_internal * r) + n_input;
// decrement sources on the contract index too, to prevent invalid genomes
const new_source = (source) => source >= contract_idx ? source-1 : source;
// decrement sinks only after the contract index
const new_sink = (sink) => sink > contract_idx ? sink-1 : sink;
const new_genes = genome.genes.map(([source, sink, weight]) => [
new_source(source),
new_sink(sink),
weight,
]);
return {
...genome,
n_internal: n_internal-1,
genes: new_genes,
};
}
// append a newly generated gene to the end of the genome
export function mut_genome_insert(genome, weight_max, r1=Math.random(), r2=Math.random(), r3=Math.random()) {
const { n_input, n_internal, n_output } = genome;
const source = Math.floor((n_input + n_internal) * r1);
const sink = Math.floor((n_internal + n_output) * r2) + n_input;
const weight = weight_max * ((2*r3)-1);
return {
...genome,
genes: [...genome.genes, [source, sink, weight]],
};
}
// delete a gene from the genome
export function mut_genome_delete(genome, r) {
const del_idx = Math.floor(r * genome.genes.length);
const genes = genome.genes.filter((_, idx) => idx != del_idx);
return { ...genome, genes };
}
function mut_gene(
[n_input, n_internal, n_output, genome],
weight_max, r1, r2, r3
) {
const gene_idx = Math.floor(genome.length * r1);
const mod = random_choice(['source', 'sink', 'weight'], r2);
let new_gene;
if (mod == 'source') {
new_gene = mut_gene_source(
n_input, n_internal, n_output,
genome[gene_idx],
r3
);
} else if (mod == 'sink') {
new_gene = mut_gene_sink(
n_input, n_internal, n_output,
genome[gene_idx],
r3
);
} else {
new_gene = mut_gene_weight(
weight_max, genome[gene_idx], r3
);
}
const new_genome = genome.map((gene, idx) => {
if (idx == gene_idx) { return new_gene; }
return gene;
});
return [
n_input, n_internal, n_output, new_genome
];
}
export function mutate_genome(obj, weight_max) {
const mut = random_choice([
'gene', 'gene', 'gene',
'gene', 'gene', 'gene',
'gene', 'gene', 'gene',
'insert', 'delete',
'insert', 'delete',
'expand', 'contract',
], Math.random());
if (mut == 'gene') {
return mut_gene(
obj, weight_max,
Math.random(), Math.random(), Math.random()
);
} else if (mut == 'insert') {
return mut_genome_insert(
obj, weight_max,
Math.random(), Math.random(), Math.random()
);
} else if (mut == 'delete') {
return mut_genome_delete(obj, Math.random());
} else if (mut == 'expand') {
return mut_genome_expand(obj, Math.random());
} else if (mut == 'contract') {
return mut_genome_contract(obj, Math.random());
} else {
throw new Error(`bad mut value: ${mut}`);
}
}

View file

@ -1,282 +0,0 @@
'use strict';
// genome structure
// {
// genes: gene[]
// n_input, n_internal, n_output
// }
import {
validate_genome,
parse_genome,
mut_gene_source,
mut_gene_sink,
mut_gene_weight,
mut_genome_expand,
mut_genome_contract,
mut_genome_insert,
mut_genome_delete,
} from './genome';
test('genome validation', () => {
expect(validate_genome({
n_input: 0, n_internal: 1, n_output: 0,
genes: [[0, 0, 1.0]],
})).toBe(true);
expect(validate_genome({
n_input: 2, n_internal: 0, n_output: 1,
genes: [[0, 2, 1]],
})).toBe(true);
expect(validate_genome({
n_input: 2, n_internal: 0, n_output: 1,
genes: [[2, 0, 1]],
})).toBe(false);
expect(validate_genome({
n_input: 2, n_internal: 0, n_output: 1,
genes: [[2, 2, 1]],
})).toBe(false);
expect(validate_genome({
n_input: 2, n_internal: 1, n_output: 1,
genes: [[3, 2, 1]],
})).toBe(false);
});
test('parse a genome into a neural net', () => {
const n = parse_genome({
n_input: 1, n_internal: 1, n_output: 1,
genes: [
[0, 1, 1],
[1, 1, 1],
[1, 2, 1]
]
});
expect(n.input_count).toBe(1);
expect(n.output_count).toBe(1);
expect(n.compute([2], [-1])).toEqual([
[ Math.tanh( Math.tanh( 2-1 ) ) ],
[ Math.tanh( 2-1 ) ],
]);
});
test('mutate gene source', () => {
const n_input = 3;
const n_internal = 4;
const n_output = 5;
expect(mut_gene_source(
n_input, n_internal, n_output,
[0, 4, 0],
0.0
)).toEqual([0, 4, 0]);
expect(mut_gene_source(
n_input, n_internal, n_output,
[0, 4, 0],
1.0
)).toEqual([1, 4, 0]);
expect(mut_gene_source(
n_input, n_internal, n_output,
[6, 4, 0],
0.0
)).toEqual([5, 4, 0]);
expect(mut_gene_source(
n_input, n_internal, n_output,
[6, 4, 0],
1.0
)).toEqual([6, 4, 0]);
});
test('mutate gene sink', () => {
const n_input = 3;
const n_internal = 4;
const n_output = 5;
expect(mut_gene_sink(
n_input, n_internal, n_output,
[0, 7, 0],
0.0
)).toEqual([0, 7, 0]);
expect(mut_gene_sink(
n_input, n_internal, n_output,
[0, 7, 0],
1.0
)).toEqual([0, 8, 0]);
expect(mut_gene_sink(
n_input, n_internal, n_output,
[6, 11, 0],
0.0
)).toEqual([6, 10, 0]);
expect(mut_gene_sink(
n_input, n_internal, n_output,
[6, 11, 0],
1.0
)).toEqual([6, 11, 0]);
});
test('mutate gene weight', () => {
const weight_max = 4.0;
expect(mut_gene_weight(
weight_max, [0, 0, 1], 0.0
)).toEqual([0, 0, (2 - 4)/3]);
expect(mut_gene_weight(
weight_max, [0, 0, -4], 1.0
)).toEqual([0, 0, (-8 + 4)/3]);
expect(mut_gene_weight(
weight_max, [0, 0, 3], 0.5
)).toEqual([0, 0, (6+0)/3]);
});
test('expand genome', () => {
const n_input = 1;
const n_internal = 3;
const n_output = 1;
const genome = {
n_input, n_internal, n_output,
genes: [
[0, 1, 0],
[1, 2, 0],
[2, 3, 0],
[3, 4, 0],
],
};
expect(mut_genome_expand(genome, 0.0)).toEqual({
n_input, n_internal: n_internal+1, n_output,
genes: [
[0, 2, 0],
[2, 3, 0],
[3, 4, 0],
[4, 5, 0],
],
});
expect(mut_genome_expand(genome, 0.5)).toEqual({
n_input, n_internal: n_internal+1, n_output,
genes: [
[0, 1, 0],
[1, 3, 0],
[3, 4, 0],
[4, 5, 0],
],
});
expect(mut_genome_expand(genome, 0.99)).toEqual({
n_input, n_internal: n_internal+1, n_output,
genes: [
[0, 1, 0],
[1, 2, 0],
[2, 4, 0],
[4, 5, 0],
],
});
});
test('contract genome', () => {
const n_input = 1;
const n_internal = 3;
const n_output = 1;
const genome = {
n_input, n_internal, n_output,
genes: [
[0, 1, 0],
[1, 2, 1],
[2, 3, 2],
[3, 4, 3],
],
};
expect(mut_genome_contract(genome, 0.0)).toEqual({
n_input, n_internal: n_internal-1, n_output,
genes: [
[0, 1, 0],
[0, 1, 1],
[1, 2, 2],
[2, 3, 3],
],
});
expect(mut_genome_contract(genome, 0.5)).toEqual({
n_input, n_internal: n_internal-1, n_output,
genes: [
[0, 1, 0],
[1, 2, 1],
[1, 2, 2],
[2, 3, 3],
],
});
expect(mut_genome_contract(genome, 0.99)).toEqual({
n_input, n_internal: n_internal-1, n_output,
genes: [
[0, 1, 0],
[1, 2, 1],
[2, 3, 2],
[2, 3, 3],
],
});
});
test('insert new genes', () => {
const n_input = 1;
const n_internal = 2;
const n_output = 1;
const weight_max = 4;
expect(mut_genome_insert({
n_input, n_internal, n_output,
genes: []
}, weight_max, 0, 0.5, 0)).toEqual({
n_input, n_internal, n_output,
genes: [[0, 2, -4]]
});
expect(mut_genome_insert({
n_input, n_internal, n_output,
genes: [[0, 2, -4]]
}, weight_max, 0.99, 0, 1)).toEqual({
n_input, n_internal, n_output,
genes: [[0, 2, -4], [2, 1, 4]]
});
});
test('remove genes', () => {
const n_input = 0;
const n_output = 0;
const n_internal = 3;
const genome = {
n_input, n_internal, n_output,
genes: [[0, 1, 0], [1, 2, 0]],
};
expect(mut_genome_delete(genome, 0.0)).toEqual({
n_input, n_internal, n_output,
genes: [[1, 2, 0]],
});
expect(mut_genome_delete(genome, 0.99)).toEqual({
n_input, n_internal, n_output,
genes: [[0, 1, 0]],
});
});

View file

@ -1,37 +0,0 @@
import {
get_size,
mut_genome_insert, mutate_genome,
} from './genome.js';
const recurse = (f, x0, n) => {
if (n == 0) {
return x0;
} else {
return f(recurse(f, x0, n-1));
}
};
const n_input = 5;
const n_output = 5;
const [_1, _2, _3, genome] = recurse(
s => mut_genome_insert(
s, 4,
Math.random(), Math.random(), Math.random()
),
[n_input, 10, n_output, []],
20);
const n_internal = get_size(n_input, n_output, genome) - n_input - n_output;
console.log([n_input, n_internal, n_output, genome]);
const mutation = recurse(
s => mutate_genome(s, 4),
[n_input, n_internal, n_output, genome],
40
);
console.log(mutation);

View file

@ -1,14 +0,0 @@
<!doctype html>
<html>
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>nerine</title>
<style>
</style>
<script type="module" src="ui/index.js"></script>
</head>
<body>
<canvas id="canvas" width="95vw" height="95vw">
</body>
</html>

View file

@ -1,37 +0,0 @@
mind
====
This module is used to create arbitrary stateful neural networks.
The only export is the following function:
```
network(input_count, internal_count, output_count, weight_max = 4 : number)
```
This function returns an object that represents a neural network with `input_count` input neurons,
`internal_count` internal (and stateful) neurons, and `output_count` output neurons. `max_weight` determines
the maximum absolute value allowed for connection weights.
A network object has two methods:
```
connect(source, sink, weight : number)
```
This method returns a new network object that is an exact copy of the original, except with a new
connection between the neuron indexed by `source` and `sink`, with weight `weight`. Neuron indices
are zero-indexed, and span all neurons, first the inputs, then the internal neurons, and finally the outputs.
An error will be thrown if `source` is in the range of output neurons or if `sink` is in the range of input
neurons.
```
compute(inputs, state : array[number])
```
This method returns a tuple `[output, newState]`, where `output` is an array of `output_count` values
corresponding to the output neuron's computed values, and `newState` is the new state of the internal neurons.
`input` must be an array of numbers with length equal to `input_count`, and `state` must be an array of numbers
with length equal to `internal_count` or an error will be raised.

View file

@ -1,183 +0,0 @@
'use strict';
import { create } from '../util.js';
const DEFAULT_WEIGHT_MAX = 4;
// prototype for network objects
const network_proto = {
connect: function(source, sink, weight) {
return network_connect(this, source, sink, weight);
},
compute: function(inputs, state) {
return network_compute(this, inputs, state);
},
};
// create a new network
export function network(input_count, internal_count, output_count, weight_max = 4) {
const count = input_count + internal_count + output_count;
const n = create({
input_count,
output_count,
adjacency: new Array(count).fill([]),
weight: [],
}, network_proto);
return n;
}
// check index is an input
function is_input(n, index) {
return index < n.input_count;
}
// check if index is an output
function is_output(n, index) {
return index >= (n.adjacency.length - n.output_count);
}
// check if index is a hidden neuron
function is_hidden(n, index) {
return (!is_input(n, index)) && (!is_output(n, index));
}
// returns a new network with an edge between the given nodes
// with the given weight
export function network_connect(n, source, sink, weight) {
if (is_input(n, sink)) {
// inputs cannot be sinks
throw new Error(`attempt to use input as sink (${source} -> ${sink})`);
}
if (is_output(n, source)) {
// outputs cannot be sources
throw new Error(`attempt to use output as source (${source} -> ${sink})`);
}
return create({
...n,
adjacency: n.adjacency.map((row, i) => {
if (i === source && i === sink) {
// self-loop
return [...row, 2];
} else if (i === source) {
return [...row, 1];
} else if (i === sink) {
return [...row, -1];
} else {
return [...row, 0];
}
}),
weight: [...n.weight, weight],
}, network_proto);
}
// gets the indices of the edges incident on the given adjacency list
function incident_edges(n, adj) {
const incident = adj
.map((edge, index) => (edge < 0) || (edge === 2) ? index : null)
.filter(index => index !== null);
return incident;
}
// get the indices of the ends of an edge
// in the case of self-loops, both values are the same
function edge_ends(n, edge) {
const ends = n.adjacency
.map((adj, index) => adj[edge] !== 0 ? index : null)
.filter(index => index != null);
ends.sort((a, b) => n.adjacency[a][edge] < n.adjacency[b][edge] ? -1 : 1);
if (ends.length === 1) {
return { source: ends[0], sink: ends[0] };
} else if (ends.length === 2) {
return { source: ends[1], sink: ends[0] };
} else {
throw new Error("something bad happened with the ends");
}
}
// recursively get the value of a node from the input nodes,
// optionally caching the computed values
function get_value(n, index, input, prev, cache) {
// check if value is cached
if (cache !== undefined && cache[index]) {
return cache[index];
}
// check if value is input
if (is_input(n, index)) {
return input[index];
}
const adj = n.adjacency[index]; // get adjacency list
const incident = incident_edges(n, adj); // get incident edges
const weight = incident.map(x => n.weight[x]); // edge weights
const sources = incident // get ancestor nodes
.map(x => edge_ends(n, x).source);
// get the value of each ancestor
const values = sources
.map(x => x === index // if the ancestor is this node
? prev[x - n.input_count] // then the value is the previous value
: get_value(n, x, input, prev, cache)); // else recurse
const sum = values // compute the weighted sum of the values
.reduce((acc, x, i) => acc + (weight[i] * x), 0);
if (sum !== sum) { // NaN test
console.log(n);
console.log(sources);
console.log(input);
throw new Error(`failed to get output for index ${index}`);
}
// compute result
const value = Math.tanh(sum);
// !!! impure caching !!!
// cache result
if (cache !== undefined) {
cache[index] = value;
}
return value;
}
// compute a network's output and new hidden state
// given the input and previous hidden state
export function network_compute(n, input, state) {
// validate input
if (input.length !== n.input_count) {
throw new Error("incorrect number of input elements");
}
// validate state
const hidden_count = n.adjacency.length - n.input_count - n.output_count;
if (state.length !== hidden_count) {
throw new Error("incorrect number of state elements");
}
// !!! impure caching !!!
const value_cache = {};
const result = Object.freeze(n.adjacency
.map((x, i) => is_output(n, i) ? i : null) // output index or null
.filter(i => i !== null) // remove nulls
.map(x => get_value(n, x, input, state, value_cache)) // map to computed value
);
const newstate = Object.freeze(n.adjacency
.map((x, i) => is_hidden(n, i) ? i : null) // hidden index or null
.filter(i => i !== null) // remove nulls
.map(x => get_value(n, x, input, state, value_cache)) // map to computed value (using cache)
);
return Object.freeze([result, newstate]);
}

View file

@ -1,235 +0,0 @@
'use strict';
import { network } from './topology';
test('basic network functionality', () => {
const n = network(0, 5, 0);
expect(n).toEqual({
input_count: 0,
output_count: 0,
adjacency: [ [], [], [], [], [] ],
weight: [],
});
expect(() => n.adjacency = []).toThrow();
expect(() => n.weight = []).toThrow();
const nn = n.connect(0, 1, -2);
expect(nn).toEqual({
input_count: 0,
output_count: 0,
adjacency: [
[ 1 ],
[ -1 ],
[ 0 ],
[ 0 ],
[ 0 ]
],
weight: [ -2 ],
});
expect(() => nn.adjacency = []).toThrow();
expect(() => nn.weight = []).toThrow();
const nnn = nn.connect(2, 4, 3);
expect(nnn).toEqual({
input_count: 0,
output_count: 0,
adjacency: [
[ 1, 0 ],
[ -1, 0 ],
[ 0, 1 ],
[ 0, 0 ],
[ 0, -1 ]
],
weight: [ -2, 3 ],
});
expect(() => nnn.adjacency = []).toThrow();
expect(() => nnn.weight = []).toThrow();
});
test(
'networks are restricted from sinking to inputs or sourcing from outputs',
() => {
const n = network(2, 2, 2);
expect(n.connect(1,2,0)).toEqual({
input_count: 2,
output_count: 2,
adjacency: [
[ 0 ],
[ 1 ],
[ -1 ],
[ 0 ],
[ 0 ],
[ 0 ],
],
weight: [ 0 ],
});
expect(() => n.connect(2, 1, 0)).toThrow();
expect(n.connect(3, 4, 2)).toEqual({
input_count: 2,
output_count: 2,
adjacency: [
[ 0 ],
[ 0 ],
[ 0 ],
[ 1 ],
[ -1 ],
[ 0 ],
],
weight: [ 2 ],
});
expect(() => n.connect(4, 3, 2)).toThrow();
});
test('self-connections work correctly', () => {
const n = network(0, 1, 0).connect(0, 0, 2.0);
expect(n).toEqual({
input_count: 0,
output_count: 0,
adjacency: [
[ 2 ],
],
weight: [ 2 ],
});
});
test('network computations', () => {
const n = network(1, 0, 1).connect(0, 1, 2.0);
const input = [ -0.5 ];
const state = [];
const result = n.compute(input, state);
expect(result).toEqual([
[ Math.tanh(-0.5 * 2.0) ],
[],
]);
expect(input).toEqual([ -0.5 ]);
expect(state).toEqual([]);
expect(() => result[0] = 'hi').toThrow();
expect(() => result[0].push('hi')).toThrow();
expect(() => result[1] = 'hi').toThrow();
expect(() => result[1].push('hi')).toThrow();
});
test('multiple input network', () => {
const n = network(4, 0, 1)
.connect(0, 4, -1.0)
.connect(1, 4, -2.0)
.connect(2, 4, 1.0)
.connect(3, 4, 2.0)
expect(n.compute([1, 2, 3, 5], [])).toEqual([
[ Math.tanh(
(-1.0 * 1) +
(-2.0 * 2) +
(1.0 * 3) +
(2.0 * 5))],
[],
]);
});
test('multiple outputs', () => {
const n = network(4, 0, 2)
.connect(0, 4, -1)
.connect(1, 4, 1)
.connect(2, 5, -1)
.connect(3, 5, 1);
expect(n.compute([1,2,3,5], [])).toEqual([
[ Math.tanh(2-1), Math.tanh(5-3) ],
[],
]);
});
test('hidden neurons', () => {
const n = network(4, 2, 1)
.connect(0, 4, -1)
.connect(1, 4, 1)
.connect(2, 5, -1)
.connect(3, 5, 1)
.connect(4, 6, -1)
.connect(5, 6, 1);
expect(n.compute([1,2,3,5], [ 0, 0 ])).toEqual([
[ Math.tanh( Math.tanh(5-3) - Math.tanh(2-1) ) ],
[ Math.tanh(2-1), Math.tanh(5-3) ],
]);
});
test('arbitrary hidden neurons', () => {
const n = network(1, 2, 1)
.connect(0, 1, 1)
.connect(1, 2, -1)
.connect(2, 3, 2)
const [output, state] = n.compute([1], [0, 0]);
expect(output).toEqual([
Math.tanh(
2*Math.tanh(
-1*Math.tanh(1)
)
)
]);
expect(state).toEqual([
Math.tanh(1),
Math.tanh( -Math.tanh(1) ),
]);
});
test('memory', () => {
const n = network(0, 1, 1).connect(0, 0, -0.5).connect(0, 1, 2);
expect(n.compute([], [1])).toEqual([
[ Math.tanh( 2 * Math.tanh( -0.5 * 1 ) ) ],
[ Math.tanh( -0.5 * 1) ],
]);
});
test('memory and input', () => {
const n = network(1, 1, 1)
.connect(0, 1, 1)
.connect(1, 1, 1)
.connect(1, 2, 1);
expect(n.compute([2], [-1])).toEqual([
[ Math.tanh( Math.tanh( 2-1 ) ) ],
[ Math.tanh( 2-1 ) ],
]);
});
test('input and state must be the correct size', () => {
const n = network(2, 1, 1)
.connect(0, 2, 1)
.connect(1, 2, 1)
.connect(2, 3, 1);
// wrong input size
expect(() => n.compute([], [4])).toThrow();
expect(() => n.compute([1], [4])).toThrow();
expect(() => n.compute([1, 1, 1], [4])).toThrow();
// wrong state size
expect(() => n.compute([1, 1], [])).toThrow();
expect(() => n.compute([1, 1], [4, 4])).toThrow();
// prove correct sizes work
n.compute([1, 1], [4]);
});

8
src/nerine.h Normal file
View file

@ -0,0 +1,8 @@
#ifndef NERINE_H
#define NERINE_H
int add(int a, int b) {
return -1;
}
#endif

Some files were not shown because too many files have changed in this diff Show more