1 #ifndef _CS_H 2 #define _CS_H 3 #include <stdlib.h> 4 #include <limits.h> 5 #include <math.h> 6 #include <stdio.h> 7 #include <stddef.h> 8 #ifdef MATLAB_MEX_FILE 9 #include "mex.h" 10 #endif 11 #define CS_VER 3 /* CSparse Version */ 12 #define CS_SUBVER 2 13 #define CS_SUBSUB 0 14 #define CS_DATE "Sept 12, 2017" /* CSparse release date */ 15 #define CS_COPYRIGHT "Copyright (c) Timothy A. Davis, 2006-2016" 16 17 #ifdef MATLAB_MEX_FILE 18 #undef csi 19 #define csi mwSignedIndex 20 #endif 21 #ifndef csi 22 #define csi ptrdiff_t 23 #endif 24 25 /* --- primary CSparse routines and data structures ------------------------- */ 26 typedef struct cs_sparse /* matrix in compressed-column or triplet form */ 27 { 28 csi nzmax ; /* maximum number of entries */ 29 csi m ; /* number of rows */ 30 csi n ; /* number of columns */ 31 csi *p ; /* column pointers (size n+1) or col indices (size nzmax) */ 32 csi *i ; /* row indices, size nzmax */ 33 double *x ; /* numerical values, size nzmax */ 34 csi nz ; /* # of entries in triplet matrix, -1 for compressed-col */ 35 } cs ; 36 37 cs *cs_add (const cs *A, const cs *B, double alpha, double beta) ; 38 csi cs_cholsol (csi order, const cs *A, double *b) ; 39 cs *cs_compress (const cs *T) ; 40 csi cs_dupl (cs *A) ; 41 csi cs_entry (cs *T, csi i, csi j, double x) ; 42 csi cs_gaxpy (const cs *A, const double *x, double *y) ; 43 cs *cs_load (FILE *f) ; 44 csi cs_lusol (csi order, const cs *A, double *b, double tol) ; 45 cs *cs_multiply (const cs *A, const cs *B) ; 46 double cs_norm (const cs *A) ; 47 csi cs_print (const cs *A, csi brief) ; 48 csi cs_qrsol (csi order, const cs *A, double *b) ; 49 cs *cs_transpose (const cs *A, csi values) ; 50 /* utilities */ 51 void *cs_calloc (csi n, size_t size) ; 52 void *cs_free (void *p) ; 53 void *cs_realloc (void *p, csi n, size_t size, csi *ok) ; 54 cs *cs_spalloc (csi m, csi n, csi nzmax, csi values, csi triplet) ; 55 cs *cs_spfree (cs *A) ; 56 csi cs_sprealloc (cs *A, csi nzmax) ; 57 void *cs_malloc (csi n, size_t size) ; 58 59 /* --- secondary CSparse routines and data structures ----------------------- */ 60 typedef struct cs_symbolic /* symbolic Cholesky, LU, or QR analysis */ 61 { 62 csi *pinv ; /* inverse row perm. for QR, fill red. perm for Chol */ 63 csi *q ; /* fill-reducing column permutation for LU and QR */ 64 csi *parent ; /* elimination tree for Cholesky and QR */ 65 csi *cp ; /* column pointers for Cholesky, row counts for QR */ 66 csi *leftmost ; /* leftmost[i] = min(find(A(i,:))), for QR */ 67 csi m2 ; /* # of rows for QR, after adding fictitious rows */ 68 double lnz ; /* # entries in L for LU or Cholesky; in V for QR */ 69 double unz ; /* # entries in U for LU; in R for QR */ 70 } css ; 71 72 typedef struct cs_numeric /* numeric Cholesky, LU, or QR factorization */ 73 { 74 cs *L ; /* L for LU and Cholesky, V for QR */ 75 cs *U ; /* U for LU, R for QR, not used for Cholesky */ 76 csi *pinv ; /* partial pivoting for LU */ 77 double *B ; /* beta [0..n-1] for QR */ 78 } csn ; 79 80 typedef struct cs_dmperm_results /* cs_dmperm or cs_scc output */ 81 { 82 csi *p ; /* size m, row permutation */ 83 csi *q ; /* size n, column permutation */ 84 csi *r ; /* size nb+1, block k is rows r[k] to r[k+1]-1 in A(p,q) */ 85 csi *s ; /* size nb+1, block k is cols s[k] to s[k+1]-1 in A(p,q) */ 86 csi nb ; /* # of blocks in fine dmperm decomposition */ 87 csi rr [5] ; /* coarse row decomposition */ 88 csi cc [5] ; /* coarse column decomposition */ 89 } csd ; 90 91 csi *cs_amd (csi order, const cs *A) ; 92 csn *cs_chol (const cs *A, const css *S) ; 93 csd *cs_dmperm (const cs *A, csi seed) ; 94 csi cs_droptol (cs *A, double tol) ; 95 csi cs_dropzeros (cs *A) ; 96 csi cs_happly (const cs *V, csi i, double beta, double *x) ; 97 csi cs_ipvec (const csi *p, const double *b, double *x, csi n) ; 98 csi cs_lsolve (const cs *L, double *x) ; 99 csi cs_ltsolve (const cs *L, double *x) ; 100 csn *cs_lu (const cs *A, const css *S, double tol) ; 101 cs *cs_permute (const cs *A, const csi *pinv, const csi *q, csi values) ; 102 csi *cs_pinv (const csi *p, csi n) ; 103 csi cs_pvec (const csi *p, const double *b, double *x, csi n) ; 104 csn *cs_qr (const cs *A, const css *S) ; 105 css *cs_schol (csi order, const cs *A) ; 106 css *cs_sqr (csi order, const cs *A, csi qr) ; 107 cs *cs_symperm (const cs *A, const csi *pinv, csi values) ; 108 csi cs_updown (cs *L, csi sigma, const cs *C, const csi *parent) ; 109 csi cs_usolve (const cs *U, double *x) ; 110 csi cs_utsolve (const cs *U, double *x) ; 111 /* utilities */ 112 css *cs_sfree (css *S) ; 113 csn *cs_nfree (csn *N) ; 114 csd *cs_dfree (csd *D) ; 115 116 /* --- tertiary CSparse routines -------------------------------------------- */ 117 csi *cs_counts (const cs *A, const csi *parent, const csi *post, csi ata) ; 118 double cs_cumsum (csi *p, csi *c, csi n) ; 119 csi cs_dfs (csi j, cs *G, csi top, csi *xi, csi *pstack, const csi *pinv) ; 120 csi cs_ereach (const cs *A, csi k, const csi *parent, csi *s, csi *w) ; 121 csi *cs_etree (const cs *A, csi ata) ; 122 csi cs_fkeep (cs *A, csi (*fkeep) (csi, csi, double, void *), void *other) ; 123 double cs_house (double *x, double *beta, csi n) ; 124 csi cs_leaf (csi i, csi j, const csi *first, csi *maxfirst, csi *prevleaf, 125 csi *ancestor, csi *jleaf) ; 126 csi *cs_maxtrans (const cs *A, csi seed) ; 127 csi *cs_post (const csi *parent, csi n) ; 128 csi *cs_randperm (csi n, csi seed) ; 129 csi cs_reach (cs *G, const cs *B, csi k, csi *xi, const csi *pinv) ; 130 csi cs_scatter (const cs *A, csi j, double beta, csi *w, double *x, csi mark, 131 cs *C, csi nz) ; 132 csd *cs_scc (cs *A) ; 133 csi cs_spsolve (cs *G, const cs *B, csi k, csi *xi, double *x, 134 const csi *pinv, csi lo) ; 135 csi cs_tdfs (csi j, csi k, csi *head, const csi *next, csi *post, 136 csi *stack) ; 137 /* utilities */ 138 csd *cs_dalloc (csi m, csi n) ; 139 csd *cs_ddone (csd *D, cs *C, void *w, csi ok) ; 140 cs *cs_done (cs *C, void *w, void *x, csi ok) ; 141 csi *cs_idone (csi *p, cs *C, void *w, csi ok) ; 142 csn *cs_ndone (csn *N, cs *C, void *w, void *x, csi ok) ; 143 144 #define CS_MAX(a,b) (((a) > (b)) ? (a) : (b)) 145 #define CS_MIN(a,b) (((a) < (b)) ? (a) : (b)) 146 #define CS_FLIP(i) (-(i)-2) 147 #define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i)) 148 #define CS_MARKED(w,j) (w [j] < 0) 149 #define CS_MARK(w,j) { w [j] = CS_FLIP (w [j]) ; } 150 #define CS_CSC(A) (A && (A->nz == -1)) 151 #define CS_TRIPLET(A) (A && (A->nz >= 0)) 152 #endif 153