1 /* ========================================================================== */ 2 /* CXSparse/Include/cs.h file */ 3 /* ========================================================================== */ 4 5 /* This is the CXSparse/Include/cs.h file. It has the same name (cs.h) as 6 the CSparse/Include/cs.h file. The 'make install' for SuiteSparse installs 7 CXSparse, and this file, instead of CSparse. The two packages have the same 8 cs.h include filename, because CXSparse is a superset of CSparse. Any user 9 program that uses CSparse can rely on CXSparse instead, with no change to the 10 user code. The #include "cs.h" line will work for both versions, in user 11 code, and the function names and user-visible typedefs from CSparse all 12 appear in CXSparse. For experimenting and changing the package itself, I 13 recommend using CSparse since it's simpler and easier to modify. For 14 using the package in production codes, I recommend CXSparse since it has 15 more features (support for complex matrices, and both int and long 16 versions). 17 */ 18 19 /* ========================================================================== */ 20 21 #ifndef _CXS_H 22 #define _CXS_H 23 #include <stdlib.h> 24 #include <limits.h> 25 #include <math.h> 26 #include <stdio.h> 27 #ifdef MATLAB_MEX_FILE 28 #include "mex.h" 29 #endif 30 31 32 #ifdef __cplusplus 33 #ifndef NCOMPLEX 34 #include <complex> 35 typedef std::complex<double> cs_complex_t ; 36 #endif 37 extern "C" { 38 #else 39 #ifndef NCOMPLEX 40 #include <complex.h> 41 #define cs_complex_t double _Complex 42 #endif 43 #endif 44 45 #define CS_VER 3 /* CXSparse Version */ 46 #define CS_SUBVER 1 47 #define CS_SUBSUB 9 48 #define CS_DATE "May 4, 2016" /* CXSparse release date */ 49 #define CS_COPYRIGHT "Copyright (c) Timothy A. Davis, 2006-2016" 50 #define CXSPARSE 51 52 #include "SuiteSparse_config.h" 53 #define cs_long_t SuiteSparse_long 54 #define cs_long_t_id SuiteSparse_long_id 55 #define cs_long_t_max SuiteSparse_long_max 56 57 /* -------------------------------------------------------------------------- */ 58 /* double/int version of CXSparse */ 59 /* -------------------------------------------------------------------------- */ 60 61 /* --- primary CSparse routines and data structures ------------------------- */ 62 63 typedef struct cs_di_sparse /* matrix in compressed-column or triplet form */ 64 { 65 int nzmax ; /* maximum number of entries */ 66 int m ; /* number of rows */ 67 int n ; /* number of columns */ 68 int *p ; /* column pointers (size n+1) or col indices (size nzmax) */ 69 int *i ; /* row indices, size nzmax */ 70 double *x ; /* numerical values, size nzmax */ 71 int nz ; /* # of entries in triplet matrix, -1 for compressed-col */ 72 } cs_di ; 73 74 cs_di *cs_di_add (const cs_di *A, const cs_di *B, double alpha, double beta) ; 75 int cs_di_cholsol (int order, const cs_di *A, double *b) ; 76 int cs_di_dupl (cs_di *A) ; 77 int cs_di_entry (cs_di *T, int i, int j, double x) ; 78 int cs_di_lusol (int order, const cs_di *A, double *b, double tol) ; 79 int cs_di_gaxpy (const cs_di *A, const double *x, double *y) ; 80 cs_di *cs_di_multiply (const cs_di *A, const cs_di *B) ; 81 int cs_di_qrsol (int order, const cs_di *A, double *b) ; 82 cs_di *cs_di_transpose (const cs_di *A, int values) ; 83 cs_di *cs_di_compress (const cs_di *T) ; 84 double cs_di_norm (const cs_di *A) ; 85 int cs_di_print (const cs_di *A, int brief) ; 86 cs_di *cs_di_load (FILE *f) ; 87 88 /* utilities */ 89 void *cs_di_calloc (int n, size_t size) ; 90 void *cs_di_free (void *p) ; 91 void *cs_di_realloc (void *p, int n, size_t size, int *ok) ; 92 cs_di *cs_di_spalloc (int m, int n, int nzmax, int values, int t) ; 93 cs_di *cs_di_spfree (cs_di *A) ; 94 int cs_di_sprealloc (cs_di *A, int nzmax) ; 95 void *cs_di_malloc (int n, size_t size) ; 96 97 /* --- secondary CSparse routines and data structures ----------------------- */ 98 99 typedef struct cs_di_symbolic /* symbolic Cholesky, LU, or QR analysis */ 100 { 101 int *pinv ; /* inverse row perm. for QR, fill red. perm for Chol */ 102 int *q ; /* fill-reducing column permutation for LU and QR */ 103 int *parent ; /* elimination tree for Cholesky and QR */ 104 int *cp ; /* column pointers for Cholesky, row counts for QR */ 105 int *leftmost ; /* leftmost[i] = min(find(A(i,:))), for QR */ 106 int m2 ; /* # of rows for QR, after adding fictitious rows */ 107 double lnz ; /* # entries in L for LU or Cholesky; in V for QR */ 108 double unz ; /* # entries in U for LU; in R for QR */ 109 } cs_dis ; 110 111 typedef struct cs_di_numeric /* numeric Cholesky, LU, or QR factorization */ 112 { 113 cs_di *L ; /* L for LU and Cholesky, V for QR */ 114 cs_di *U ; /* U for LU, r for QR, not used for Cholesky */ 115 int *pinv ; /* partial pivoting for LU */ 116 double *B ; /* beta [0..n-1] for QR */ 117 } cs_din ; 118 119 typedef struct cs_di_dmperm_results /* cs_di_dmperm or cs_di_scc output */ 120 { 121 int *p ; /* size m, row permutation */ 122 int *q ; /* size n, column permutation */ 123 int *r ; /* size nb+1, block k is rows r[k] to r[k+1]-1 in A(p,q) */ 124 int *s ; /* size nb+1, block k is cols s[k] to s[k+1]-1 in A(p,q) */ 125 int nb ; /* # of blocks in fine dmperm decomposition */ 126 int rr [5] ; /* coarse row decomposition */ 127 int cc [5] ; /* coarse column decomposition */ 128 } cs_did ; 129 130 int *cs_di_amd (int order, const cs_di *A) ; 131 cs_din *cs_di_chol (const cs_di *A, const cs_dis *S) ; 132 cs_did *cs_di_dmperm (const cs_di *A, int seed) ; 133 int cs_di_droptol (cs_di *A, double tol) ; 134 int cs_di_dropzeros (cs_di *A) ; 135 int cs_di_happly (const cs_di *V, int i, double beta, double *x) ; 136 int cs_di_ipvec (const int *p, const double *b, double *x, int n) ; 137 int cs_di_lsolve (const cs_di *L, double *x) ; 138 int cs_di_ltsolve (const cs_di *L, double *x) ; 139 cs_din *cs_di_lu (const cs_di *A, const cs_dis *S, double tol) ; 140 cs_di *cs_di_permute (const cs_di *A, const int *pinv, const int *q, 141 int values) ; 142 int *cs_di_pinv (const int *p, int n) ; 143 int cs_di_pvec (const int *p, const double *b, double *x, int n) ; 144 cs_din *cs_di_qr (const cs_di *A, const cs_dis *S) ; 145 cs_dis *cs_di_schol (int order, const cs_di *A) ; 146 cs_dis *cs_di_sqr (int order, const cs_di *A, int qr) ; 147 cs_di *cs_di_symperm (const cs_di *A, const int *pinv, int values) ; 148 int cs_di_usolve (const cs_di *U, double *x) ; 149 int cs_di_utsolve (const cs_di *U, double *x) ; 150 int cs_di_updown (cs_di *L, int sigma, const cs_di *C, const int *parent) ; 151 152 /* utilities */ 153 cs_dis *cs_di_sfree (cs_dis *S) ; 154 cs_din *cs_di_nfree (cs_din *N) ; 155 cs_did *cs_di_dfree (cs_did *D) ; 156 157 /* --- tertiary CSparse routines -------------------------------------------- */ 158 159 int *cs_di_counts (const cs_di *A, const int *parent, const int *post, 160 int ata) ; 161 double cs_di_cumsum (int *p, int *c, int n) ; 162 int cs_di_dfs (int j, cs_di *G, int top, int *xi, int *pstack, 163 const int *pinv) ; 164 int *cs_di_etree (const cs_di *A, int ata) ; 165 int cs_di_fkeep (cs_di *A, int (*fkeep) (int, int, double, void *), 166 void *other) ; 167 double cs_di_house (double *x, double *beta, int n) ; 168 int *cs_di_maxtrans (const cs_di *A, int seed) ; 169 int *cs_di_post (const int *parent, int n) ; 170 cs_did *cs_di_scc (cs_di *A) ; 171 int cs_di_scatter (const cs_di *A, int j, double beta, int *w, double *x, 172 int mark, cs_di *C, int nz) ; 173 int cs_di_tdfs (int j, int k, int *head, const int *next, int *post, 174 int *stack) ; 175 int cs_di_leaf (int i, int j, const int *first, int *maxfirst, int *prevleaf, 176 int *ancestor, int *jleaf) ; 177 int cs_di_reach (cs_di *G, const cs_di *B, int k, int *xi, const int *pinv) ; 178 int cs_di_spsolve (cs_di *L, const cs_di *B, int k, int *xi, double *x, 179 const int *pinv, int lo) ; 180 int cs_di_ereach (const cs_di *A, int k, const int *parent, int *s, int *w) ; 181 int *cs_di_randperm (int n, int seed) ; 182 183 /* utilities */ 184 cs_did *cs_di_dalloc (int m, int n) ; 185 cs_di *cs_di_done (cs_di *C, void *w, void *x, int ok) ; 186 int *cs_di_idone (int *p, cs_di *C, void *w, int ok) ; 187 cs_din *cs_di_ndone (cs_din *N, cs_di *C, void *w, void *x, int ok) ; 188 cs_did *cs_di_ddone (cs_did *D, cs_di *C, void *w, int ok) ; 189 190 191 /* -------------------------------------------------------------------------- */ 192 /* double/cs_long_t version of CXSparse */ 193 /* -------------------------------------------------------------------------- */ 194 195 /* --- primary CSparse routines and data structures ------------------------- */ 196 197 typedef struct cs_dl_sparse /* matrix in compressed-column or triplet form */ 198 { 199 cs_long_t nzmax ; /* maximum number of entries */ 200 cs_long_t m ; /* number of rows */ 201 cs_long_t n ; /* number of columns */ 202 cs_long_t *p ; /* column pointers (size n+1) or col indlces (size nzmax) */ 203 cs_long_t *i ; /* row indices, size nzmax */ 204 double *x ; /* numerical values, size nzmax */ 205 cs_long_t nz ; /* # of entries in triplet matrix, -1 for compressed-col */ 206 } cs_dl ; 207 208 cs_dl *cs_dl_add (const cs_dl *A, const cs_dl *B, double alpha, double beta) ; 209 cs_long_t cs_dl_cholsol (cs_long_t order, const cs_dl *A, double *b) ; 210 cs_long_t cs_dl_dupl (cs_dl *A) ; 211 cs_long_t cs_dl_entry (cs_dl *T, cs_long_t i, cs_long_t j, double x) ; 212 cs_long_t cs_dl_lusol (cs_long_t order, const cs_dl *A, double *b, double tol) ; 213 cs_long_t cs_dl_gaxpy (const cs_dl *A, const double *x, double *y) ; 214 cs_dl *cs_dl_multiply (const cs_dl *A, const cs_dl *B) ; 215 cs_long_t cs_dl_qrsol (cs_long_t order, const cs_dl *A, double *b) ; 216 cs_dl *cs_dl_transpose (const cs_dl *A, cs_long_t values) ; 217 cs_dl *cs_dl_compress (const cs_dl *T) ; 218 double cs_dl_norm (const cs_dl *A) ; 219 cs_long_t cs_dl_print (const cs_dl *A, cs_long_t brief) ; 220 cs_dl *cs_dl_load (FILE *f) ; 221 222 /* utilities */ 223 void *cs_dl_calloc (cs_long_t n, size_t size) ; 224 void *cs_dl_free (void *p) ; 225 void *cs_dl_realloc (void *p, cs_long_t n, size_t size, cs_long_t *ok) ; 226 cs_dl *cs_dl_spalloc (cs_long_t m, cs_long_t n, cs_long_t nzmax, cs_long_t values, 227 cs_long_t t) ; 228 cs_dl *cs_dl_spfree (cs_dl *A) ; 229 cs_long_t cs_dl_sprealloc (cs_dl *A, cs_long_t nzmax) ; 230 void *cs_dl_malloc (cs_long_t n, size_t size) ; 231 232 /* --- secondary CSparse routines and data structures ----------------------- */ 233 234 typedef struct cs_dl_symbolic /* symbolic Cholesky, LU, or QR analysis */ 235 { 236 cs_long_t *pinv ; /* inverse row perm. for QR, fill red. perm for Chol */ 237 cs_long_t *q ; /* fill-reducing column permutation for LU and QR */ 238 cs_long_t *parent ; /* elimination tree for Cholesky and QR */ 239 cs_long_t *cp ; /* column pointers for Cholesky, row counts for QR */ 240 cs_long_t *leftmost ; /* leftmost[i] = min(find(A(i,:))), for QR */ 241 cs_long_t m2 ; /* # of rows for QR, after adding fictitious rows */ 242 double lnz ; /* # entries in L for LU or Cholesky; in V for QR */ 243 double unz ; /* # entries in U for LU; in R for QR */ 244 } cs_dls ; 245 246 typedef struct cs_dl_numeric /* numeric Cholesky, LU, or QR factorization */ 247 { 248 cs_dl *L ; /* L for LU and Cholesky, V for QR */ 249 cs_dl *U ; /* U for LU, r for QR, not used for Cholesky */ 250 cs_long_t *pinv ; /* partial pivoting for LU */ 251 double *B ; /* beta [0..n-1] for QR */ 252 } cs_dln ; 253 254 typedef struct cs_dl_dmperm_results /* cs_dl_dmperm or cs_dl_scc output */ 255 { 256 cs_long_t *p ; /* size m, row permutation */ 257 cs_long_t *q ; /* size n, column permutation */ 258 cs_long_t *r ; /* size nb+1, block k is rows r[k] to r[k+1]-1 in A(p,q) */ 259 cs_long_t *s ; /* size nb+1, block k is cols s[k] to s[k+1]-1 in A(p,q) */ 260 cs_long_t nb ; /* # of blocks in fine dmperm decomposition */ 261 cs_long_t rr [5] ; /* coarse row decomposition */ 262 cs_long_t cc [5] ; /* coarse column decomposition */ 263 } cs_dld ; 264 265 cs_long_t *cs_dl_amd (cs_long_t order, const cs_dl *A) ; 266 cs_dln *cs_dl_chol (const cs_dl *A, const cs_dls *S) ; 267 cs_dld *cs_dl_dmperm (const cs_dl *A, cs_long_t seed) ; 268 cs_long_t cs_dl_droptol (cs_dl *A, double tol) ; 269 cs_long_t cs_dl_dropzeros (cs_dl *A) ; 270 cs_long_t cs_dl_happly (const cs_dl *V, cs_long_t i, double beta, double *x) ; 271 cs_long_t cs_dl_ipvec (const cs_long_t *p, const double *b, double *x, cs_long_t n) ; 272 cs_long_t cs_dl_lsolve (const cs_dl *L, double *x) ; 273 cs_long_t cs_dl_ltsolve (const cs_dl *L, double *x) ; 274 cs_dln *cs_dl_lu (const cs_dl *A, const cs_dls *S, double tol) ; 275 cs_dl *cs_dl_permute (const cs_dl *A, const cs_long_t *pinv, const cs_long_t *q, 276 cs_long_t values) ; 277 cs_long_t *cs_dl_pinv (const cs_long_t *p, cs_long_t n) ; 278 cs_long_t cs_dl_pvec (const cs_long_t *p, const double *b, double *x, cs_long_t n) ; 279 cs_dln *cs_dl_qr (const cs_dl *A, const cs_dls *S) ; 280 cs_dls *cs_dl_schol (cs_long_t order, const cs_dl *A) ; 281 cs_dls *cs_dl_sqr (cs_long_t order, const cs_dl *A, cs_long_t qr) ; 282 cs_dl *cs_dl_symperm (const cs_dl *A, const cs_long_t *pinv, cs_long_t values) ; 283 cs_long_t cs_dl_usolve (const cs_dl *U, double *x) ; 284 cs_long_t cs_dl_utsolve (const cs_dl *U, double *x) ; 285 cs_long_t cs_dl_updown (cs_dl *L, cs_long_t sigma, const cs_dl *C, 286 const cs_long_t *parent) ; 287 288 /* utilities */ 289 cs_dls *cs_dl_sfree (cs_dls *S) ; 290 cs_dln *cs_dl_nfree (cs_dln *N) ; 291 cs_dld *cs_dl_dfree (cs_dld *D) ; 292 293 /* --- tertiary CSparse routines -------------------------------------------- */ 294 295 cs_long_t *cs_dl_counts (const cs_dl *A, const cs_long_t *parent, 296 const cs_long_t *post, cs_long_t ata) ; 297 double cs_dl_cumsum (cs_long_t *p, cs_long_t *c, cs_long_t n) ; 298 cs_long_t cs_dl_dfs (cs_long_t j, cs_dl *G, cs_long_t top, cs_long_t *xi, 299 cs_long_t *pstack, const cs_long_t *pinv) ; 300 cs_long_t *cs_dl_etree (const cs_dl *A, cs_long_t ata) ; 301 cs_long_t cs_dl_fkeep (cs_dl *A, 302 cs_long_t (*fkeep) (cs_long_t, cs_long_t, double, void *), void *other) ; 303 double cs_dl_house (double *x, double *beta, cs_long_t n) ; 304 cs_long_t *cs_dl_maxtrans (const cs_dl *A, cs_long_t seed) ; 305 cs_long_t *cs_dl_post (const cs_long_t *parent, cs_long_t n) ; 306 cs_dld *cs_dl_scc (cs_dl *A) ; 307 cs_long_t cs_dl_scatter (const cs_dl *A, cs_long_t j, double beta, cs_long_t *w, 308 double *x, cs_long_t mark,cs_dl *C, cs_long_t nz) ; 309 cs_long_t cs_dl_tdfs (cs_long_t j, cs_long_t k, cs_long_t *head, const cs_long_t *next, 310 cs_long_t *post, cs_long_t *stack) ; 311 cs_long_t cs_dl_leaf (cs_long_t i, cs_long_t j, const cs_long_t *first, 312 cs_long_t *maxfirst, cs_long_t *prevleaf, cs_long_t *ancestor, cs_long_t *jleaf) ; 313 cs_long_t cs_dl_reach (cs_dl *G, const cs_dl *B, cs_long_t k, cs_long_t *xi, 314 const cs_long_t *pinv) ; 315 cs_long_t cs_dl_spsolve (cs_dl *L, const cs_dl *B, cs_long_t k, cs_long_t *xi, 316 double *x, const cs_long_t *pinv, cs_long_t lo) ; 317 cs_long_t cs_dl_ereach (const cs_dl *A, cs_long_t k, const cs_long_t *parent, 318 cs_long_t *s, cs_long_t *w) ; 319 cs_long_t *cs_dl_randperm (cs_long_t n, cs_long_t seed) ; 320 321 /* utilities */ 322 cs_dld *cs_dl_dalloc (cs_long_t m, cs_long_t n) ; 323 cs_dl *cs_dl_done (cs_dl *C, void *w, void *x, cs_long_t ok) ; 324 cs_long_t *cs_dl_idone (cs_long_t *p, cs_dl *C, void *w, cs_long_t ok) ; 325 cs_dln *cs_dl_ndone (cs_dln *N, cs_dl *C, void *w, void *x, cs_long_t ok) ; 326 cs_dld *cs_dl_ddone (cs_dld *D, cs_dl *C, void *w, cs_long_t ok) ; 327 328 329 /* -------------------------------------------------------------------------- */ 330 /* complex/int version of CXSparse */ 331 /* -------------------------------------------------------------------------- */ 332 333 #ifndef NCOMPLEX 334 335 /* --- primary CSparse routines and data structures ------------------------- */ 336 337 typedef struct cs_ci_sparse /* matrix in compressed-column or triplet form */ 338 { 339 int nzmax ; /* maximum number of entries */ 340 int m ; /* number of rows */ 341 int n ; /* number of columns */ 342 int *p ; /* column pointers (size n+1) or col indices (size nzmax) */ 343 int *i ; /* row indices, size nzmax */ 344 cs_complex_t *x ; /* numerical values, size nzmax */ 345 int nz ; /* # of entries in triplet matrix, -1 for compressed-col */ 346 } cs_ci ; 347 348 cs_ci *cs_ci_add (const cs_ci *A, const cs_ci *B, cs_complex_t alpha, 349 cs_complex_t beta) ; 350 int cs_ci_cholsol (int order, const cs_ci *A, cs_complex_t *b) ; 351 int cs_ci_dupl (cs_ci *A) ; 352 int cs_ci_entry (cs_ci *T, int i, int j, cs_complex_t x) ; 353 int cs_ci_lusol (int order, const cs_ci *A, cs_complex_t *b, double tol) ; 354 int cs_ci_gaxpy (const cs_ci *A, const cs_complex_t *x, cs_complex_t *y) ; 355 cs_ci *cs_ci_multiply (const cs_ci *A, const cs_ci *B) ; 356 int cs_ci_qrsol (int order, const cs_ci *A, cs_complex_t *b) ; 357 cs_ci *cs_ci_transpose (const cs_ci *A, int values) ; 358 cs_ci *cs_ci_compress (const cs_ci *T) ; 359 double cs_ci_norm (const cs_ci *A) ; 360 int cs_ci_print (const cs_ci *A, int brief) ; 361 cs_ci *cs_ci_load (FILE *f) ; 362 363 /* utilities */ 364 void *cs_ci_calloc (int n, size_t size) ; 365 void *cs_ci_free (void *p) ; 366 void *cs_ci_realloc (void *p, int n, size_t size, int *ok) ; 367 cs_ci *cs_ci_spalloc (int m, int n, int nzmax, int values, int t) ; 368 cs_ci *cs_ci_spfree (cs_ci *A) ; 369 int cs_ci_sprealloc (cs_ci *A, int nzmax) ; 370 void *cs_ci_malloc (int n, size_t size) ; 371 372 /* --- secondary CSparse routines and data structures ----------------------- */ 373 374 typedef struct cs_ci_symbolic /* symbolic Cholesky, LU, or QR analysis */ 375 { 376 int *pinv ; /* inverse row perm. for QR, fill red. perm for Chol */ 377 int *q ; /* fill-reducing column permutation for LU and QR */ 378 int *parent ; /* elimination tree for Cholesky and QR */ 379 int *cp ; /* column pointers for Cholesky, row counts for QR */ 380 int *leftmost ; /* leftmost[i] = min(find(A(i,:))), for QR */ 381 int m2 ; /* # of rows for QR, after adding fictitious rows */ 382 double lnz ; /* # entries in L for LU or Cholesky; in V for QR */ 383 double unz ; /* # entries in U for LU; in R for QR */ 384 } cs_cis ; 385 386 typedef struct cs_ci_numeric /* numeric Cholesky, LU, or QR factorization */ 387 { 388 cs_ci *L ; /* L for LU and Cholesky, V for QR */ 389 cs_ci *U ; /* U for LU, r for QR, not used for Cholesky */ 390 int *pinv ; /* partial pivoting for LU */ 391 double *B ; /* beta [0..n-1] for QR */ 392 } cs_cin ; 393 394 typedef struct cs_ci_dmperm_results /* cs_ci_dmperm or cs_ci_scc output */ 395 { 396 int *p ; /* size m, row permutation */ 397 int *q ; /* size n, column permutation */ 398 int *r ; /* size nb+1, block k is rows r[k] to r[k+1]-1 in A(p,q) */ 399 int *s ; /* size nb+1, block k is cols s[k] to s[k+1]-1 in A(p,q) */ 400 int nb ; /* # of blocks in fine dmperm decomposition */ 401 int rr [5] ; /* coarse row decomposition */ 402 int cc [5] ; /* coarse column decomposition */ 403 } cs_cid ; 404 405 int *cs_ci_amd (int order, const cs_ci *A) ; 406 cs_cin *cs_ci_chol (const cs_ci *A, const cs_cis *S) ; 407 cs_cid *cs_ci_dmperm (const cs_ci *A, int seed) ; 408 int cs_ci_droptol (cs_ci *A, double tol) ; 409 int cs_ci_dropzeros (cs_ci *A) ; 410 int cs_ci_happly (const cs_ci *V, int i, double beta, cs_complex_t *x) ; 411 int cs_ci_ipvec (const int *p, const cs_complex_t *b, cs_complex_t *x, int n) ; 412 int cs_ci_lsolve (const cs_ci *L, cs_complex_t *x) ; 413 int cs_ci_ltsolve (const cs_ci *L, cs_complex_t *x) ; 414 cs_cin *cs_ci_lu (const cs_ci *A, const cs_cis *S, double tol) ; 415 cs_ci *cs_ci_permute (const cs_ci *A, const int *pinv, const int *q, 416 int values) ; 417 int *cs_ci_pinv (const int *p, int n) ; 418 int cs_ci_pvec (const int *p, const cs_complex_t *b, cs_complex_t *x, int n) ; 419 cs_cin *cs_ci_qr (const cs_ci *A, const cs_cis *S) ; 420 cs_cis *cs_ci_schol (int order, const cs_ci *A) ; 421 cs_cis *cs_ci_sqr (int order, const cs_ci *A, int qr) ; 422 cs_ci *cs_ci_symperm (const cs_ci *A, const int *pinv, int values) ; 423 int cs_ci_usolve (const cs_ci *U, cs_complex_t *x) ; 424 int cs_ci_utsolve (const cs_ci *U, cs_complex_t *x) ; 425 int cs_ci_updown (cs_ci *L, int sigma, const cs_ci *C, const int *parent) ; 426 427 /* utilities */ 428 cs_cis *cs_ci_sfree (cs_cis *S) ; 429 cs_cin *cs_ci_nfree (cs_cin *N) ; 430 cs_cid *cs_ci_dfree (cs_cid *D) ; 431 432 /* --- tertiary CSparse routines -------------------------------------------- */ 433 434 int *cs_ci_counts (const cs_ci *A, const int *parent, const int *post, 435 int ata) ; 436 double cs_ci_cumsum (int *p, int *c, int n) ; 437 int cs_ci_dfs (int j, cs_ci *G, int top, int *xi, int *pstack, 438 const int *pinv) ; 439 int *cs_ci_etree (const cs_ci *A, int ata) ; 440 int cs_ci_fkeep (cs_ci *A, int (*fkeep) (int, int, cs_complex_t, void *), 441 void *other) ; 442 cs_complex_t cs_ci_house (cs_complex_t *x, double *beta, int n) ; 443 int *cs_ci_maxtrans (const cs_ci *A, int seed) ; 444 int *cs_ci_post (const int *parent, int n) ; 445 cs_cid *cs_ci_scc (cs_ci *A) ; 446 int cs_ci_scatter (const cs_ci *A, int j, cs_complex_t beta, int *w, 447 cs_complex_t *x, int mark,cs_ci *C, int nz) ; 448 int cs_ci_tdfs (int j, int k, int *head, const int *next, int *post, 449 int *stack) ; 450 int cs_ci_leaf (int i, int j, const int *first, int *maxfirst, int *prevleaf, 451 int *ancestor, int *jleaf) ; 452 int cs_ci_reach (cs_ci *G, const cs_ci *B, int k, int *xi, const int *pinv) ; 453 int cs_ci_spsolve (cs_ci *L, const cs_ci *B, int k, int *xi, 454 cs_complex_t *x, const int *pinv, int lo) ; 455 int cs_ci_ereach (const cs_ci *A, int k, const int *parent, int *s, int *w) ; 456 int *cs_ci_randperm (int n, int seed) ; 457 458 /* utilities */ 459 cs_cid *cs_ci_dalloc (int m, int n) ; 460 cs_ci *cs_ci_done (cs_ci *C, void *w, void *x, int ok) ; 461 int *cs_ci_idone (int *p, cs_ci *C, void *w, int ok) ; 462 cs_cin *cs_ci_ndone (cs_cin *N, cs_ci *C, void *w, void *x, int ok) ; 463 cs_cid *cs_ci_ddone (cs_cid *D, cs_ci *C, void *w, int ok) ; 464 465 466 /* -------------------------------------------------------------------------- */ 467 /* complex/cs_long_t version of CXSparse */ 468 /* -------------------------------------------------------------------------- */ 469 470 /* --- primary CSparse routines and data structures ------------------------- */ 471 472 typedef struct cs_cl_sparse /* matrix in compressed-column or triplet form */ 473 { 474 cs_long_t nzmax ; /* maximum number of entries */ 475 cs_long_t m ; /* number of rows */ 476 cs_long_t n ; /* number of columns */ 477 cs_long_t *p ; /* column pointers (size n+1) or col indlces (size nzmax) */ 478 cs_long_t *i ; /* row indices, size nzmax */ 479 cs_complex_t *x ; /* numerical values, size nzmax */ 480 cs_long_t nz ; /* # of entries in triplet matrix, -1 for compressed-col */ 481 } cs_cl ; 482 483 cs_cl *cs_cl_add (const cs_cl *A, const cs_cl *B, cs_complex_t alpha, 484 cs_complex_t beta) ; 485 cs_long_t cs_cl_cholsol (cs_long_t order, const cs_cl *A, cs_complex_t *b) ; 486 cs_long_t cs_cl_dupl (cs_cl *A) ; 487 cs_long_t cs_cl_entry (cs_cl *T, cs_long_t i, cs_long_t j, cs_complex_t x) ; 488 cs_long_t cs_cl_lusol (cs_long_t order, const cs_cl *A, cs_complex_t *b, 489 double tol) ; 490 cs_long_t cs_cl_gaxpy (const cs_cl *A, const cs_complex_t *x, cs_complex_t *y) ; 491 cs_cl *cs_cl_multiply (const cs_cl *A, const cs_cl *B) ; 492 cs_long_t cs_cl_qrsol (cs_long_t order, const cs_cl *A, cs_complex_t *b) ; 493 cs_cl *cs_cl_transpose (const cs_cl *A, cs_long_t values) ; 494 cs_cl *cs_cl_compress (const cs_cl *T) ; 495 double cs_cl_norm (const cs_cl *A) ; 496 cs_long_t cs_cl_print (const cs_cl *A, cs_long_t brief) ; 497 cs_cl *cs_cl_load (FILE *f) ; 498 499 /* utilities */ 500 void *cs_cl_calloc (cs_long_t n, size_t size) ; 501 void *cs_cl_free (void *p) ; 502 void *cs_cl_realloc (void *p, cs_long_t n, size_t size, cs_long_t *ok) ; 503 cs_cl *cs_cl_spalloc (cs_long_t m, cs_long_t n, cs_long_t nzmax, cs_long_t values, 504 cs_long_t t) ; 505 cs_cl *cs_cl_spfree (cs_cl *A) ; 506 cs_long_t cs_cl_sprealloc (cs_cl *A, cs_long_t nzmax) ; 507 void *cs_cl_malloc (cs_long_t n, size_t size) ; 508 509 /* --- secondary CSparse routines and data structures ----------------------- */ 510 511 typedef struct cs_cl_symbolic /* symbolic Cholesky, LU, or QR analysis */ 512 { 513 cs_long_t *pinv ; /* inverse row perm. for QR, fill red. perm for Chol */ 514 cs_long_t *q ; /* fill-reducing column permutation for LU and QR */ 515 cs_long_t *parent ; /* elimination tree for Cholesky and QR */ 516 cs_long_t *cp ; /* column pointers for Cholesky, row counts for QR */ 517 cs_long_t *leftmost ; /* leftmost[i] = min(find(A(i,:))), for QR */ 518 cs_long_t m2 ; /* # of rows for QR, after adding fictitious rows */ 519 double lnz ; /* # entries in L for LU or Cholesky; in V for QR */ 520 double unz ; /* # entries in U for LU; in R for QR */ 521 } cs_cls ; 522 523 typedef struct cs_cl_numeric /* numeric Cholesky, LU, or QR factorization */ 524 { 525 cs_cl *L ; /* L for LU and Cholesky, V for QR */ 526 cs_cl *U ; /* U for LU, r for QR, not used for Cholesky */ 527 cs_long_t *pinv ; /* partial pivoting for LU */ 528 double *B ; /* beta [0..n-1] for QR */ 529 } cs_cln ; 530 531 typedef struct cs_cl_dmperm_results /* cs_cl_dmperm or cs_cl_scc output */ 532 { 533 cs_long_t *p ; /* size m, row permutation */ 534 cs_long_t *q ; /* size n, column permutation */ 535 cs_long_t *r ; /* size nb+1, block k is rows r[k] to r[k+1]-1 in A(p,q) */ 536 cs_long_t *s ; /* size nb+1, block k is cols s[k] to s[k+1]-1 in A(p,q) */ 537 cs_long_t nb ; /* # of blocks in fine dmperm decomposition */ 538 cs_long_t rr [5] ; /* coarse row decomposition */ 539 cs_long_t cc [5] ; /* coarse column decomposition */ 540 } cs_cld ; 541 542 cs_long_t *cs_cl_amd (cs_long_t order, const cs_cl *A) ; 543 cs_cln *cs_cl_chol (const cs_cl *A, const cs_cls *S) ; 544 cs_cld *cs_cl_dmperm (const cs_cl *A, cs_long_t seed) ; 545 cs_long_t cs_cl_droptol (cs_cl *A, double tol) ; 546 cs_long_t cs_cl_dropzeros (cs_cl *A) ; 547 cs_long_t cs_cl_happly (const cs_cl *V, cs_long_t i, double beta, cs_complex_t *x) ; 548 cs_long_t cs_cl_ipvec (const cs_long_t *p, const cs_complex_t *b, 549 cs_complex_t *x, cs_long_t n) ; 550 cs_long_t cs_cl_lsolve (const cs_cl *L, cs_complex_t *x) ; 551 cs_long_t cs_cl_ltsolve (const cs_cl *L, cs_complex_t *x) ; 552 cs_cln *cs_cl_lu (const cs_cl *A, const cs_cls *S, double tol) ; 553 cs_cl *cs_cl_permute (const cs_cl *A, const cs_long_t *pinv, const cs_long_t *q, 554 cs_long_t values) ; 555 cs_long_t *cs_cl_pinv (const cs_long_t *p, cs_long_t n) ; 556 cs_long_t cs_cl_pvec (const cs_long_t *p, const cs_complex_t *b, 557 cs_complex_t *x, cs_long_t n) ; 558 cs_cln *cs_cl_qr (const cs_cl *A, const cs_cls *S) ; 559 cs_cls *cs_cl_schol (cs_long_t order, const cs_cl *A) ; 560 cs_cls *cs_cl_sqr (cs_long_t order, const cs_cl *A, cs_long_t qr) ; 561 cs_cl *cs_cl_symperm (const cs_cl *A, const cs_long_t *pinv, cs_long_t values) ; 562 cs_long_t cs_cl_usolve (const cs_cl *U, cs_complex_t *x) ; 563 cs_long_t cs_cl_utsolve (const cs_cl *U, cs_complex_t *x) ; 564 cs_long_t cs_cl_updown (cs_cl *L, cs_long_t sigma, const cs_cl *C, 565 const cs_long_t *parent) ; 566 567 /* utilities */ 568 cs_cls *cs_cl_sfree (cs_cls *S) ; 569 cs_cln *cs_cl_nfree (cs_cln *N) ; 570 cs_cld *cs_cl_dfree (cs_cld *D) ; 571 572 /* --- tertiary CSparse routines -------------------------------------------- */ 573 574 cs_long_t *cs_cl_counts (const cs_cl *A, const cs_long_t *parent, 575 const cs_long_t *post, cs_long_t ata) ; 576 double cs_cl_cumsum (cs_long_t *p, cs_long_t *c, cs_long_t n) ; 577 cs_long_t cs_cl_dfs (cs_long_t j, cs_cl *G, cs_long_t top, cs_long_t *xi, 578 cs_long_t *pstack, const cs_long_t *pinv) ; 579 cs_long_t *cs_cl_etree (const cs_cl *A, cs_long_t ata) ; 580 cs_long_t cs_cl_fkeep (cs_cl *A, 581 cs_long_t (*fkeep) (cs_long_t, cs_long_t, cs_complex_t, void *), void *other) ; 582 cs_complex_t cs_cl_house (cs_complex_t *x, double *beta, cs_long_t n) ; 583 cs_long_t *cs_cl_maxtrans (const cs_cl *A, cs_long_t seed) ; 584 cs_long_t *cs_cl_post (const cs_long_t *parent, cs_long_t n) ; 585 cs_cld *cs_cl_scc (cs_cl *A) ; 586 cs_long_t cs_cl_scatter (const cs_cl *A, cs_long_t j, cs_complex_t beta, 587 cs_long_t *w, cs_complex_t *x, cs_long_t mark,cs_cl *C, cs_long_t nz) ; 588 cs_long_t cs_cl_tdfs (cs_long_t j, cs_long_t k, cs_long_t *head, const cs_long_t *next, 589 cs_long_t *post, cs_long_t *stack) ; 590 cs_long_t cs_cl_leaf (cs_long_t i, cs_long_t j, const cs_long_t *first, 591 cs_long_t *maxfirst, cs_long_t *prevleaf, cs_long_t *ancestor, cs_long_t *jleaf) ; 592 cs_long_t cs_cl_reach (cs_cl *G, const cs_cl *B, cs_long_t k, cs_long_t *xi, 593 const cs_long_t *pinv) ; 594 cs_long_t cs_cl_spsolve (cs_cl *L, const cs_cl *B, cs_long_t k, cs_long_t *xi, 595 cs_complex_t *x, const cs_long_t *pinv, cs_long_t lo) ; 596 cs_long_t cs_cl_ereach (const cs_cl *A, cs_long_t k, const cs_long_t *parent, 597 cs_long_t *s, cs_long_t *w) ; 598 cs_long_t *cs_cl_randperm (cs_long_t n, cs_long_t seed) ; 599 600 /* utilities */ 601 cs_cld *cs_cl_dalloc (cs_long_t m, cs_long_t n) ; 602 cs_cl *cs_cl_done (cs_cl *C, void *w, void *x, cs_long_t ok) ; 603 cs_long_t *cs_cl_idone (cs_long_t *p, cs_cl *C, void *w, cs_long_t ok) ; 604 cs_cln *cs_cl_ndone (cs_cln *N, cs_cl *C, void *w, void *x, cs_long_t ok) ; 605 cs_cld *cs_cl_ddone (cs_cld *D, cs_cl *C, void *w, cs_long_t ok) ; 606 607 #endif 608 609 /* -------------------------------------------------------------------------- */ 610 /* Macros for constructing each version of CSparse */ 611 /* -------------------------------------------------------------------------- */ 612 613 #ifdef CS_LONG 614 #define CS_INT cs_long_t 615 #define CS_INT_MAX cs_long_t_max 616 #define CS_ID cs_long_t_id 617 #ifdef CS_COMPLEX 618 #define CS_ENTRY cs_complex_t 619 #define CS_NAME(nm) cs_cl ## nm 620 #define cs cs_cl 621 #else 622 #define CS_ENTRY double 623 #define CS_NAME(nm) cs_dl ## nm 624 #define cs cs_dl 625 #endif 626 #else 627 #define CS_INT int 628 #define CS_INT_MAX INT_MAX 629 #define CS_ID "%d" 630 #ifdef CS_COMPLEX 631 #define CS_ENTRY cs_complex_t 632 #define CS_NAME(nm) cs_ci ## nm 633 #define cs cs_ci 634 #else 635 #define CS_ENTRY double 636 #define CS_NAME(nm) cs_di ## nm 637 #define cs cs_di 638 #endif 639 #endif 640 641 #ifdef CS_COMPLEX 642 #define CS_REAL(x) creal(x) 643 #define CS_IMAG(x) cimag(x) 644 #define CS_CONJ(x) conj(x) 645 #define CS_ABS(x) cabs(x) 646 #else 647 #define CS_REAL(x) (x) 648 #define CS_IMAG(x) (0.) 649 #define CS_CONJ(x) (x) 650 #define CS_ABS(x) fabs(x) 651 #endif 652 653 #define CS_MAX(a,b) (((a) > (b)) ? (a) : (b)) 654 #define CS_MIN(a,b) (((a) < (b)) ? (a) : (b)) 655 #define CS_FLIP(i) (-(i)-2) 656 #define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i)) 657 #define CS_MARKED(w,j) (w [j] < 0) 658 #define CS_MARK(w,j) { w [j] = CS_FLIP (w [j]) ; } 659 #define CS_CSC(A) (A && (A->nz == -1)) 660 #define CS_TRIPLET(A) (A && (A->nz >= 0)) 661 662 /* --- primary CSparse routines and data structures ------------------------- */ 663 664 #define cs_add CS_NAME (_add) 665 #define cs_cholsol CS_NAME (_cholsol) 666 #define cs_dupl CS_NAME (_dupl) 667 #define cs_entry CS_NAME (_entry) 668 #define cs_lusol CS_NAME (_lusol) 669 #define cs_gaxpy CS_NAME (_gaxpy) 670 #define cs_multiply CS_NAME (_multiply) 671 #define cs_qrsol CS_NAME (_qrsol) 672 #define cs_transpose CS_NAME (_transpose) 673 #define cs_compress CS_NAME (_compress) 674 #define cs_norm CS_NAME (_norm) 675 #define cs_print CS_NAME (_print) 676 #define cs_load CS_NAME (_load) 677 678 /* utilities */ 679 #define cs_calloc CS_NAME (_calloc) 680 #define cs_free CS_NAME (_free) 681 #define cs_realloc CS_NAME (_realloc) 682 #define cs_spalloc CS_NAME (_spalloc) 683 #define cs_spfree CS_NAME (_spfree) 684 #define cs_sprealloc CS_NAME (_sprealloc) 685 #define cs_malloc CS_NAME (_malloc) 686 687 /* --- secondary CSparse routines and data structures ----------------------- */ 688 #define css CS_NAME (s) 689 #define csn CS_NAME (n) 690 #define csd CS_NAME (d) 691 692 #define cs_amd CS_NAME (_amd) 693 #define cs_chol CS_NAME (_chol) 694 #define cs_dmperm CS_NAME (_dmperm) 695 #define cs_droptol CS_NAME (_droptol) 696 #define cs_dropzeros CS_NAME (_dropzeros) 697 #define cs_happly CS_NAME (_happly) 698 #define cs_ipvec CS_NAME (_ipvec) 699 #define cs_lsolve CS_NAME (_lsolve) 700 #define cs_ltsolve CS_NAME (_ltsolve) 701 #define cs_lu CS_NAME (_lu) 702 #define cs_permute CS_NAME (_permute) 703 #define cs_pinv CS_NAME (_pinv) 704 #define cs_pvec CS_NAME (_pvec) 705 #define cs_qr CS_NAME (_qr) 706 #define cs_schol CS_NAME (_schol) 707 #define cs_sqr CS_NAME (_sqr) 708 #define cs_symperm CS_NAME (_symperm) 709 #define cs_usolve CS_NAME (_usolve) 710 #define cs_utsolve CS_NAME (_utsolve) 711 #define cs_updown CS_NAME (_updown) 712 713 /* utilities */ 714 #define cs_sfree CS_NAME (_sfree) 715 #define cs_nfree CS_NAME (_nfree) 716 #define cs_dfree CS_NAME (_dfree) 717 718 /* --- tertiary CSparse routines -------------------------------------------- */ 719 #define cs_counts CS_NAME (_counts) 720 #define cs_cumsum CS_NAME (_cumsum) 721 #define cs_dfs CS_NAME (_dfs) 722 #define cs_etree CS_NAME (_etree) 723 #define cs_fkeep CS_NAME (_fkeep) 724 #define cs_house CS_NAME (_house) 725 #define cs_invmatch CS_NAME (_invmatch) 726 #define cs_maxtrans CS_NAME (_maxtrans) 727 #define cs_post CS_NAME (_post) 728 #define cs_scc CS_NAME (_scc) 729 #define cs_scatter CS_NAME (_scatter) 730 #define cs_tdfs CS_NAME (_tdfs) 731 #define cs_reach CS_NAME (_reach) 732 #define cs_spsolve CS_NAME (_spsolve) 733 #define cs_ereach CS_NAME (_ereach) 734 #define cs_randperm CS_NAME (_randperm) 735 #define cs_leaf CS_NAME (_leaf) 736 737 /* utilities */ 738 #define cs_dalloc CS_NAME (_dalloc) 739 #define cs_done CS_NAME (_done) 740 #define cs_idone CS_NAME (_idone) 741 #define cs_ndone CS_NAME (_ndone) 742 #define cs_ddone CS_NAME (_ddone) 743 744 /* -------------------------------------------------------------------------- */ 745 /* Conversion routines */ 746 /* -------------------------------------------------------------------------- */ 747 748 #ifndef NCOMPLEX 749 cs_di *cs_i_real (cs_ci *A, int real) ; 750 cs_ci *cs_i_complex (cs_di *A, int real) ; 751 cs_dl *cs_l_real (cs_cl *A, cs_long_t real) ; 752 cs_cl *cs_l_complex (cs_dl *A, cs_long_t real) ; 753 #endif 754 755 #ifdef __cplusplus 756 } 757 #endif 758 #endif 759