Matrix  \$Rev: 3071 \$ at \$LastChangedDate: 2015-03-26 15:35:47 +0100 (Thu, 26 Mar 2015) \$
cs.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* C = alpha*A + beta*B */
3 cs *cs_add (const cs *A, const cs *B, double alpha, double beta)
4 {
5  csi p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values ;
6  double *x, *Bx, *Cx ;
7  cs *C ;
8  if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ; /* check inputs */
9  if (A->m != B->m || A->n != B->n) return (NULL) ;
10  m = A->m ; anz = A->p [A->n] ;
11  n = B->n ; Bp = B->p ; Bx = B->x ; bnz = Bp [n] ;
12  w = cs_calloc (m, sizeof (csi)) ; /* get workspace */
13  values = (A->x != NULL) && (Bx != NULL) ;
14  x = values ? cs_malloc (m, sizeof (double)) : NULL ; /* get workspace */
15  C = cs_spalloc (m, n, anz + bnz, values, 0) ; /* allocate result*/
16  if (!C || !w || (values && !x)) return (cs_done (C, w, x, 0)) ;
17  Cp = C->p ; Ci = C->i ; Cx = C->x ;
18  for (j = 0 ; j < n ; j++)
19  {
20  Cp [j] = nz ; /* column j of C starts here */
21  nz = cs_scatter (A, j, alpha, w, x, j+1, C, nz) ; /* alpha*A(:,j)*/
22  nz = cs_scatter (B, j, beta, w, x, j+1, C, nz) ; /* beta*B(:,j) */
23  if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Ci [p]] ;
24  }
25  Cp [n] = nz ; /* finalize the last column of C */
26  cs_sprealloc (C, 0) ; /* remove extra space from C */
27  return (cs_done (C, w, x, 1)) ; /* success; free workspace, return C */
28 }
29 /* clear w */
30 static csi cs_wclear (csi mark, csi lemax, csi *w, csi n)
31 {
32  csi k ;
33  if (mark < 2 || (mark + lemax < 0))
34  {
35  for (k = 0 ; k < n ; k++) if (w [k] != 0) w [k] = 1 ;
36  mark = 2 ;
37  }
38  return (mark) ; /* at this point, w [0..n-1] < mark holds */
39 }
40
41 /* keep off-diagonal entries; drop diagonal entries */
42 static csi cs_diag (csi i, csi j, double aij, void *other) { return (i != j) ; }
43
44 /* p = amd(A+A') if symmetric is true, or amd(A'A) otherwise */
45 csi *cs_amd (csi order, const cs *A) /* order 0:natural, 1:Chol, 2:LU, 3:QR */
46 {
47  cs *C, *A2, *AT ;
48  csi *Cp, *Ci, *last, *W, *len, *nv, *next, *P, *head, *elen, *degree, *w,
49  *hhead, *ATp, *ATi, d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
50  k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
51  ok, cnz, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, n, m, t ;
52  csi h ;
53  /* --- Construct matrix C ----------------------------------------------- */
54  if (!CS_CSC (A) || order <= 0 || order > 3) return (NULL) ; /* check */
55  AT = cs_transpose (A, 0) ; /* compute A' */
56  if (!AT) return (NULL) ;
57  m = A->m ; n = A->n ;
58  dense = CS_MAX (16, 10 * sqrt ((double) n)) ; /* find dense threshold */
59  dense = CS_MIN (n-2, dense) ;
60  if (order == 1 && n == m)
61  {
62  C = cs_add (A, AT, 0, 0) ; /* C = A+A' */
63  }
64  else if (order == 2)
65  {
66  ATp = AT->p ; /* drop dense columns from AT */
67  ATi = AT->i ;
68  for (p2 = 0, j = 0 ; j < m ; j++)
69  {
70  p = ATp [j] ; /* column j of AT starts here */
71  ATp [j] = p2 ; /* new column j starts here */
72  if (ATp [j+1] - p > dense) continue ; /* skip dense col j */
73  for ( ; p < ATp [j+1] ; p++) ATi [p2++] = ATi [p] ;
74  }
75  ATp [m] = p2 ; /* finalize AT */
76  A2 = cs_transpose (AT, 0) ; /* A2 = AT' */
77  C = A2 ? cs_multiply (AT, A2) : NULL ; /* C=A'*A with no dense rows */
78  cs_spfree (A2) ;
79  }
80  else
81  {
82  C = cs_multiply (AT, A) ; /* C=A'*A */
83  }
84  cs_spfree (AT) ;
85  if (!C) return (NULL) ;
86  cs_fkeep (C, &cs_diag, NULL) ; /* drop diagonal entries */
87  Cp = C->p ;
88  cnz = Cp [n] ;
89  P = cs_malloc (n+1, sizeof (csi)) ; /* allocate result */
90  W = cs_malloc (8*(n+1), sizeof (csi)) ; /* get workspace */
91  t = cnz + cnz/5 + 2*n ; /* add elbow room to C */
92  if (!P || !W || !cs_sprealloc (C, t)) return (cs_idone (P, C, W, 0)) ;
93  len = W ; nv = W + (n+1) ; next = W + 2*(n+1) ;
94  head = W + 3*(n+1) ; elen = W + 4*(n+1) ; degree = W + 5*(n+1) ;
95  w = W + 6*(n+1) ; hhead = W + 7*(n+1) ;
96  last = P ; /* use P as workspace for last */
97  /* --- Initialize quotient graph ---------------------------------------- */
98  for (k = 0 ; k < n ; k++) len [k] = Cp [k+1] - Cp [k] ;
99  len [n] = 0 ;
100  nzmax = C->nzmax ;
101  Ci = C->i ;
102  for (i = 0 ; i <= n ; i++)
103  {
104  head [i] = -1 ; /* degree list i is empty */
105  last [i] = -1 ;
106  next [i] = -1 ;
107  hhead [i] = -1 ; /* hash list i is empty */
108  nv [i] = 1 ; /* node i is just one node */
109  w [i] = 1 ; /* node i is alive */
110  elen [i] = 0 ; /* Ek of node i is empty */
111  degree [i] = len [i] ; /* degree of node i */
112  }
113  mark = cs_wclear (0, 0, w, n) ; /* clear w */
114  elen [n] = -2 ; /* n is a dead element */
115  Cp [n] = -1 ; /* n is a root of assembly tree */
116  w [n] = 0 ; /* n is a dead element */
117  /* --- Initialize degree lists ------------------------------------------ */
118  for (i = 0 ; i < n ; i++)
119  {
120  d = degree [i] ;
121  if (d == 0) /* node i is empty */
122  {
123  elen [i] = -2 ; /* element i is dead */
124  nel++ ;
125  Cp [i] = -1 ; /* i is a root of assembly tree */
126  w [i] = 0 ;
127  }
128  else if (d > dense) /* node i is dense */
129  {
130  nv [i] = 0 ; /* absorb i into element n */
131  elen [i] = -1 ; /* node i is dead */
132  nel++ ;
133  Cp [i] = CS_FLIP (n) ;
134  nv [n]++ ;
135  }
136  else
137  {
138  if (head [d] != -1) last [head [d]] = i ;
139  next [i] = head [d] ; /* put node i in degree list d */
140  head [d] = i ;
141  }
142  }
143  while (nel < n) /* while (selecting pivots) do */
144  {
145  /* --- Select node of minimum approximate degree -------------------- */
146  for (k = -1 ; mindeg < n && (k = head [mindeg]) == -1 ; mindeg++) ;
147  if (next [k] != -1) last [next [k]] = -1 ;
148  head [mindeg] = next [k] ; /* remove k from degree list */
149  elenk = elen [k] ; /* elenk = |Ek| */
150  nvk = nv [k] ; /* # of nodes k represents */
151  nel += nvk ; /* nv[k] nodes of A eliminated */
152  /* --- Garbage collection ------------------------------------------- */
153  if (elenk > 0 && cnz + mindeg >= nzmax)
154  {
155  for (j = 0 ; j < n ; j++)
156  {
157  if ((p = Cp [j]) >= 0) /* j is a live node or element */
158  {
159  Cp [j] = Ci [p] ; /* save first entry of object */
160  Ci [p] = CS_FLIP (j) ; /* first entry is now CS_FLIP(j) */
161  }
162  }
163  for (q = 0, p = 0 ; p < cnz ; ) /* scan all of memory */
164  {
165  if ((j = CS_FLIP (Ci [p++])) >= 0) /* found object j */
166  {
167  Ci [q] = Cp [j] ; /* restore first entry of object */
168  Cp [j] = q++ ; /* new pointer to object j */
169  for (k3 = 0 ; k3 < len [j]-1 ; k3++) Ci [q++] = Ci [p++] ;
170  }
171  }
172  cnz = q ; /* Ci [cnz...nzmax-1] now free */
173  }
174  /* --- Construct new element ---------------------------------------- */
175  dk = 0 ;
176  nv [k] = -nvk ; /* flag k as in Lk */
177  p = Cp [k] ;
178  pk1 = (elenk == 0) ? p : cnz ; /* do in place if elen[k] == 0 */
179  pk2 = pk1 ;
180  for (k1 = 1 ; k1 <= elenk + 1 ; k1++)
181  {
182  if (k1 > elenk)
183  {
184  e = k ; /* search the nodes in k */
185  pj = p ; /* list of nodes starts at Ci[pj]*/
186  ln = len [k] - elenk ; /* length of list of nodes in k */
187  }
188  else
189  {
190  e = Ci [p++] ; /* search the nodes in e */
191  pj = Cp [e] ;
192  ln = len [e] ; /* length of list of nodes in e */
193  }
194  for (k2 = 1 ; k2 <= ln ; k2++)
195  {
196  i = Ci [pj++] ;
197  if ((nvi = nv [i]) <= 0) continue ; /* node i dead, or seen */
198  dk += nvi ; /* degree[Lk] += size of node i */
199  nv [i] = -nvi ; /* negate nv[i] to denote i in Lk*/
200  Ci [pk2++] = i ; /* place i in Lk */
201  if (next [i] != -1) last [next [i]] = last [i] ;
202  if (last [i] != -1) /* remove i from degree list */
203  {
204  next [last [i]] = next [i] ;
205  }
206  else
207  {
208  head [degree [i]] = next [i] ;
209  }
210  }
211  if (e != k)
212  {
213  Cp [e] = CS_FLIP (k) ; /* absorb e into k */
214  w [e] = 0 ; /* e is now a dead element */
215  }
216  }
217  if (elenk != 0) cnz = pk2 ; /* Ci [cnz...nzmax] is free */
218  degree [k] = dk ; /* external degree of k - |Lk\i| */
219  Cp [k] = pk1 ; /* element k is in Ci[pk1..pk2-1] */
220  len [k] = pk2 - pk1 ;
221  elen [k] = -2 ; /* k is now an element */
222  /* --- Find set differences ----------------------------------------- */
223  mark = cs_wclear (mark, lemax, w, n) ; /* clear w if necessary */
224  for (pk = pk1 ; pk < pk2 ; pk++) /* scan 1: find |Le\Lk| */
225  {
226  i = Ci [pk] ;
227  if ((eln = elen [i]) <= 0) continue ;/* skip if elen[i] empty */
228  nvi = -nv [i] ; /* nv [i] was negated */
229  wnvi = mark - nvi ;
230  for (p = Cp [i] ; p <= Cp [i] + eln - 1 ; p++) /* scan Ei */
231  {
232  e = Ci [p] ;
233  if (w [e] >= mark)
234  {
235  w [e] -= nvi ; /* decrement |Le\Lk| */
236  }
237  else if (w [e] != 0) /* ensure e is a live element */
238  {
239  w [e] = degree [e] + wnvi ; /* 1st time e seen in scan 1 */
240  }
241  }
242  }
243  /* --- Degree update ------------------------------------------------ */
244  for (pk = pk1 ; pk < pk2 ; pk++) /* scan2: degree update */
245  {
246  i = Ci [pk] ; /* consider node i in Lk */
247  p1 = Cp [i] ;
248  p2 = p1 + elen [i] - 1 ;
249  pn = p1 ;
250  for (h = 0, d = 0, p = p1 ; p <= p2 ; p++) /* scan Ei */
251  {
252  e = Ci [p] ;
253  if (w [e] != 0) /* e is an unabsorbed element */
254  {
255  dext = w [e] - mark ; /* dext = |Le\Lk| */
256  if (dext > 0)
257  {
258  d += dext ; /* sum up the set differences */
259  Ci [pn++] = e ; /* keep e in Ei */
260  h += e ; /* compute the hash of node i */
261  }
262  else
263  {
264  Cp [e] = CS_FLIP (k) ; /* aggressive absorb. e->k */
265  w [e] = 0 ; /* e is a dead element */
266  }
267  }
268  }
269  elen [i] = pn - p1 + 1 ; /* elen[i] = |Ei| */
270  p3 = pn ;
271  p4 = p1 + len [i] ;
272  for (p = p2 + 1 ; p < p4 ; p++) /* prune edges in Ai */
273  {
274  j = Ci [p] ;
275  if ((nvj = nv [j]) <= 0) continue ; /* node j dead or in Lk */
276  d += nvj ; /* degree(i) += |j| */
277  Ci [pn++] = j ; /* place j in node list of i */
278  h += j ; /* compute hash for node i */
279  }
280  if (d == 0) /* check for mass elimination */
281  {
282  Cp [i] = CS_FLIP (k) ; /* absorb i into k */
283  nvi = -nv [i] ;
284  dk -= nvi ; /* |Lk| -= |i| */
285  nvk += nvi ; /* |k| += nv[i] */
286  nel += nvi ;
287  nv [i] = 0 ;
288  elen [i] = -1 ; /* node i is dead */
289  }
290  else
291  {
292  degree [i] = CS_MIN (degree [i], d) ; /* update degree(i) */
293  Ci [pn] = Ci [p3] ; /* move first node to end */
294  Ci [p3] = Ci [p1] ; /* move 1st el. to end of Ei */
295  Ci [p1] = k ; /* add k as 1st element in of Ei */
296  len [i] = pn - p1 + 1 ; /* new len of adj. list of node i */
297  h = ((h<0) ? (-h):h) % n ; /* finalize hash of i */
298  next [i] = hhead [h] ; /* place i in hash bucket */
299  hhead [h] = i ;
300  last [i] = h ; /* save hash of i in last[i] */
301  }
302  } /* scan2 is done */
303  degree [k] = dk ; /* finalize |Lk| */
304  lemax = CS_MAX (lemax, dk) ;
305  mark = cs_wclear (mark+lemax, lemax, w, n) ; /* clear w */
306  /* --- Supernode detection ------------------------------------------ */
307  for (pk = pk1 ; pk < pk2 ; pk++)
308  {
309  i = Ci [pk] ;
310  if (nv [i] >= 0) continue ; /* skip if i is dead */
311  h = last [i] ; /* scan hash bucket of node i */
312  i = hhead [h] ;
313  hhead [h] = -1 ; /* hash bucket will be empty */
314  for ( ; i != -1 && next [i] != -1 ; i = next [i], mark++)
315  {
316  ln = len [i] ;
317  eln = elen [i] ;
318  for (p = Cp [i]+1 ; p <= Cp [i] + ln-1 ; p++) w [Ci [p]] = mark;
319  jlast = i ;
320  for (j = next [i] ; j != -1 ; ) /* compare i with all j */
321  {
322  ok = (len [j] == ln) && (elen [j] == eln) ;
323  for (p = Cp [j] + 1 ; ok && p <= Cp [j] + ln - 1 ; p++)
324  {
325  if (w [Ci [p]] != mark) ok = 0 ; /* compare i and j*/
326  }
327  if (ok) /* i and j are identical */
328  {
329  Cp [j] = CS_FLIP (i) ; /* absorb j into i */
330  nv [i] += nv [j] ;
331  nv [j] = 0 ;
332  elen [j] = -1 ; /* node j is dead */
333  j = next [j] ; /* delete j from hash bucket */
334  next [jlast] = j ;
335  }
336  else
337  {
338  jlast = j ; /* j and i are different */
339  j = next [j] ;
340  }
341  }
342  }
343  }
344  /* --- Finalize new element------------------------------------------ */
345  for (p = pk1, pk = pk1 ; pk < pk2 ; pk++) /* finalize Lk */
346  {
347  i = Ci [pk] ;
348  if ((nvi = -nv [i]) <= 0) continue ;/* skip if i is dead */
349  nv [i] = nvi ; /* restore nv[i] */
350  d = degree [i] + dk - nvi ; /* compute external degree(i) */
351  d = CS_MIN (d, n - nel - nvi) ;
352  if (head [d] != -1) last [head [d]] = i ;
353  next [i] = head [d] ; /* put i back in degree list */
354  last [i] = -1 ;
355  head [d] = i ;
356  mindeg = CS_MIN (mindeg, d) ; /* find new minimum degree */
357  degree [i] = d ;
358  Ci [p++] = i ; /* place i in Lk */
359  }
360  nv [k] = nvk ; /* # nodes absorbed into k */
361  if ((len [k] = p-pk1) == 0) /* length of adj list of element k*/
362  {
363  Cp [k] = -1 ; /* k is a root of the tree */
364  w [k] = 0 ; /* k is now a dead element */
365  }
366  if (elenk != 0) cnz = p ; /* free unused space in Lk */
367  }
368  /* --- Postordering ----------------------------------------------------- */
369  for (i = 0 ; i < n ; i++) Cp [i] = CS_FLIP (Cp [i]) ;/* fix assembly tree */
370  for (j = 0 ; j <= n ; j++) head [j] = -1 ;
371  for (j = n ; j >= 0 ; j--) /* place unordered nodes in lists */
372  {
373  if (nv [j] > 0) continue ; /* skip if j is an element */
374  next [j] = head [Cp [j]] ; /* place j in list of its parent */
375  head [Cp [j]] = j ;
376  }
377  for (e = n ; e >= 0 ; e--) /* place elements in lists */
378  {
379  if (nv [e] <= 0) continue ; /* skip unless e is an element */
380  if (Cp [e] != -1)
381  {
382  next [e] = head [Cp [e]] ; /* place e in list of its parent */
383  head [Cp [e]] = e ;
384  }
385  }
386  for (k = 0, i = 0 ; i <= n ; i++) /* postorder the assembly tree */
387  {
388  if (Cp [i] == -1) k = cs_tdfs (i, k, head, next, P, w) ;
389  }
390  return (cs_idone (P, C, W, 1)) ;
391 }
392 /* L = chol (A, [pinv parent cp]), pinv is optional */
393 csn *cs_chol (const cs *A, const css *S)
394 {
395  double d, lki, *Lx, *x, *Cx ;
396  csi top, i, p, k, n, *Li, *Lp, *cp, *pinv, *s, *c, *parent, *Cp, *Ci ;
397  cs *L, *C, *E ;
398  csn *N ;
399  if (!CS_CSC (A) || !S || !S->cp || !S->parent) return (NULL) ;
400  n = A->n ;
401  N = cs_calloc (1, sizeof (csn)) ; /* allocate result */
402  c = cs_malloc (2*n, sizeof (csi)) ; /* get csi workspace */
403  x = cs_malloc (n, sizeof (double)) ; /* get double workspace */
404  cp = S->cp ; pinv = S->pinv ; parent = S->parent ;
405  C = pinv ? cs_symperm (A, pinv, 1) : ((cs *) A) ;
406  E = pinv ? C : NULL ; /* E is alias for A, or a copy E=A(p,p) */
407  if (!N || !c || !x || !C) return (cs_ndone (N, E, c, x, 0)) ;
408  s = c + n ;
409  Cp = C->p ; Ci = C->i ; Cx = C->x ;
410  N->L = L = cs_spalloc (n, n, cp [n], 1, 0) ; /* allocate result */
411  if (!L) return (cs_ndone (N, E, c, x, 0)) ;
412  Lp = L->p ; Li = L->i ; Lx = L->x ;
413  for (k = 0 ; k < n ; k++) Lp [k] = c [k] = cp [k] ;
414  for (k = 0 ; k < n ; k++) /* compute L(k,:) for L*L' = C */
415  {
416  /* --- Nonzero pattern of L(k,:) ------------------------------------ */
417  top = cs_ereach (C, k, parent, s, c) ; /* find pattern of L(k,:) */
418  x [k] = 0 ; /* x (0:k) is now zero */
419  for (p = Cp [k] ; p < Cp [k+1] ; p++) /* x = full(triu(C(:,k))) */
420  {
421  if (Ci [p] <= k) x [Ci [p]] = Cx [p] ;
422  }
423  d = x [k] ; /* d = C(k,k) */
424  x [k] = 0 ; /* clear x for k+1st iteration */
425  /* --- Triangular solve --------------------------------------------- */
426  for ( ; top < n ; top++) /* solve L(0:k-1,0:k-1) * x = C(:,k) */
427  {
428  i = s [top] ; /* s [top..n-1] is pattern of L(k,:) */
429  lki = x [i] / Lx [Lp [i]] ; /* L(k,i) = x (i) / L(i,i) */
430  x [i] = 0 ; /* clear x for k+1st iteration */
431  for (p = Lp [i] + 1 ; p < c [i] ; p++)
432  {
433  x [Li [p]] -= Lx [p] * lki ;
434  }
435  d -= lki * lki ; /* d = d - L(k,i)*L(k,i) */
436  p = c [i]++ ;
437  Li [p] = k ; /* store L(k,i) in column i */
438  Lx [p] = lki ;
439  }
440  /* --- Compute L(k,k) ----------------------------------------------- */
441  if (d <= 0) return (cs_ndone (N, E, c, x, 0)) ; /* not pos def */
442  p = c [k]++ ;
443  Li [p] = k ; /* store L(k,k) = sqrt (d) in column k */
444  Lx [p] = sqrt (d) ;
445  }
446  Lp [n] = cp [n] ; /* finalize L */
447  return (cs_ndone (N, E, c, x, 1)) ; /* success: free E,s,x; return N */
448 }
449 /* x=A\b where A is symmetric positive definite; b overwritten with solution */
450 csi cs_cholsol (csi order, const cs *A, double *b)
451 {
452  double *x ;
453  css *S ;
454  csn *N ;
455  csi n, ok ;
456  if (!CS_CSC (A) || !b) return (0) ; /* check inputs */
457  n = A->n ;
458  S = cs_schol (order, A) ; /* ordering and symbolic analysis */
459  N = cs_chol (A, S) ; /* numeric Cholesky factorization */
460  x = cs_malloc (n, sizeof (double)) ; /* get workspace */
461  ok = (S && N && x) ;
462  if (ok)
463  {
464  cs_ipvec (S->pinv, b, x, n) ; /* x = P*b */
465  cs_lsolve (N->L, x) ; /* x = L\x */
466  cs_ltsolve (N->L, x) ; /* x = L'\x */
467  cs_pvec (S->pinv, x, b, n) ; /* b = P'*x */
468  }
469  cs_free (x) ;
470  cs_sfree (S) ;
471  cs_nfree (N) ;
472  return (ok) ;
473 }
474 /* C = compressed-column form of a triplet matrix T */
475 cs *cs_compress (const cs *T)
476 {
477  csi m, n, nz, p, k, *Cp, *Ci, *w, *Ti, *Tj ;
478  double *Cx, *Tx ;
479  cs *C ;
480  if (!CS_TRIPLET (T)) return (NULL) ; /* check inputs */
481  m = T->m ; n = T->n ; Ti = T->i ; Tj = T->p ; Tx = T->x ; nz = T->nz ;
482  C = cs_spalloc (m, n, nz, Tx != NULL, 0) ; /* allocate result */
483  w = cs_calloc (n, sizeof (csi)) ; /* get workspace */
484  if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */
485  Cp = C->p ; Ci = C->i ; Cx = C->x ;
486  for (k = 0 ; k < nz ; k++) w [Tj [k]]++ ; /* column counts */
487  cs_cumsum (Cp, w, n) ; /* column pointers */
488  for (k = 0 ; k < nz ; k++)
489  {
490  Ci [p = w [Tj [k]]++] = Ti [k] ; /* A(i,j) is the pth entry in C */
491  if (Cx) Cx [p] = Tx [k] ;
492  }
493  return (cs_done (C, w, NULL, 1)) ; /* success; free w and return C */
494 }
495 /* column counts of LL'=A or LL'=A'A, given parent & post ordering */
497 #define NEXT(J) (ata ? next [J] : -1)
498 static void init_ata (cs *AT, const csi *post, csi *w, csi **head, csi **next)
499 {
500  csi i, k, p, m = AT->n, n = AT->m, *ATp = AT->p, *ATi = AT->i ;
501  *head = w+4*n, *next = w+5*n+1 ;
502  for (k = 0 ; k < n ; k++) w [post [k]] = k ; /* invert post */
503  for (i = 0 ; i < m ; i++)
504  {
505  for (k = n, p = ATp[i] ; p < ATp[i+1] ; p++) k = CS_MIN (k, w [ATi[p]]);
506  (*next) [i] = (*head) [k] ; /* place row i in linked list k */
507  (*head) [k] = i ;
508  }
509 }
510 csi *cs_counts (const cs *A, const csi *parent, const csi *post, csi ata)
511 {
512  csi i, j, k, n, m, J, s, p, q, jleaf, *ATp, *ATi, *maxfirst, *prevleaf,
513  *ancestor, *head = NULL, *next = NULL, *colcount, *w, *first, *delta ;
514  cs *AT ;
515  if (!CS_CSC (A) || !parent || !post) return (NULL) ; /* check inputs */
516  m = A->m ; n = A->n ;
517  s = 4*n + (ata ? (n+m+1) : 0) ;
518  delta = colcount = cs_malloc (n, sizeof (csi)) ; /* allocate result */
519  w = cs_malloc (s, sizeof (csi)) ; /* get workspace */
520  AT = cs_transpose (A, 0) ; /* AT = A' */
521  if (!AT || !colcount || !w) return (cs_idone (colcount, AT, w, 0)) ;
522  ancestor = w ; maxfirst = w+n ; prevleaf = w+2*n ; first = w+3*n ;
523  for (k = 0 ; k < s ; k++) w [k] = -1 ; /* clear workspace w [0..s-1] */
524  for (k = 0 ; k < n ; k++) /* find first [j] */
525  {
526  j = post [k] ;
527  delta [j] = (first [j] == -1) ? 1 : 0 ; /* delta[j]=1 if j is a leaf */
528  for ( ; j != -1 && first [j] == -1 ; j = parent [j]) first [j] = k ;
529  }
530  ATp = AT->p ; ATi = AT->i ;
531  if (ata) init_ata (AT, post, w, &head, &next) ;
532  for (i = 0 ; i < n ; i++) ancestor [i] = i ; /* each node in its own set */
533  for (k = 0 ; k < n ; k++)
534  {
535  j = post [k] ; /* j is the kth node in postordered etree */
536  if (parent [j] != -1) delta [parent [j]]-- ; /* j is not a root */
537  for (J = HEAD (k,j) ; J != -1 ; J = NEXT (J)) /* J=j for LL'=A case */
538  {
539  for (p = ATp [J] ; p < ATp [J+1] ; p++)
540  {
541  i = ATi [p] ;
542  q = cs_leaf (i, j, first, maxfirst, prevleaf, ancestor, &jleaf);
543  if (jleaf >= 1) delta [j]++ ; /* A(i,j) is in skeleton */
544  if (jleaf == 2) delta [q]-- ; /* account for overlap in q */
545  }
546  }
547  if (parent [j] != -1) ancestor [j] = parent [j] ;
548  }
549  for (j = 0 ; j < n ; j++) /* sum up delta's of each child */
550  {
551  if (parent [j] != -1) colcount [parent [j]] += colcount [j] ;
552  }
553  return (cs_idone (colcount, AT, w, 1)) ; /* success: free workspace */
554 }
555 /* p [0..n] = cumulative sum of c [0..n-1], and then copy p [0..n-1] into c */
556 double cs_cumsum (csi *p, csi *c, csi n)
557 {
558  csi i, nz = 0 ;
559  double nz2 = 0 ;
560  if (!p || !c) return (-1) ; /* check inputs */
561  for (i = 0 ; i < n ; i++)
562  {
563  p [i] = nz ;
564  nz += c [i] ;
565  nz2 += c [i] ; /* also in double to avoid csi overflow */
566  c [i] = p [i] ; /* also copy p[0..n-1] back into c[0..n-1]*/
567  }
568  p [n] = nz ;
569  return (nz2) ; /* return sum (c [0..n-1]) */
570 }
571 /* depth-first-search of the graph of a matrix, starting at node j */
572 csi cs_dfs (csi j, cs *G, csi top, csi *xi, csi *pstack, const csi *pinv)
573 {
574  csi i, p, p2, done, jnew, head = 0, *Gp, *Gi ;
575  if (!CS_CSC (G) || !xi || !pstack) return (-1) ; /* check inputs */
576  Gp = G->p ; Gi = G->i ;
577  xi [0] = j ; /* initialize the recursion stack */
579  {
580  j = xi [head] ; /* get j from the top of the recursion stack */
581  jnew = pinv ? (pinv [j]) : j ;
582  if (!CS_MARKED (Gp, j))
583  {
584  CS_MARK (Gp, j) ; /* mark node j as visited */
585  pstack [head] = (jnew < 0) ? 0 : CS_UNFLIP (Gp [jnew]) ;
586  }
587  done = 1 ; /* node j done if no unvisited neighbors */
588  p2 = (jnew < 0) ? 0 : CS_UNFLIP (Gp [jnew+1]) ;
589  for (p = pstack [head] ; p < p2 ; p++) /* examine all neighbors of j */
590  {
591  i = Gi [p] ; /* consider neighbor node i */
592  if (CS_MARKED (Gp, i)) continue ; /* skip visited node i */
593  pstack [head] = p ; /* pause depth-first search of node j */
594  xi [++head] = i ; /* start dfs at node i */
595  done = 0 ; /* node j is not done */
596  break ; /* break, to start dfs (i) */
597  }
598  if (done) /* depth-first search at node j is done */
599  {
600  head-- ; /* remove j from the recursion stack */
601  xi [--top] = j ; /* and place in the output stack */
602  }
603  }
604  return (top) ;
605 }
606 /* breadth-first search for coarse decomposition (C0,C1,R1 or R0,R3,C3) */
607 static csi cs_bfs (const cs *A, csi n, csi *wi, csi *wj, csi *queue,
608  const csi *imatch, const csi *jmatch, csi mark)
609 {
610  csi *Ap, *Ai, head = 0, tail = 0, j, i, p, j2 ;
611  cs *C ;
612  for (j = 0 ; j < n ; j++) /* place all unmatched nodes in queue */
613  {
614  if (imatch [j] >= 0) continue ; /* skip j if matched */
615  wj [j] = 0 ; /* j in set C0 (R0 if transpose) */
616  queue [tail++] = j ; /* place unmatched col j in queue */
617  }
618  if (tail == 0) return (1) ; /* quick return if no unmatched nodes */
619  C = (mark == 1) ? ((cs *) A) : cs_transpose (A, 0) ;
620  if (!C) return (0) ; /* bfs of C=A' to find R3,C3 from R0 */
621  Ap = C->p ; Ai = C->i ;
622  while (head < tail) /* while queue is not empty */
623  {
624  j = queue [head++] ; /* get the head of the queue */
625  for (p = Ap [j] ; p < Ap [j+1] ; p++)
626  {
627  i = Ai [p] ;
628  if (wi [i] >= 0) continue ; /* skip if i is marked */
629  wi [i] = mark ; /* i in set R1 (C3 if transpose) */
630  j2 = jmatch [i] ; /* traverse alternating path to j2 */
631  if (wj [j2] >= 0) continue ;/* skip j2 if it is marked */
632  wj [j2] = mark ; /* j2 in set C1 (R3 if transpose) */
633  queue [tail++] = j2 ; /* add j2 to queue */
634  }
635  }
636  if (mark != 1) cs_spfree (C) ; /* free A' if it was created */
637  return (1) ;
638 }
639
640 /* collect matched rows and columns into p and q */
641 static void cs_matched (csi n, const csi *wj, const csi *imatch, csi *p, csi *q,
642  csi *cc, csi *rr, csi set, csi mark)
643 {
644  csi kc = cc [set], j ;
645  csi kr = rr [set-1] ;
646  for (j = 0 ; j < n ; j++)
647  {
648  if (wj [j] != mark) continue ; /* skip if j is not in C set */
649  p [kr++] = imatch [j] ;
650  q [kc++] = j ;
651  }
652  cc [set+1] = kc ;
653  rr [set] = kr ;
654 }
655
656 /* collect unmatched rows into the permutation vector p */
657 static void cs_unmatched (csi m, const csi *wi, csi *p, csi *rr, csi set)
658 {
659  csi i, kr = rr [set] ;
660  for (i = 0 ; i < m ; i++) if (wi [i] == 0) p [kr++] = i ;
661  rr [set+1] = kr ;
662 }
663
664 /* return 1 if row i is in R2 */
665 static csi cs_rprune (csi i, csi j, double aij, void *other)
666 {
667  csi *rr = (csi *) other ;
668  return (i >= rr [1] && i < rr [2]) ;
669 }
670
671 /* Given A, compute coarse and then fine dmperm */
672 csd *cs_dmperm (const cs *A, csi seed)
673 {
674  csi m, n, i, j, k, cnz, nc, *jmatch, *imatch, *wi, *wj, *pinv, *Cp, *Ci,
675  *ps, *rs, nb1, nb2, *p, *q, *cc, *rr, *r, *s, ok ;
676  cs *C ;
677  csd *D, *scc ;
678  /* --- Maximum matching ------------------------------------------------- */
679  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
680  m = A->m ; n = A->n ;
681  D = cs_dalloc (m, n) ; /* allocate result */
682  if (!D) return (NULL) ;
683  p = D->p ; q = D->q ; r = D->r ; s = D->s ; cc = D->cc ; rr = D->rr ;
684  jmatch = cs_maxtrans (A, seed) ; /* max transversal */
685  imatch = jmatch + m ; /* imatch = inverse of jmatch */
686  if (!jmatch) return (cs_ddone (D, NULL, jmatch, 0)) ;
687  /* --- Coarse decomposition --------------------------------------------- */
688  wi = r ; wj = s ; /* use r and s as workspace */
689  for (j = 0 ; j < n ; j++) wj [j] = -1 ; /* unmark all cols for bfs */
690  for (i = 0 ; i < m ; i++) wi [i] = -1 ; /* unmark all rows for bfs */
691  cs_bfs (A, n, wi, wj, q, imatch, jmatch, 1) ; /* find C1, R1 from C0*/
692  ok = cs_bfs (A, m, wj, wi, p, jmatch, imatch, 3) ; /* find R3, C3 from R0*/
693  if (!ok) return (cs_ddone (D, NULL, jmatch, 0)) ;
694  cs_unmatched (n, wj, q, cc, 0) ; /* unmatched set C0 */
695  cs_matched (n, wj, imatch, p, q, cc, rr, 1, 1) ; /* set R1 and C1 */
696  cs_matched (n, wj, imatch, p, q, cc, rr, 2, -1) ; /* set R2 and C2 */
697  cs_matched (n, wj, imatch, p, q, cc, rr, 3, 3) ; /* set R3 and C3 */
698  cs_unmatched (m, wi, p, rr, 3) ; /* unmatched set R0 */
699  cs_free (jmatch) ;
700  /* --- Fine decomposition ----------------------------------------------- */
701  pinv = cs_pinv (p, m) ; /* pinv=p' */
702  if (!pinv) return (cs_ddone (D, NULL, NULL, 0)) ;
703  C = cs_permute (A, pinv, q, 0) ;/* C=A(p,q) (it will hold A(R2,C2)) */
704  cs_free (pinv) ;
705  if (!C) return (cs_ddone (D, NULL, NULL, 0)) ;
706  Cp = C->p ;
707  nc = cc [3] - cc [2] ; /* delete cols C0, C1, and C3 from C */
708  if (cc [2] > 0) for (j = cc [2] ; j <= cc [3] ; j++) Cp [j-cc[2]] = Cp [j] ;
709  C->n = nc ;
710  if (rr [2] - rr [1] < m) /* delete rows R0, R1, and R3 from C */
711  {
712  cs_fkeep (C, cs_rprune, rr) ;
713  cnz = Cp [nc] ;
714  Ci = C->i ;
715  if (rr [1] > 0) for (k = 0 ; k < cnz ; k++) Ci [k] -= rr [1] ;
716  }
717  C->m = nc ;
718  scc = cs_scc (C) ; /* find strongly connected components of C*/
719  if (!scc) return (cs_ddone (D, C, NULL, 0)) ;
720  /* --- Combine coarse and fine decompositions --------------------------- */
721  ps = scc->p ; /* C(ps,ps) is the permuted matrix */
722  rs = scc->r ; /* kth block is rs[k]..rs[k+1]-1 */
723  nb1 = scc->nb ; /* # of blocks of A(R2,C2) */
724  for (k = 0 ; k < nc ; k++) wj [k] = q [ps [k] + cc [2]] ;
725  for (k = 0 ; k < nc ; k++) q [k + cc [2]] = wj [k] ;
726  for (k = 0 ; k < nc ; k++) wi [k] = p [ps [k] + rr [1]] ;
727  for (k = 0 ; k < nc ; k++) p [k + rr [1]] = wi [k] ;
728  nb2 = 0 ; /* create the fine block partitions */
729  r [0] = s [0] = 0 ;
730  if (cc [2] > 0) nb2++ ; /* leading coarse block A (R1, [C0 C1]) */
731  for (k = 0 ; k < nb1 ; k++) /* coarse block A (R2,C2) */
732  {
733  r [nb2] = rs [k] + rr [1] ; /* A (R2,C2) splits into nb1 fine blocks */
734  s [nb2] = rs [k] + cc [2] ;
735  nb2++ ;
736  }
737  if (rr [2] < m)
738  {
739  r [nb2] = rr [2] ; /* trailing coarse block A ([R3 R0], C3) */
740  s [nb2] = cc [3] ;
741  nb2++ ;
742  }
743  r [nb2] = m ;
744  s [nb2] = n ;
745  D->nb = nb2 ;
746  cs_dfree (scc) ;
747  return (cs_ddone (D, C, NULL, 1)) ;
748 }
749 static csi cs_tol (csi i, csi j, double aij, void *tol)
750 {
751  return (fabs (aij) > *((double *) tol)) ;
752 }
753 csi cs_droptol (cs *A, double tol)
754 {
755  return (cs_fkeep (A, &cs_tol, &tol)) ; /* keep all large entries */
756 }
757 static csi cs_nonzero (csi i, csi j, double aij, void *other)
758 {
759  return (aij != 0) ;
760 }
762 {
763  return (cs_fkeep (A, &cs_nonzero, NULL)) ; /* keep all nonzero entries */
764 }
765 /* remove duplicate entries from A */
767 {
768  csi i, j, p, q, nz = 0, n, m, *Ap, *Ai, *w ;
769  double *Ax ;
770  if (!CS_CSC (A)) return (0) ; /* check inputs */
771  m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
772  w = cs_malloc (m, sizeof (csi)) ; /* get workspace */
773  if (!w) return (0) ; /* out of memory */
774  for (i = 0 ; i < m ; i++) w [i] = -1 ; /* row i not yet seen */
775  for (j = 0 ; j < n ; j++)
776  {
777  q = nz ; /* column j will start at q */
778  for (p = Ap [j] ; p < Ap [j+1] ; p++)
779  {
780  i = Ai [p] ; /* A(i,j) is nonzero */
781  if (w [i] >= q)
782  {
783  Ax [w [i]] += Ax [p] ; /* A(i,j) is a duplicate */
784  }
785  else
786  {
787  w [i] = nz ; /* record where row i occurs */
788  Ai [nz] = i ; /* keep A(i,j) */
789  Ax [nz++] = Ax [p] ;
790  }
791  }
792  Ap [j] = q ; /* record start of column j */
793  }
794  Ap [n] = nz ; /* finalize A */
795  cs_free (w) ; /* free workspace */
796  return (cs_sprealloc (A, 0)) ; /* remove extra space from A */
797 }
798 /* add an entry to a triplet matrix; return 1 if ok, 0 otherwise */
799 csi cs_entry (cs *T, csi i, csi j, double x)
800 {
801  if (!CS_TRIPLET (T) || i < 0 || j < 0) return (0) ; /* check inputs */
802  if (T->nz >= T->nzmax && !cs_sprealloc (T,2*(T->nzmax))) return (0) ;
803  if (T->x) T->x [T->nz] = x ;
804  T->i [T->nz] = i ;
805  T->p [T->nz++] = j ;
806  T->m = CS_MAX (T->m, i+1) ;
807  T->n = CS_MAX (T->n, j+1) ;
808  return (1) ;
809 }
810 /* find nonzero pattern of Cholesky L(k,1:k-1) using etree and triu(A(:,k)) */
811 csi cs_ereach (const cs *A, csi k, const csi *parent, csi *s, csi *w)
812 {
813  csi i, p, n, len, top, *Ap, *Ai ;
814  if (!CS_CSC (A) || !parent || !s || !w) return (-1) ; /* check inputs */
815  top = n = A->n ; Ap = A->p ; Ai = A->i ;
816  CS_MARK (w, k) ; /* mark node k as visited */
817  for (p = Ap [k] ; p < Ap [k+1] ; p++)
818  {
819  i = Ai [p] ; /* A(i,k) is nonzero */
820  if (i > k) continue ; /* only use upper triangular part of A */
821  for (len = 0 ; !CS_MARKED (w,i) ; i = parent [i]) /* traverse up etree*/
822  {
823  s [len++] = i ; /* L(k,i) is nonzero */
824  CS_MARK (w, i) ; /* mark i as visited */
825  }
826  while (len > 0) s [--top] = s [--len] ; /* push path onto stack */
827  }
828  for (p = top ; p < n ; p++) CS_MARK (w, s [p]) ; /* unmark all nodes */
829  CS_MARK (w, k) ; /* unmark node k */
830  return (top) ; /* s [top..n-1] contains pattern of L(k,:)*/
831 }
832 /* compute the etree of A (using triu(A), or A'A without forming A'A */
833 csi *cs_etree (const cs *A, csi ata)
834 {
835  csi i, k, p, m, n, inext, *Ap, *Ai, *w, *parent, *ancestor, *prev ;
836  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
837  m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ;
838  parent = cs_malloc (n, sizeof (csi)) ; /* allocate result */
839  w = cs_malloc (n + (ata ? m : 0), sizeof (csi)) ; /* get workspace */
840  if (!w || !parent) return (cs_idone (parent, NULL, w, 0)) ;
841  ancestor = w ; prev = w + n ;
842  if (ata) for (i = 0 ; i < m ; i++) prev [i] = -1 ;
843  for (k = 0 ; k < n ; k++)
844  {
845  parent [k] = -1 ; /* node k has no parent yet */
846  ancestor [k] = -1 ; /* nor does k have an ancestor */
847  for (p = Ap [k] ; p < Ap [k+1] ; p++)
848  {
849  i = ata ? (prev [Ai [p]]) : (Ai [p]) ;
850  for ( ; i != -1 && i < k ; i = inext) /* traverse from i to k */
851  {
852  inext = ancestor [i] ; /* inext = ancestor of i */
853  ancestor [i] = k ; /* path compression */
854  if (inext == -1) parent [i] = k ; /* no anc., parent is k */
855  }
856  if (ata) prev [Ai [p]] = k ;
857  }
858  }
859  return (cs_idone (parent, NULL, w, 1)) ;
860 }
861 /* drop entries for which fkeep(A(i,j)) is false; return nz if OK, else -1 */
862 csi cs_fkeep (cs *A, csi (*fkeep) (csi, csi, double, void *), void *other)
863 {
864  csi j, p, nz = 0, n, *Ap, *Ai ;
865  double *Ax ;
866  if (!CS_CSC (A) || !fkeep) return (-1) ; /* check inputs */
867  n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
868  for (j = 0 ; j < n ; j++)
869  {
870  p = Ap [j] ; /* get current location of col j */
871  Ap [j] = nz ; /* record new location of col j */
872  for ( ; p < Ap [j+1] ; p++)
873  {
874  if (fkeep (Ai [p], j, Ax ? Ax [p] : 1, other))
875  {
876  if (Ax) Ax [nz] = Ax [p] ; /* keep A(i,j) */
877  Ai [nz++] = Ai [p] ;
878  }
879  }
880  }
881  Ap [n] = nz ; /* finalize A */
882  cs_sprealloc (A, 0) ; /* remove extra space from A */
883  return (nz) ;
884 }
885 /* y = A*x+y */
886 csi cs_gaxpy (const cs *A, const double *x, double *y)
887 {
888  csi p, j, n, *Ap, *Ai ;
889  double *Ax ;
890  if (!CS_CSC (A) || !x || !y) return (0) ; /* check inputs */
891  n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
892  for (j = 0 ; j < n ; j++)
893  {
894  for (p = Ap [j] ; p < Ap [j+1] ; p++)
895  {
896  y [Ai [p]] += Ax [p] * x [j] ;
897  }
898  }
899  return (1) ;
900 }
901 /* apply the ith Householder vector to x */
902 csi cs_happly (const cs *V, csi i, double beta, double *x)
903 {
904  csi p, *Vp, *Vi ;
905  double *Vx, tau = 0 ;
906  if (!CS_CSC (V) || !x) return (0) ; /* check inputs */
907  Vp = V->p ; Vi = V->i ; Vx = V->x ;
908  for (p = Vp [i] ; p < Vp [i+1] ; p++) /* tau = v'*x */
909  {
910  tau += Vx [p] * x [Vi [p]] ;
911  }
912  tau *= beta ; /* tau = beta*(v'*x) */
913  for (p = Vp [i] ; p < Vp [i+1] ; p++) /* x = x - v*tau */
914  {
915  x [Vi [p]] -= Vx [p] * tau ;
916  }
917  return (1) ;
918 }
919 /* create a Householder reflection [v,beta,s]=house(x), overwrite x with v,
920  * where (I-beta*v*v')*x = s*e1. See Algo 5.1.1, Golub & Van Loan, 3rd ed. */
921 double cs_house (double *x, double *beta, csi n)
922 {
923  double s, sigma = 0 ;
924  csi i ;
925  if (!x || !beta) return (-1) ; /* check inputs */
926  for (i = 1 ; i < n ; i++) sigma += x [i] * x [i] ;
927  if (sigma == 0)
928  {
929  s = fabs (x [0]) ; /* s = |x(0)| */
930  (*beta) = (x [0] <= 0) ? 2 : 0 ;
931  x [0] = 1 ;
932  }
933  else
934  {
935  s = sqrt (x [0] * x [0] + sigma) ; /* s = norm (x) */
936  x [0] = (x [0] <= 0) ? (x [0] - s) : (-sigma / (x [0] + s)) ;
937  (*beta) = -1. / (s * x [0]) ;
938  }
939  return (s) ;
940 }
941 /* x(p) = b, for dense vectors x and b; p=NULL denotes identity */
942 csi cs_ipvec (const csi *p, const double *b, double *x, csi n)
943 {
944  csi k ;
945  if (!x || !b) return (0) ; /* check inputs */
946  for (k = 0 ; k < n ; k++) x [p ? p [k] : k] = b [k] ;
947  return (1) ;
948 }
949 /* consider A(i,j), node j in ith row subtree and return lca(jprev,j) */
950 csi cs_leaf (csi i, csi j, const csi *first, csi *maxfirst, csi *prevleaf,
951  csi *ancestor, csi *jleaf)
952 {
953  csi q, s, sparent, jprev ;
954  if (!first || !maxfirst || !prevleaf || !ancestor || !jleaf) return (-1) ;
955  *jleaf = 0 ;
956  if (i <= j || first [j] <= maxfirst [i]) return (-1) ; /* j not a leaf */
957  maxfirst [i] = first [j] ; /* update max first[j] seen so far */
958  jprev = prevleaf [i] ; /* jprev = previous leaf of ith subtree */
959  prevleaf [i] = j ;
960  *jleaf = (jprev == -1) ? 1: 2 ; /* j is first or subsequent leaf */
961  if (*jleaf == 1) return (i) ; /* if 1st leaf, q = root of ith subtree */
962  for (q = jprev ; q != ancestor [q] ; q = ancestor [q]) ;
963  for (s = jprev ; s != q ; s = sparent)
964  {
965  sparent = ancestor [s] ; /* path compression */
966  ancestor [s] = q ;
967  }
968  return (q) ; /* q = least common ancester (jprev,j) */
969 }
970 /* load a triplet matrix from a file */
972 {
973  double i, j ; /* use double for integers to avoid csi conflicts */
974  double x ;
975  cs *T ;
976  if (!f) return (NULL) ; /* check inputs */
977  T = cs_spalloc (0, 0, 1, 1, 1) ; /* allocate result */
978  while (fscanf (f, "%lg %lg %lg\n", &i, &j, &x) == 3)
979  {
980  if (!cs_entry (T, (csi) i, (csi) j, x)) return (cs_spfree (T)) ;
981  }
982  return (T) ;
983 }
984 /* solve Lx=b where x and b are dense. x=b on input, solution on output. */
985 csi cs_lsolve (const cs *L, double *x)
986 {
987  csi p, j, n, *Lp, *Li ;
988  double *Lx ;
989  if (!CS_CSC (L) || !x) return (0) ; /* check inputs */
990  n = L->n ; Lp = L->p ; Li = L->i ; Lx = L->x ;
991  for (j = 0 ; j < n ; j++)
992  {
993  x [j] /= Lx [Lp [j]] ;
994  for (p = Lp [j]+1 ; p < Lp [j+1] ; p++)
995  {
996  x [Li [p]] -= Lx [p] * x [j] ;
997  }
998  }
999  return (1) ;
1000 }
1001 /* solve L'x=b where x and b are dense. x=b on input, solution on output. */
1002 csi cs_ltsolve (const cs *L, double *x)
1003 {
1004  csi p, j, n, *Lp, *Li ;
1005  double *Lx ;
1006  if (!CS_CSC (L) || !x) return (0) ; /* check inputs */
1007  n = L->n ; Lp = L->p ; Li = L->i ; Lx = L->x ;
1008  for (j = n-1 ; j >= 0 ; j--)
1009  {
1010  for (p = Lp [j]+1 ; p < Lp [j+1] ; p++)
1011  {
1012  x [j] -= Lx [p] * x [Li [p]] ;
1013  }
1014  x [j] /= Lx [Lp [j]] ;
1015  }
1016  return (1) ;
1017 }
1018 /* [L,U,pinv]=lu(A, [q lnz unz]). lnz and unz can be guess */
1019 csn *cs_lu (const cs *A, const css *S, double tol)
1020 {
1021  cs *L, *U ;
1022  csn *N ;
1023  double pivot, *Lx, *Ux, *x, a, t ;
1024  csi *Lp, *Li, *Up, *Ui, *pinv, *xi, *q, n, ipiv, k, top, p, i, col, lnz,unz;
1025  if (!CS_CSC (A) || !S) return (NULL) ; /* check inputs */
1026  n = A->n ;
1027  q = S->q ; lnz = S->lnz ; unz = S->unz ;
1028  x = cs_malloc (n, sizeof (double)) ; /* get double workspace */
1029  xi = cs_malloc (2*n, sizeof (csi)) ; /* get csi workspace */
1030  N = cs_calloc (1, sizeof (csn)) ; /* allocate result */
1031  if (!x || !xi || !N) return (cs_ndone (N, NULL, xi, x, 0)) ;
1032  N->L = L = cs_spalloc (n, n, lnz, 1, 0) ; /* allocate result L */
1033  N->U = U = cs_spalloc (n, n, unz, 1, 0) ; /* allocate result U */
1034  N->pinv = pinv = cs_malloc (n, sizeof (csi)) ; /* allocate result pinv */
1035  if (!L || !U || !pinv) return (cs_ndone (N, NULL, xi, x, 0)) ;
1036  Lp = L->p ; Up = U->p ;
1037  for (i = 0 ; i < n ; i++) x [i] = 0 ; /* clear workspace */
1038  for (i = 0 ; i < n ; i++) pinv [i] = -1 ; /* no rows pivotal yet */
1039  for (k = 0 ; k <= n ; k++) Lp [k] = 0 ; /* no cols of L yet */
1040  lnz = unz = 0 ;
1041  for (k = 0 ; k < n ; k++) /* compute L(:,k) and U(:,k) */
1042  {
1043  /* --- Triangular solve --------------------------------------------- */
1044  Lp [k] = lnz ; /* L(:,k) starts here */
1045  Up [k] = unz ; /* U(:,k) starts here */
1046  if ((lnz + n > L->nzmax && !cs_sprealloc (L, 2*L->nzmax + n)) ||
1047  (unz + n > U->nzmax && !cs_sprealloc (U, 2*U->nzmax + n)))
1048  {
1049  return (cs_ndone (N, NULL, xi, x, 0)) ;
1050  }
1051  Li = L->i ; Lx = L->x ; Ui = U->i ; Ux = U->x ;
1052  col = q ? (q [k]) : k ;
1053  top = cs_spsolve (L, A, col, xi, x, pinv, 1) ; /* x = L\A(:,col) */
1054  /* --- Find pivot --------------------------------------------------- */
1055  ipiv = -1 ;
1056  a = -1 ;
1057  for (p = top ; p < n ; p++)
1058  {
1059  i = xi [p] ; /* x(i) is nonzero */
1060  if (pinv [i] < 0) /* row i is not yet pivotal */
1061  {
1062  if ((t = fabs (x [i])) > a)
1063  {
1064  a = t ; /* largest pivot candidate so far */
1065  ipiv = i ;
1066  }
1067  }
1068  else /* x(i) is the entry U(pinv[i],k) */
1069  {
1070  Ui [unz] = pinv [i] ;
1071  Ux [unz++] = x [i] ;
1072  }
1073  }
1074  if (ipiv == -1 || a <= 0) return (cs_ndone (N, NULL, xi, x, 0)) ;
1075  if (pinv [col] < 0 && fabs (x [col]) >= a*tol) ipiv = col ;
1076  /* --- Divide by pivot ---------------------------------------------- */
1077  pivot = x [ipiv] ; /* the chosen pivot */
1078  Ui [unz] = k ; /* last entry in U(:,k) is U(k,k) */
1079  Ux [unz++] = pivot ;
1080  pinv [ipiv] = k ; /* ipiv is the kth pivot row */
1081  Li [lnz] = ipiv ; /* first entry in L(:,k) is L(k,k) = 1 */
1082  Lx [lnz++] = 1 ;
1083  for (p = top ; p < n ; p++) /* L(k+1:n,k) = x / pivot */
1084  {
1085  i = xi [p] ;
1086  if (pinv [i] < 0) /* x(i) is an entry in L(:,k) */
1087  {
1088  Li [lnz] = i ; /* save unpermuted row in L */
1089  Lx [lnz++] = x [i] / pivot ; /* scale pivot column */
1090  }
1091  x [i] = 0 ; /* x [0..n-1] = 0 for next k */
1092  }
1093  }
1094  /* --- Finalize L and U ------------------------------------------------- */
1095  Lp [n] = lnz ;
1096  Up [n] = unz ;
1097  Li = L->i ; /* fix row indices of L for final pinv */
1098  for (p = 0 ; p < lnz ; p++) Li [p] = pinv [Li [p]] ;
1099  cs_sprealloc (L, 0) ; /* remove extra space from L and U */
1100  cs_sprealloc (U, 0) ;
1101  return (cs_ndone (N, NULL, xi, x, 1)) ; /* success */
1102 }
1103 /* x=A\b where A is unsymmetric; b overwritten with solution */
1104 csi cs_lusol (csi order, const cs *A, double *b, double tol)
1105 {
1106  double *x ;
1107  css *S ;
1108  csn *N ;
1109  csi n, ok ;
1110  if (!CS_CSC (A) || !b) return (0) ; /* check inputs */
1111  n = A->n ;
1112  S = cs_sqr (order, A, 0) ; /* ordering and symbolic analysis */
1113  N = cs_lu (A, S, tol) ; /* numeric LU factorization */
1114  x = cs_malloc (n, sizeof (double)) ; /* get workspace */
1115  ok = (S && N && x) ;
1116  if (ok)
1117  {
1118  cs_ipvec (N->pinv, b, x, n) ; /* x = b(p) */
1119  cs_lsolve (N->L, x) ; /* x = L\x */
1120  cs_usolve (N->U, x) ; /* x = U\x */
1121  cs_ipvec (S->q, x, b, n) ; /* b(q) = x */
1122  }
1123  cs_free (x) ;
1124  cs_sfree (S) ;
1125  cs_nfree (N) ;
1126  return (ok) ;
1127 }
1128 #ifdef MATLAB_MEX_FILE
1129 #define malloc mxMalloc
1130 #define free mxFree
1131 #define realloc mxRealloc
1132 #define calloc mxCalloc
1133 #endif
1134
1135 /* wrapper for malloc */
1136 void *cs_malloc (csi n, size_t size)
1137 {
1138  return (malloc (CS_MAX (n,1) * size)) ;
1139 }
1140
1141 /* wrapper for calloc */
1142 void *cs_calloc (csi n, size_t size)
1143 {
1144  return (calloc (CS_MAX (n,1), size)) ;
1145 }
1146
1148 void *cs_free (void *p)
1149 {
1150  if (p) free (p) ; /* free p if it is not already NULL */
1151  return (NULL) ; /* return NULL to simplify the use of cs_free */
1152 }
1153
1154 /* wrapper for realloc */
1155 void *cs_realloc (void *p, csi n, size_t size, csi *ok)
1156 {
1157  void *pnew ;
1158  pnew = realloc (p, CS_MAX (n,1) * size) ; /* realloc the block */
1159  *ok = (pnew != NULL) ; /* realloc fails if pnew is NULL */
1160  return ((*ok) ? pnew : p) ; /* return original p if failure */
1161 }
1162 /* find an augmenting path starting at column k and extend the match if found */
1163 static void cs_augment (csi k, const cs *A, csi *jmatch, csi *cheap, csi *w,
1164  csi *js, csi *is, csi *ps)
1165 {
1166  csi found = 0, p, i = -1, *Ap = A->p, *Ai = A->i, head = 0, j ;
1167  js [0] = k ; /* start with just node k in jstack */
1169  {
1170  /* --- Start (or continue) depth-first-search at node j ------------- */
1171  j = js [head] ; /* get j from top of jstack */
1172  if (w [j] != k) /* 1st time j visited for kth path */
1173  {
1174  w [j] = k ; /* mark j as visited for kth path */
1175  for (p = cheap [j] ; p < Ap [j+1] && !found ; p++)
1176  {
1177  i = Ai [p] ; /* try a cheap assignment (i,j) */
1178  found = (jmatch [i] == -1) ;
1179  }
1180  cheap [j] = p ; /* start here next time j is traversed*/
1181  if (found)
1182  {
1183  is [head] = i ; /* column j matched with row i */
1184  break ; /* end of augmenting path */
1185  }
1186  ps [head] = Ap [j] ; /* no cheap match: start dfs for j */
1187  }
1188  /* --- Depth-first-search of neighbors of j ------------------------- */
1189  for (p = ps [head] ; p < Ap [j+1] ; p++)
1190  {
1191  i = Ai [p] ; /* consider row i */
1192  if (w [jmatch [i]] == k) continue ; /* skip jmatch [i] if marked */
1193  ps [head] = p + 1 ; /* pause dfs of node j */
1194  is [head] = i ; /* i will be matched with j if found */
1195  js [++head] = jmatch [i] ; /* start dfs at column jmatch [i] */
1196  break ;
1197  }
1198  if (p == Ap [j+1]) head-- ; /* node j is done; pop from stack */
1199  } /* augment the match if path found: */
1200  if (found) for (p = head ; p >= 0 ; p--) jmatch [is [p]] = js [p] ;
1201 }
1202
1203 /* find a maximum transveral */
1204 csi *cs_maxtrans (const cs *A, csi seed) /*[jmatch [0..m-1]; imatch [0..n-1]]*/
1205 {
1206  csi i, j, k, n, m, p, n2 = 0, m2 = 0, *Ap, *jimatch, *w, *cheap, *js, *is,
1207  *ps, *Ai, *Cp, *jmatch, *imatch, *q ;
1208  cs *C ;
1209  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
1210  n = A->n ; m = A->m ; Ap = A->p ; Ai = A->i ;
1211  w = jimatch = cs_calloc (m+n, sizeof (csi)) ; /* allocate result */
1212  if (!jimatch) return (NULL) ;
1213  for (k = 0, j = 0 ; j < n ; j++) /* count nonempty rows and columns */
1214  {
1215  n2 += (Ap [j] < Ap [j+1]) ;
1216  for (p = Ap [j] ; p < Ap [j+1] ; p++)
1217  {
1218  w [Ai [p]] = 1 ;
1219  k += (j == Ai [p]) ; /* count entries already on diagonal */
1220  }
1221  }
1222  if (k == CS_MIN (m,n)) /* quick return if diagonal zero-free */
1223  {
1224  jmatch = jimatch ; imatch = jimatch + m ;
1225  for (i = 0 ; i < k ; i++) jmatch [i] = i ;
1226  for ( ; i < m ; i++) jmatch [i] = -1 ;
1227  for (j = 0 ; j < k ; j++) imatch [j] = j ;
1228  for ( ; j < n ; j++) imatch [j] = -1 ;
1229  return (cs_idone (jimatch, NULL, NULL, 1)) ;
1230  }
1231  for (i = 0 ; i < m ; i++) m2 += w [i] ;
1232  C = (m2 < n2) ? cs_transpose (A,0) : ((cs *) A) ; /* transpose if needed */
1233  if (!C) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, NULL, 0)) ;
1234  n = C->n ; m = C->m ; Cp = C->p ;
1235  jmatch = (m2 < n2) ? jimatch + n : jimatch ;
1236  imatch = (m2 < n2) ? jimatch : jimatch + m ;
1237  w = cs_malloc (5*n, sizeof (csi)) ; /* get workspace */
1238  if (!w) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 0)) ;
1239  cheap = w + n ; js = w + 2*n ; is = w + 3*n ; ps = w + 4*n ;
1240  for (j = 0 ; j < n ; j++) cheap [j] = Cp [j] ; /* for cheap assignment */
1241  for (j = 0 ; j < n ; j++) w [j] = -1 ; /* all columns unflagged */
1242  for (i = 0 ; i < m ; i++) jmatch [i] = -1 ; /* nothing matched yet */
1243  q = cs_randperm (n, seed) ; /* q = random permutation */
1244  for (k = 0 ; k < n ; k++) /* augment, starting at column q[k] */
1245  {
1246  cs_augment (q ? q [k]: k, C, jmatch, cheap, w, js, is, ps) ;
1247  }
1248  cs_free (q) ;
1249  for (j = 0 ; j < n ; j++) imatch [j] = -1 ; /* find row match */
1250  for (i = 0 ; i < m ; i++) if (jmatch [i] >= 0) imatch [jmatch [i]] = i ;
1251  return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 1)) ;
1252 }
1253 /* C = A*B */
1254 cs *cs_multiply (const cs *A, const cs *B)
1255 {
1256  csi p, j, nz = 0, anz, *Cp, *Ci, *Bp, m, n, bnz, *w, values, *Bi ;
1257  double *x, *Bx, *Cx ;
1258  cs *C ;
1259  if (!CS_CSC (A) || !CS_CSC (B)) return (NULL) ; /* check inputs */
1260  if (A->n != B->m) return (NULL) ;
1261  m = A->m ; anz = A->p [A->n] ;
1262  n = B->n ; Bp = B->p ; Bi = B->i ; Bx = B->x ; bnz = Bp [n] ;
1263  w = cs_calloc (m, sizeof (csi)) ; /* get workspace */
1264  values = (A->x != NULL) && (Bx != NULL) ;
1265  x = values ? cs_malloc (m, sizeof (double)) : NULL ; /* get workspace */
1266  C = cs_spalloc (m, n, anz + bnz, values, 0) ; /* allocate result */
1267  if (!C || !w || (values && !x)) return (cs_done (C, w, x, 0)) ;
1268  Cp = C->p ;
1269  for (j = 0 ; j < n ; j++)
1270  {
1271  if (nz + m > C->nzmax && !cs_sprealloc (C, 2*(C->nzmax)+m))
1272  {
1273  return (cs_done (C, w, x, 0)) ; /* out of memory */
1274  }
1275  Ci = C->i ; Cx = C->x ; /* C->i and C->x may be reallocated */
1276  Cp [j] = nz ; /* column j of C starts here */
1277  for (p = Bp [j] ; p < Bp [j+1] ; p++)
1278  {
1279  nz = cs_scatter (A, Bi [p], Bx ? Bx [p] : 1, w, x, j+1, C, nz) ;
1280  }
1281  if (values) for (p = Cp [j] ; p < nz ; p++) Cx [p] = x [Ci [p]] ;
1282  }
1283  Cp [n] = nz ; /* finalize the last column of C */
1284  cs_sprealloc (C, 0) ; /* remove extra space from C */
1285  return (cs_done (C, w, x, 1)) ; /* success; free workspace, return C */
1286 }
1287 /* 1-norm of a sparse matrix = max (sum (abs (A))), largest column sum */
1288 double cs_norm (const cs *A)
1289 {
1290  csi p, j, n, *Ap ;
1291  double *Ax, norm = 0, s ;
1292  if (!CS_CSC (A) || !A->x) return (-1) ; /* check inputs */
1293  n = A->n ; Ap = A->p ; Ax = A->x ;
1294  for (j = 0 ; j < n ; j++)
1295  {
1296  for (s = 0, p = Ap [j] ; p < Ap [j+1] ; p++) s += fabs (Ax [p]) ;
1297  norm = CS_MAX (norm, s) ;
1298  }
1299  return (norm) ;
1300 }
1301 /* C = A(p,q) where p and q are permutations of 0..m-1 and 0..n-1. */
1302 cs *cs_permute (const cs *A, const csi *pinv, const csi *q, csi values)
1303 {
1304  csi t, j, k, nz = 0, m, n, *Ap, *Ai, *Cp, *Ci ;
1305  double *Cx, *Ax ;
1306  cs *C ;
1307  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
1308  m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
1309  C = cs_spalloc (m, n, Ap [n], values && Ax != NULL, 0) ; /* alloc result */
1310  if (!C) return (cs_done (C, NULL, NULL, 0)) ; /* out of memory */
1311  Cp = C->p ; Ci = C->i ; Cx = C->x ;
1312  for (k = 0 ; k < n ; k++)
1313  {
1314  Cp [k] = nz ; /* column k of C is column q[k] of A */
1315  j = q ? (q [k]) : k ;
1316  for (t = Ap [j] ; t < Ap [j+1] ; t++)
1317  {
1318  if (Cx) Cx [nz] = Ax [t] ; /* row i of A is row pinv[i] of C */
1319  Ci [nz++] = pinv ? (pinv [Ai [t]]) : Ai [t] ;
1320  }
1321  }
1322  Cp [n] = nz ; /* finalize the last column of C */
1323  return (cs_done (C, NULL, NULL, 1)) ;
1324 }
1325 /* pinv = p', or p = pinv' */
1326 csi *cs_pinv (csi const *p, csi n)
1327 {
1328  csi k, *pinv ;
1329  if (!p) return (NULL) ; /* p = NULL denotes identity */
1330  pinv = cs_malloc (n, sizeof (csi)) ; /* allocate result */
1331  if (!pinv) return (NULL) ; /* out of memory */
1332  for (k = 0 ; k < n ; k++) pinv [p [k]] = k ;/* invert the permutation */
1333  return (pinv) ; /* return result */
1334 }
1335 /* post order a forest */
1336 csi *cs_post (const csi *parent, csi n)
1337 {
1338  csi j, k = 0, *post, *w, *head, *next, *stack ;
1339  if (!parent) return (NULL) ; /* check inputs */
1340  post = cs_malloc (n, sizeof (csi)) ; /* allocate result */
1341  w = cs_malloc (3*n, sizeof (csi)) ; /* get workspace */
1342  if (!w || !post) return (cs_idone (post, NULL, w, 0)) ;
1343  head = w ; next = w + n ; stack = w + 2*n ;
1344  for (j = 0 ; j < n ; j++) head [j] = -1 ; /* empty linked lists */
1345  for (j = n-1 ; j >= 0 ; j--) /* traverse nodes in reverse order*/
1346  {
1347  if (parent [j] == -1) continue ; /* j is a root */
1348  next [j] = head [parent [j]] ; /* add j to list of its parent */
1349  head [parent [j]] = j ;
1350  }
1351  for (j = 0 ; j < n ; j++)
1352  {
1353  if (parent [j] != -1) continue ; /* skip j if it is not a root */
1354  k = cs_tdfs (j, k, head, next, post, stack) ;
1355  }
1356  return (cs_idone (post, NULL, w, 1)) ; /* success; free w, return post */
1357 }
1358 /* print a sparse matrix; use %g for integers to avoid differences with csi */
1359 csi cs_print (const cs *A, csi brief)
1360 {
1361  csi p, j, m, n, nzmax, nz, *Ap, *Ai ;
1362  double *Ax ;
1363  if (!A) { Rprintf ("(null)\n") ; return (0) ; }
1364  m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
1365  nzmax = A->nzmax ; nz = A->nz ;
1366  Rprintf ("CSparse Version %d.%d.%d, %s. %s\n", CS_VER, CS_SUBVER,
1368  if (nz < 0)
1369  {
1370  Rprintf ("%g-by-%g, nzmax: %g nnz: %g, 1-norm: %g\n", (double) m,
1371  (double) n, (double) nzmax, (double) (Ap [n]), cs_norm (A)) ;
1372  for (j = 0 ; j < n ; j++)
1373  {
1374  Rprintf (" col %g : locations %g to %g\n", (double) j,
1375  (double) (Ap [j]), (double) (Ap [j+1]-1)) ;
1376  for (p = Ap [j] ; p < Ap [j+1] ; p++)
1377  {
1378  Rprintf (" %g : %g\n", (double) (Ai [p]), Ax ? Ax [p] : 1) ;
1379  if (brief && p > 20) { Rprintf (" ...\n") ; return (1) ; }
1380  }
1381  }
1382  }
1383  else
1384  {
1385  Rprintf ("triplet: %g-by-%g, nzmax: %g nnz: %g\n", (double) m,
1386  (double) n, (double) nzmax, (double) nz) ;
1387  for (p = 0 ; p < nz ; p++)
1388  {
1389  Rprintf (" %g %g : %g\n", (double) (Ai [p]), (double) (Ap [p]),
1390  Ax ? Ax [p] : 1) ;
1391  if (brief && p > 20) { Rprintf (" ...\n") ; return (1) ; }
1392  }
1393  }
1394  return (1) ;
1395 }
1396 /* x = b(p), for dense vectors x and b; p=NULL denotes identity */
1397 csi cs_pvec (const csi *p, const double *b, double *x, csi n)
1398 {
1399  csi k ;
1400  if (!x || !b) return (0) ; /* check inputs */
1401  for (k = 0 ; k < n ; k++) x [k] = b [p ? p [k] : k] ;
1402  return (1) ;
1403 }
1404 /* sparse QR factorization [V,beta,pinv,R] = qr (A) */
1405 csn *cs_qr (const cs *A, const css *S)
1406 {
1407  double *Rx, *Vx, *Ax, *x, *Beta ;
1408  csi i, k, p, m, n, vnz, p1, top, m2, len, col, rnz, *s, *leftmost, *Ap, *Ai,
1409  *parent, *Rp, *Ri, *Vp, *Vi, *w, *pinv, *q ;
1410  cs *R, *V ;
1411  csn *N ; // the result
1412  if (!CS_CSC (A) || !S) return (NULL) ;
1413  m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
1414  q = S->q ; parent = S->parent ; pinv = S->pinv ; m2 = S->m2 ;
1415  vnz = S->lnz ; rnz = S->unz ; leftmost = S->leftmost ;
1416  w = cs_malloc (m2+n, sizeof (csi)) ; /* get csi workspace */
1417  x = cs_malloc (m2, sizeof (double)) ; /* get double workspace */
1418  N = cs_calloc (1, sizeof (csn)) ; /* allocate result */
1419  if (!w || !x || !N) return (cs_ndone (N, NULL, w, x, 0)) ;
1420  s = w + m2 ; /* s is size n */
1421  for (k = 0 ; k < m2 ; k++) x [k] = 0 ; /* clear workspace x */
1422  N->L = V = cs_spalloc (m2, n, vnz, 1, 0) ; /* allocate result V */
1423  N->U = R = cs_spalloc (m2, n, rnz, 1, 0) ; /* allocate result R */
1424  N->B = Beta = cs_malloc (n, sizeof (double)) ; /* allocate result Beta */
1425  if (!R || !V || !Beta) return (cs_ndone (N, NULL, w, x, 0)) ;
1426  Rp = R->p ; Ri = R->i ; Rx = R->x ;
1427  Vp = V->p ; Vi = V->i ; Vx = V->x ;
1428  for (i = 0 ; i < m2 ; i++) w [i] = -1 ; /* clear w, to mark nodes */
1429  rnz = 0 ; vnz = 0 ;
1430  for (k = 0 ; k < n ; k++) /* compute V and R */
1431  {
1432  Rp [k] = rnz ; /* R(:,k) starts here */
1433  Vp [k] = p1 = vnz ; /* V(:,k) starts here */
1434  w [k] = k ; /* add V(k,k) to pattern of V */
1435  Vi [vnz++] = k ;
1436  top = n ;
1437  col = q ? q [k] : k ;
1438  for (p = Ap [col] ; p < Ap [col+1] ; p++) /* find R(:,k) pattern */
1439  {
1440  i = leftmost [Ai [p]] ; /* i = min(find(A(i,q))) */
1441  for (len = 0 ; w [i] != k ; i = parent [i]) /* traverse up to k */
1442  {
1443  s [len++] = i ;
1444  w [i] = k ;
1445  }
1446  while (len > 0) s [--top] = s [--len] ; /* push path on stack */
1447  i = pinv [Ai [p]] ; /* i = permuted row of A(:,col) */
1448  x [i] = Ax [p] ; /* x (i) = A(:,col) */
1449  if (i > k && w [i] < k) /* pattern of V(:,k) = x (k+1:m) */
1450  {
1451  Vi [vnz++] = i ; /* add i to pattern of V(:,k) */
1452  w [i] = k ;
1453  }
1454  }
1455  for (p = top ; p < n ; p++) /* for each i in pattern of R(:,k) */
1456  {
1457  i = s [p] ; /* R(i,k) is nonzero */
1458  cs_happly (V, i, Beta [i], x) ; /* apply (V(i),Beta(i)) to x */
1459  Ri [rnz] = i ; /* R(i,k) = x(i) */
1460  Rx [rnz++] = x [i] ;
1461  x [i] = 0 ;
1462  if (parent [i] == k) vnz = cs_scatter (V, i, 0, w, NULL, k, V, vnz);
1463  }
1464  for (p = p1 ; p < vnz ; p++) /* gather V(:,k) = x */
1465  {
1466  Vx [p] = x [Vi [p]] ;
1467  x [Vi [p]] = 0 ;
1468  }
1469  Ri [rnz] = k ; /* R(k,k) = norm (x) */
1470  Rx [rnz++] = cs_house (Vx+p1, Beta+k, vnz-p1) ; /* [v,beta]=house(x) */
1471  }
1472  Rp [n] = rnz ; /* finalize R */
1473  Vp [n] = vnz ; /* finalize V */
1474  return (cs_ndone (N, NULL, w, x, 1)) ; /* success */
1475 }
1476 /* x=A\b where A can be rectangular; b overwritten with solution */
1477 csi cs_qrsol (csi order, const cs *A, double *b)
1478 {
1479  double *x ;
1480  css *S ;
1481  csn *N ;
1482  cs *AT = NULL ;
1483  csi k, m, n, ok ;
1484  if (!CS_CSC (A) || !b) return (0) ; /* check inputs */
1485  n = A->n ;
1486  m = A->m ;
1487  if (m >= n)
1488  {
1489  S = cs_sqr (order, A, 1) ; /* ordering and symbolic analysis */
1490  N = cs_qr (A, S) ; /* numeric QR factorization */
1491  x = cs_calloc (S ? S->m2 : 1, sizeof (double)) ; /* get workspace */
1492  ok = (S && N && x) ;
1493  if (ok)
1494  {
1495  cs_ipvec (S->pinv, b, x, m) ; /* x(0:m-1) = b(p(0:m-1) */
1496  for (k = 0 ; k < n ; k++) /* apply Householder refl. to x */
1497  {
1498  cs_happly (N->L, k, N->B [k], x) ;
1499  }
1500  cs_usolve (N->U, x) ; /* x = R\x */
1501  cs_ipvec (S->q, x, b, n) ; /* b(q(0:n-1)) = x(0:n-1) */
1502  }
1503  }
1504  else
1505  {
1506  AT = cs_transpose (A, 1) ; /* Ax=b is underdetermined */
1507  S = cs_sqr (order, AT, 1) ; /* ordering and symbolic analysis */
1508  N = cs_qr (AT, S) ; /* numeric QR factorization of A' */
1509  x = cs_calloc (S ? S->m2 : 1, sizeof (double)) ; /* get workspace */
1510  ok = (AT && S && N && x) ;
1511  if (ok)
1512  {
1513  cs_pvec (S->q, b, x, m) ; /* x(q(0:m-1)) = b(0:m-1) */
1514  cs_utsolve (N->U, x) ; /* x = R'\x */
1515  for (k = m-1 ; k >= 0 ; k--) /* apply Householder refl. to x */
1516  {
1517  cs_happly (N->L, k, N->B [k], x) ;
1518  }
1519  cs_pvec (S->pinv, x, b, n) ; /* b(0:n-1) = x(p(0:n-1)) */
1520  }
1521  }
1522  cs_free (x) ;
1523  cs_sfree (S) ;
1524  cs_nfree (N) ;
1525  cs_spfree (AT) ;
1526  return (ok) ;
1527 }
1528 /* return a random permutation vector, the identity perm, or p = n-1:-1:0.
1529  * seed = -1 means p = n-1:-1:0. seed = 0 means p = identity. otherwise
1530  * p = random permutation. */
1532 {
1533  csi *p, k, j, t ;
1534  if (seed == 0) return (NULL) ; /* return p = NULL (identity) */
1535  p = cs_malloc (n, sizeof (csi)) ; /* allocate result */
1536  if (!p) return (NULL) ; /* out of memory */
1537  for (k = 0 ; k < n ; k++) p [k] = n-k-1 ;
1538  if (seed == -1) return (p) ; /* return reverse permutation */
1539  GetRNGstate();/* <- for R package Matrix
1540  srand (seed) ; .* get new random number seed */
1541  for (k = 0 ; k < n ; k++)
1542  {
1543  j = k + (int)(unif_rand() * (n-k)); // j = rand integer in range k to n-1
1544  t = p [j] ; /* swap p[k] and p[j] */
1545  p [j] = p [k] ;
1546  p [k] = t ;
1547  }
1548  PutRNGstate(); // <- R package Matrix
1549  return (p) ;
1550 }
1551 /* xi [top...n-1] = nodes reachable from graph of G*P' via nodes in B(:,k).
1552  * xi [n...2n-1] used as workspace */
1553 csi cs_reach (cs *G, const cs *B, csi k, csi *xi, const csi *pinv)
1554 {
1555  csi p, n, top, *Bp, *Bi, *Gp ;
1556  if (!CS_CSC (G) || !CS_CSC (B) || !xi) return (-1) ; /* check inputs */
1557  n = G->n ; Bp = B->p ; Bi = B->i ; Gp = G->p ;
1558  top = n ;
1559  for (p = Bp [k] ; p < Bp [k+1] ; p++)
1560  {
1561  if (!CS_MARKED (Gp, Bi [p])) /* start a dfs at unmarked node i */
1562  {
1563  top = cs_dfs (Bi [p], G, top, xi, xi+n, pinv) ;
1564  }
1565  }
1566  for (p = top ; p < n ; p++) CS_MARK (Gp, xi [p]) ; /* restore G */
1567  return (top) ;
1568 }
1569 /* x = x + beta * A(:,j), where x is a dense vector and A(:,j) is sparse */
1570 csi cs_scatter (const cs *A, csi j, double beta, csi *w, double *x, csi mark,
1571  cs *C, csi nz)
1572 {
1573  csi i, p, *Ap, *Ai, *Ci ;
1574  double *Ax ;
1575  if (!CS_CSC (A) || !w || !CS_CSC (C)) return (-1) ; /* check inputs */
1576  Ap = A->p ; Ai = A->i ; Ax = A->x ; Ci = C->i ;
1577  for (p = Ap [j] ; p < Ap [j+1] ; p++)
1578  {
1579  i = Ai [p] ; /* A(i,j) is nonzero */
1580  if (w [i] < mark)
1581  {
1582  w [i] = mark ; /* i is new entry in column j */
1583  Ci [nz++] = i ; /* add i to pattern of C(:,j) */
1584  if (x) x [i] = beta * Ax [p] ; /* x(i) = beta*A(i,j) */
1585  }
1586  else if (x) x [i] += beta * Ax [p] ; /* i exists in C(:,j) already */
1587  }
1588  return (nz) ;
1589 }
1590 /* find the strongly connected components of a square matrix */
1591 csd *cs_scc (cs *A) /* matrix A temporarily modified, then restored */
1592 {
1593  csi n, i, k, b, nb = 0, top, *xi, *pstack, *p, *r, *Ap, *ATp, *rcopy, *Blk ;
1594  cs *AT ;
1595  csd *D ;
1596  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
1597  n = A->n ; Ap = A->p ;
1598  D = cs_dalloc (n, 0) ; /* allocate result */
1599  AT = cs_transpose (A, 0) ; /* AT = A' */
1600  xi = cs_malloc (2*n+1, sizeof (csi)) ; /* get workspace */
1601  if (!D || !AT || !xi) return (cs_ddone (D, AT, xi, 0)) ;
1602  Blk = xi ; rcopy = pstack = xi + n ;
1603  p = D->p ; r = D->r ; ATp = AT->p ;
1604  top = n ;
1605  for (i = 0 ; i < n ; i++) /* first dfs(A) to find finish times (xi) */
1606  {
1607  if (!CS_MARKED (Ap, i)) top = cs_dfs (i, A, top, xi, pstack, NULL) ;
1608  }
1609  for (i = 0 ; i < n ; i++) CS_MARK (Ap, i) ; /* restore A; unmark all nodes*/
1610  top = n ;
1611  nb = n ;
1612  for (k = 0 ; k < n ; k++) /* dfs(A') to find strongly connnected comp */
1613  {
1614  i = xi [k] ; /* get i in reverse order of finish times */
1615  if (CS_MARKED (ATp, i)) continue ; /* skip node i if already ordered */
1616  r [nb--] = top ; /* node i is the start of a component in p */
1617  top = cs_dfs (i, AT, top, p, pstack, NULL) ;
1618  }
1619  r [nb] = 0 ; /* first block starts at zero; shift r up */
1620  for (k = nb ; k <= n ; k++) r [k-nb] = r [k] ;
1621  D->nb = nb = n-nb ; /* nb = # of strongly connected components */
1622  for (b = 0 ; b < nb ; b++) /* sort each block in natural order */
1623  {
1624  for (k = r [b] ; k < r [b+1] ; k++) Blk [p [k]] = b ;
1625  }
1626  for (b = 0 ; b <= nb ; b++) rcopy [b] = r [b] ;
1627  for (i = 0 ; i < n ; i++) p [rcopy [Blk [i]]++] = i ;
1628  return (cs_ddone (D, AT, xi, 1)) ;
1629 }
1630 /* ordering and symbolic analysis for a Cholesky factorization */
1631 css *cs_schol (csi order, const cs *A)
1632 {
1633  csi n, *c, *post, *P ;
1634  cs *C ;
1635  css *S ;
1636  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
1637  n = A->n ;
1638  S = cs_calloc (1, sizeof (css)) ; /* allocate result S */
1639  if (!S) return (NULL) ; /* out of memory */
1640  P = cs_amd (order, A) ; /* P = amd(A+A'), or natural */
1641  S->pinv = cs_pinv (P, n) ; /* find inverse permutation */
1642  cs_free (P) ;
1643  if (order && !S->pinv) return (cs_sfree (S)) ;
1644  C = cs_symperm (A, S->pinv, 0) ; /* C = spones(triu(A(P,P))) */
1645  S->parent = cs_etree (C, 0) ; /* find etree of C */
1646  post = cs_post (S->parent, n) ; /* postorder the etree */
1647  c = cs_counts (C, S->parent, post, 0) ; /* find column counts of chol(C) */
1648  cs_free (post) ;
1649  cs_spfree (C) ;
1650  S->cp = cs_malloc (n+1, sizeof (csi)) ; /* allocate result S->cp */
1651  S->unz = S->lnz = cs_cumsum (S->cp, c, n) ; /* find column pointers for L */
1652  cs_free (c) ;
1653  return ((S->lnz >= 0) ? S : cs_sfree (S)) ;
1654 }
1655 /* solve Gx=b(:,k), where G is either upper (lo=0) or lower (lo=1) triangular */
1656 csi cs_spsolve (cs *G, const cs *B, csi k, csi *xi, double *x, const csi *pinv,
1657  csi lo)
1658 {
1659  csi j, J, p, q, px, top, n, *Gp, *Gi, *Bp, *Bi ;
1660  double *Gx, *Bx ;
1661  if (!CS_CSC (G) || !CS_CSC (B) || !xi || !x) return (-1) ;
1662  Gp = G->p ; Gi = G->i ; Gx = G->x ; n = G->n ;
1663  Bp = B->p ; Bi = B->i ; Bx = B->x ;
1664  top = cs_reach (G, B, k, xi, pinv) ; /* xi[top..n-1]=Reach(B(:,k)) */
1665  for (p = top ; p < n ; p++) x [xi [p]] = 0 ; /* clear x */
1666  for (p = Bp [k] ; p < Bp [k+1] ; p++) x [Bi [p]] = Bx [p] ; /* scatter B */
1667  for (px = top ; px < n ; px++)
1668  {
1669  j = xi [px] ; /* x(j) is nonzero */
1670  J = pinv ? (pinv [j]) : j ; /* j maps to col J of G */
1671  if (J < 0) continue ; /* column J is empty */
1672  x [j] /= Gx [lo ? (Gp [J]) : (Gp [J+1]-1)] ;/* x(j) /= G(j,j) */
1673  p = lo ? (Gp [J]+1) : (Gp [J]) ; /* lo: L(j,j) 1st entry */
1674  q = lo ? (Gp [J+1]) : (Gp [J+1]-1) ; /* up: U(j,j) last entry */
1675  for ( ; p < q ; p++)
1676  {
1677  x [Gi [p]] -= Gx [p] * x [j] ; /* x(i) -= G(i,j) * x(j) */
1678  }
1679  }
1681 }
1682 /* compute nnz(V) = S->lnz, S->pinv, S->leftmost, S->m2 from A and S->parent */
1683 static csi cs_vcount (const cs *A, css *S)
1684 {
1685  csi i, k, p, pa, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i, *next, *head,
1686  *tail, *nque, *pinv, *leftmost, *w, *parent = S->parent ;
1687  S->pinv = pinv = cs_malloc (m+n, sizeof (csi)) ; /* allocate pinv, */
1688  S->leftmost = leftmost = cs_malloc (m, sizeof (csi)) ; /* and leftmost */
1689  w = cs_malloc (m+3*n, sizeof (csi)) ; /* get workspace */
1690  if (!pinv || !w || !leftmost)
1691  {
1692  cs_free (w) ; /* pinv and leftmost freed later */
1693  return (0) ; /* out of memory */
1694  }
1695  next = w ; head = w + m ; tail = w + m + n ; nque = w + m + 2*n ;
1696  for (k = 0 ; k < n ; k++) head [k] = -1 ; /* queue k is empty */
1697  for (k = 0 ; k < n ; k++) tail [k] = -1 ;
1698  for (k = 0 ; k < n ; k++) nque [k] = 0 ;
1699  for (i = 0 ; i < m ; i++) leftmost [i] = -1 ;
1700  for (k = n-1 ; k >= 0 ; k--)
1701  {
1702  for (p = Ap [k] ; p < Ap [k+1] ; p++)
1703  {
1704  leftmost [Ai [p]] = k ; /* leftmost[i] = min(find(A(i,:)))*/
1705  }
1706  }
1707  for (i = m-1 ; i >= 0 ; i--) /* scan rows in reverse order */
1708  {
1709  pinv [i] = -1 ; /* row i is not yet ordered */
1710  k = leftmost [i] ;
1711  if (k == -1) continue ; /* row i is empty */
1712  if (nque [k]++ == 0) tail [k] = i ; /* first row in queue k */
1713  next [i] = head [k] ; /* put i at head of queue k */
1714  head [k] = i ;
1715  }
1716  S->lnz = 0 ;
1717  S->m2 = m ;
1718  for (k = 0 ; k < n ; k++) /* find row permutation and nnz(V)*/
1719  {
1720  i = head [k] ; /* remove row i from queue k */
1721  S->lnz++ ; /* count V(k,k) as nonzero */
1722  if (i < 0) i = S->m2++ ; /* add a fictitious row */
1723  pinv [i] = k ; /* associate row i with V(:,k) */
1724  if (--nque [k] <= 0) continue ; /* skip if V(k+1:m,k) is empty */
1725  S->lnz += nque [k] ; /* nque [k] is nnz (V(k+1:m,k)) */
1726  if ((pa = parent [k]) != -1) /* move all rows to parent of k */
1727  {
1728  if (nque [pa] == 0) tail [pa] = tail [k] ;
1729  next [tail [k]] = head [pa] ;
1730  head [pa] = next [i] ;
1731  nque [pa] += nque [k] ;
1732  }
1733  }
1734  for (i = 0 ; i < m ; i++) if (pinv [i] < 0) pinv [i] = k++ ;
1735  cs_free (w) ;
1736  return (1) ;
1737 }
1738
1739 /* symbolic ordering and analysis for QR or LU */
1740 css *cs_sqr (csi order, const cs *A, csi qr)
1741 {
1742  csi n, k, ok = 1, *post ;
1743  css *S ;
1744  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
1745  n = A->n ;
1746  S = cs_calloc (1, sizeof (css)) ; /* allocate result S */
1747  if (!S) return (NULL) ; /* out of memory */
1748  S->q = cs_amd (order, A) ; /* fill-reducing ordering */
1749  if (order && !S->q) return (cs_sfree (S)) ;
1750  if (qr) /* QR symbolic analysis */
1751  {
1752  cs *C = order ? cs_permute (A, NULL, S->q, 0) : ((cs *) A) ;
1753  S->parent = cs_etree (C, 1) ; /* etree of C'*C, where C=A(:,q) */
1754  post = cs_post (S->parent, n) ;
1755  S->cp = cs_counts (C, S->parent, post, 1) ; /* col counts chol(C'*C) */
1756  cs_free (post) ;
1757  ok = C && S->parent && S->cp && cs_vcount (C, S) ;
1758  if (ok) for (S->unz = 0, k = 0 ; k < n ; k++) S->unz += S->cp [k] ;
1759  if (order) cs_spfree (C) ;
1760  }
1761  else
1762  {
1763  S->unz = 4*(A->p [n]) + n ; /* for LU factorization only, */
1764  S->lnz = S->unz ; /* guess nnz(L) and nnz(U) */
1765  }
1766  return (ok ? S : cs_sfree (S)) ; /* return result S */
1767 }
1768 /* C = A(p,p) where A and C are symmetric the upper part stored; pinv not p */
1769 cs *cs_symperm (const cs *A, const csi *pinv, csi values)
1770 {
1771  csi i, j, p, q, i2, j2, n, *Ap, *Ai, *Cp, *Ci, *w ;
1772  double *Cx, *Ax ;
1773  cs *C ;
1774  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
1775  n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
1776  C = cs_spalloc (n, n, Ap [n], values && (Ax != NULL), 0) ; /* alloc result*/
1777  w = cs_calloc (n, sizeof (csi)) ; /* get workspace */
1778  if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */
1779  Cp = C->p ; Ci = C->i ; Cx = C->x ;
1780  for (j = 0 ; j < n ; j++) /* count entries in each column of C */
1781  {
1782  j2 = pinv ? pinv [j] : j ; /* column j of A is column j2 of C */
1783  for (p = Ap [j] ; p < Ap [j+1] ; p++)
1784  {
1785  i = Ai [p] ;
1786  if (i > j) continue ; /* skip lower triangular part of A */
1787  i2 = pinv ? pinv [i] : i ; /* row i of A is row i2 of C */
1788  w [CS_MAX (i2, j2)]++ ; /* column count of C */
1789  }
1790  }
1791  cs_cumsum (Cp, w, n) ; /* compute column pointers of C */
1792  for (j = 0 ; j < n ; j++)
1793  {
1794  j2 = pinv ? pinv [j] : j ; /* column j of A is column j2 of C */
1795  for (p = Ap [j] ; p < Ap [j+1] ; p++)
1796  {
1797  i = Ai [p] ;
1798  if (i > j) continue ; /* skip lower triangular part of A*/
1799  i2 = pinv ? pinv [i] : i ; /* row i of A is row i2 of C */
1800  Ci [q = w [CS_MAX (i2, j2)]++] = CS_MIN (i2, j2) ;
1801  if (Cx) Cx [q] = Ax [p] ;
1802  }
1803  }
1804  return (cs_done (C, w, NULL, 1)) ; /* success; free workspace, return C */
1805 }
1806 /* depth-first search and postorder of a tree rooted at node j */
1807 csi cs_tdfs (csi j, csi k, csi *head, const csi *next, csi *post, csi *stack)
1808 {
1809  csi i, p, top = 0 ;
1810  if (!head || !next || !post || !stack) return (-1) ; /* check inputs */
1811  stack [0] = j ; /* place j on the stack */
1812  while (top >= 0) /* while (stack is not empty) */
1813  {
1814  p = stack [top] ; /* p = top of stack */
1815  i = head [p] ; /* i = youngest child of p */
1816  if (i == -1)
1817  {
1818  top-- ; /* p has no unordered children left */
1819  post [k++] = p ; /* node p is the kth postordered node */
1820  }
1821  else
1822  {
1823  head [p] = next [i] ; /* remove i from children of p */
1824  stack [++top] = i ; /* start dfs on child node i */
1825  }
1826  }
1827  return (k) ;
1828 }
1829 /* C = A' */
1830 cs *cs_transpose (const cs *A, csi values)
1831 {
1832  csi p, q, j, *Cp, *Ci, n, m, *Ap, *Ai, *w ;
1833  double *Cx, *Ax ;
1834  cs *C ;
1835  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
1836  m = A->m ; n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
1837  C = cs_spalloc (n, m, Ap [n], values && Ax, 0) ; /* allocate result */
1838  w = cs_calloc (m, sizeof (csi)) ; /* get workspace */
1839  if (!C || !w) return (cs_done (C, w, NULL, 0)) ; /* out of memory */
1840  Cp = C->p ; Ci = C->i ; Cx = C->x ;
1841  for (p = 0 ; p < Ap [n] ; p++) w [Ai [p]]++ ; /* row counts */
1842  cs_cumsum (Cp, w, m) ; /* row pointers */
1843  for (j = 0 ; j < n ; j++)
1844  {
1845  for (p = Ap [j] ; p < Ap [j+1] ; p++)
1846  {
1847  Ci [q = w [Ai [p]]++] = j ; /* place A(i,j) as entry C(j,i) */
1848  if (Cx) Cx [q] = Ax [p] ;
1849  }
1850  }
1851  return (cs_done (C, w, NULL, 1)) ; /* success; free w and return C */
1852 }
1853 /* sparse Cholesky update/downdate, L*L' + sigma*w*w' (sigma = +1 or -1) */
1854 csi cs_updown (cs *L, csi sigma, const cs *C, const csi *parent)
1855 {
1856  csi n, p, f, j, *Lp, *Li, *Cp, *Ci ;
1857  double *Lx, *Cx, alpha, beta = 1, delta, gamma, w1, w2, *w, beta2 = 1 ;
1858  if (!CS_CSC (L) || !CS_CSC (C) || !parent) return (0) ; /* check inputs */
1859  Lp = L->p ; Li = L->i ; Lx = L->x ; n = L->n ;
1860  Cp = C->p ; Ci = C->i ; Cx = C->x ;
1861  if ((p = Cp [0]) >= Cp [1]) return (1) ; /* return if C empty */
1862  w = cs_malloc (n, sizeof (double)) ; /* get workspace */
1863  if (!w) return (0) ; /* out of memory */
1864  f = Ci [p] ;
1865  for ( ; p < Cp [1] ; p++) f = CS_MIN (f, Ci [p]) ; /* f = min (find (C)) */
1866  for (j = f ; j != -1 ; j = parent [j]) w [j] = 0 ; /* clear workspace w */
1867  for (p = Cp [0] ; p < Cp [1] ; p++) w [Ci [p]] = Cx [p] ; /* w = C */
1868  for (j = f ; j != -1 ; j = parent [j]) /* walk path f up to root */
1869  {
1870  p = Lp [j] ;
1871  alpha = w [j] / Lx [p] ; /* alpha = w(j) / L(j,j) */
1872  beta2 = beta*beta + sigma*alpha*alpha ;
1873  if (beta2 <= 0) break ; /* not positive definite */
1874  beta2 = sqrt (beta2) ;
1875  delta = (sigma > 0) ? (beta / beta2) : (beta2 / beta) ;
1876  gamma = sigma * alpha / (beta2 * beta) ;
1877  Lx [p] = delta * Lx [p] + ((sigma > 0) ? (gamma * w [j]) : 0) ;
1878  beta = beta2 ;
1879  for (p++ ; p < Lp [j+1] ; p++)
1880  {
1881  w1 = w [Li [p]] ;
1882  w [Li [p]] = w2 = w1 - alpha * Lx [p] ;
1883  Lx [p] = delta * Lx [p] + gamma * ((sigma > 0) ? w1 : w2) ;
1884  }
1885  }
1886  cs_free (w) ;
1887  return (beta2 > 0) ;
1888 }
1889 /* solve Ux=b where x and b are dense. x=b on input, solution on output. */
1890 csi cs_usolve (const cs *U, double *x)
1891 {
1892  csi p, j, n, *Up, *Ui ;
1893  double *Ux ;
1894  if (!CS_CSC (U) || !x) return (0) ; /* check inputs */
1895  n = U->n ; Up = U->p ; Ui = U->i ; Ux = U->x ;
1896  for (j = n-1 ; j >= 0 ; j--)
1897  {
1898  x [j] /= Ux [Up [j+1]-1] ;
1899  for (p = Up [j] ; p < Up [j+1]-1 ; p++)
1900  {
1901  x [Ui [p]] -= Ux [p] * x [j] ;
1902  }
1903  }
1904  return (1) ;
1905 }
1906 /* allocate a sparse matrix (triplet form or compressed-column form) */
1907 cs *cs_spalloc (csi m, csi n, csi nzmax, csi values, csi triplet)
1908 {
1909  cs *A = cs_calloc (1, sizeof (cs)) ; /* allocate the cs struct */
1910  if (!A) return (NULL) ; /* out of memory */
1911  A->m = m ; /* define dimensions and nzmax */
1912  A->n = n ;
1913  A->nzmax = nzmax = CS_MAX (nzmax, 1) ;
1914  A->nz = triplet ? 0 : -1 ; /* allocate triplet or comp.col */
1915  A->p = cs_malloc (triplet ? nzmax : n+1, sizeof (csi)) ;
1916  A->i = cs_malloc (nzmax, sizeof (csi)) ;
1917  A->x = values ? cs_malloc (nzmax, sizeof (double)) : NULL ;
1918  return ((!A->p || !A->i || (values && !A->x)) ? cs_spfree (A) : A) ;
1919 }
1920
1921 /* change the max # of entries sparse matrix */
1922 csi cs_sprealloc (cs *A, csi nzmax)
1923 {
1924  csi ok, oki, okj = 1, okx = 1 ;
1925  if (!A) return (0) ;
1926  if (nzmax <= 0) nzmax = (CS_CSC (A)) ? (A->p [A->n]) : A->nz ;
1927  A->i = cs_realloc (A->i, nzmax, sizeof (csi), &oki) ;
1928  if (CS_TRIPLET (A)) A->p = cs_realloc (A->p, nzmax, sizeof (csi), &okj) ;
1929  if (A->x) A->x = cs_realloc (A->x, nzmax, sizeof (double), &okx) ;
1930  ok = (oki && okj && okx) ;
1931  if (ok) A->nzmax = nzmax ;
1932  return (ok) ;
1933 }
1934
1935 /* free a sparse matrix */
1937 {
1938  if (!A) return (NULL) ; /* do nothing if A already NULL */
1939  cs_free (A->p) ;
1940  cs_free (A->i) ;
1941  cs_free (A->x) ;
1942  return ((cs *) cs_free (A)) ; /* free the cs struct and return NULL */
1943 }
1944
1945 /* free a numeric factorization */
1947 {
1948  if (!N) return (NULL) ; /* do nothing if N already NULL */
1949  cs_spfree (N->L) ;
1950  cs_spfree (N->U) ;
1951  cs_free (N->pinv) ;
1952  cs_free (N->B) ;
1953  return ((csn *) cs_free (N)) ; /* free the csn struct and return NULL */
1954 }
1955
1956 /* free a symbolic factorization */
1958 {
1959  if (!S) return (NULL) ; /* do nothing if S already NULL */
1960  cs_free (S->pinv) ;
1961  cs_free (S->q) ;
1962  cs_free (S->parent) ;
1963  cs_free (S->cp) ;
1964  cs_free (S->leftmost) ;
1965  return ((css *) cs_free (S)) ; /* free the css struct and return NULL */
1966 }
1967
1968 /* allocate a cs_dmperm or cs_scc result */
1970 {
1971  csd *D ;
1972  D = cs_calloc (1, sizeof (csd)) ;
1973  if (!D) return (NULL) ;
1974  D->p = cs_malloc (m, sizeof (csi)) ;
1975  D->r = cs_malloc (m+6, sizeof (csi)) ;
1976  D->q = cs_malloc (n, sizeof (csi)) ;
1977  D->s = cs_malloc (n+6, sizeof (csi)) ;
1978  return ((!D->p || !D->r || !D->q || !D->s) ? cs_dfree (D) : D) ;
1979 }
1980
1981 /* free a cs_dmperm or cs_scc result */
1983 {
1984  if (!D) return (NULL) ; /* do nothing if D already NULL */
1985  cs_free (D->p) ;
1986  cs_free (D->q) ;
1987  cs_free (D->r) ;
1988  cs_free (D->s) ;
1989  return ((csd *) cs_free (D)) ; /* free the csd struct and return NULL */
1990 }
1991
1992 /* free workspace and return a sparse matrix result */
1993 cs *cs_done (cs *C, void *w, void *x, csi ok)
1994 {
1995  cs_free (w) ; /* free workspace */
1996  cs_free (x) ;
1997  return (ok ? C : cs_spfree (C)) ; /* return result if OK, else free it */
1998 }
1999
2000 /* free workspace and return csi array result */
2001 csi *cs_idone (csi *p, cs *C, void *w, csi ok)
2002 {
2003  cs_spfree (C) ; /* free temporary matrix */
2004  cs_free (w) ; /* free workspace */
2005  return (ok ? p : (csi *) cs_free (p)) ; /* return result, or free it */
2006 }
2007
2008 /* free workspace and return a numeric factorization (Cholesky, LU, or QR) */
2009 csn *cs_ndone (csn *N, cs *C, void *w, void *x, csi ok)
2010 {
2011  cs_spfree (C) ; /* free temporary matrix */
2012  cs_free (w) ; /* free workspace */
2013  cs_free (x) ;
2014  return (ok ? N : cs_nfree (N)) ; /* return result if OK, else free it */
2015 }
2016
2017 /* free workspace and return a csd result */
2018 csd *cs_ddone (csd *D, cs *C, void *w, csi ok)
2019 {
2020  cs_spfree (C) ; /* free temporary matrix */
2021  cs_free (w) ; /* free workspace */
2022  return (ok ? D : cs_dfree (D)) ; /* return result if OK, else free it */
2023 }
2024 /* solve U'x=b where x and b are dense. x=b on input, solution on output. */
2025 csi cs_utsolve (const cs *U, double *x)
2026 {
2027  csi p, j, n, *Up, *Ui ;
2028  double *Ux ;
2029  if (!CS_CSC (U) || !x) return (0) ; /* check inputs */
2030  n = U->n ; Up = U->p ; Ui = U->i ; Ux = U->x ;
2031  for (j = 0 ; j < n ; j++)
2032  {
2033  for (p = Up [j] ; p < Up [j+1]-1 ; p++)
2034  {
2035  x [j] -= Ux [p] * x [Ui [p]] ;
2036  }
2037  x [j] /= Ux [Up [j+1]-1] ;
2038  }
2039  return (1) ;
2040 }
csn * cs_lu(const cs *A, const css *S, double tol)
Definition: cs.c:1019
csi cs_qrsol(csi order, const cs *A, double *b)
Definition: cs.c:1477
csd * cs_ddone(csd *D, cs *C, void *w, csi ok)
Definition: cs.c:2018
static void init_ata(cs *AT, const csi *post, csi *w, csi **head, csi **next)
Definition: cs.c:498
static csi cs_diag(csi i, csi j, double aij, void *other)
Definition: cs.c:42
csi cs_pvec(const csi *p, const double *b, double *x, csi n)
Definition: cs.c:1397
#define CS_MARKED(w, j)
Definition: cs.h:155
csi * cs_idone(csi *p, cs *C, void *w, csi ok)
Definition: cs.c:2001
css * cs_schol(csi order, const cs *A)
Definition: cs.c:1631
static void cs_augment(csi k, const cs *A, csi *jmatch, csi *cheap, csi *w, csi *js, csi *is, csi *ps)
Definition: cs.c:1163
csi cs_reach(cs *G, const cs *B, csi k, csi *xi, const csi *pinv)
Definition: cs.c:1553
csi cs_print(const cs *A, csi brief)
Definition: cs.c:1359
csn * cs_qr(const cs *A, const css *S)
Definition: cs.c:1405
csi cs_utsolve(const cs *U, double *x)
Definition: cs.c:2025
csi cs_scatter(const cs *A, csi j, double beta, csi *w, double *x, csi mark, cs *C, csi nz)
Definition: cs.c:1570
csi cs_gaxpy(const cs *A, const double *x, double *y)
Definition: cs.c:886
double cs_cumsum(csi *p, csi *c, csi n)
Definition: cs.c:556
csi * cs_amd(csi order, const cs *A)
Definition: cs.c:45
static void cs_matched(csi n, const csi *wj, const csi *imatch, csi *p, csi *q, csi *cc, csi *rr, csi set, csi mark)
Definition: cs.c:641
csi * pinv
Definition: cs.h:83
csi n
Definition: cs.h:37
csi cs_entry(cs *T, csi i, csi j, double x)
Definition: cs.c:799
void * cs_realloc(void *p, csi n, size_t size, csi *ok)
Definition: cs.c:1155
#define CS_MARK(w, j)
Definition: cs.h:156
csi cs_cholsol(csi order, const cs *A, double *b)
Definition: cs.c:450
csi * cp
Definition: cs.h:72
double * B
Definition: cs.h:84
void * cs_calloc(csi n, size_t size)
Definition: cs.c:1142
csi m2
Definition: cs.h:74
csi m
Definition: cs.h:36
csi * q
Definition: cs.h:90
cs * cs_done(cs *C, void *w, void *x, csi ok)
Definition: cs.c:1993
csd * cs_dalloc(csi m, csi n)
Definition: cs.c:1969
csi * p
Definition: cs.h:38
#define CS_SUBSUB
Definition: cs.h:18
void * cs_free(void *p)
Definition: cs.c:1148
cs * U
Definition: cs.h:82
csi cs_droptol(cs *A, double tol)
Definition: cs.c:753
csi * cs_counts(const cs *A, const csi *parent, const csi *post, csi ata)
Definition: cs.c:510
csi cs_lusol(csi order, const cs *A, double *b, double tol)
Definition: cs.c:1104
csi cs_ltsolve(const cs *L, double *x)
Definition: cs.c:1002
cs * cs_symperm(const cs *A, const csi *pinv, csi values)
Definition: cs.c:1769
csi * i
Definition: cs.h:39
#define CS_VER
Definition: cs.h:16
#define CS_SUBVER
Definition: cs.h:17
double lnz
Definition: cs.h:75
Definition: cs.h:79
#define CS_FLIP(i)
Definition: cs.h:153
csi * parent
Definition: cs.h:71
csi * cs_pinv(csi const *p, csi n)
Definition: cs.c:1326
static csi cs_bfs(const cs *A, csi n, csi *wi, csi *wj, csi *queue, const csi *imatch, const csi *jmatch, csi mark)
Definition: cs.c:607
cs * cs_permute(const cs *A, const csi *pinv, const csi *q, csi values)
Definition: cs.c:1302
#define CS_MAX(a, b)
Definition: cs.h:151
css * cs_sfree(css *S)
Definition: cs.c:1957
csi * leftmost
Definition: cs.h:73
csi cs_leaf(csi i, csi j, const csi *first, csi *maxfirst, csi *prevleaf, csi *ancestor, csi *jleaf)
Definition: cs.c:950
double unz
Definition: cs.h:76
csi rr[5]
Definition: cs.h:94
Definition: cs.c:496
double cs_house(double *x, double *beta, csi n)
Definition: cs.c:921
csi cs_ereach(const cs *A, csi k, const csi *parent, csi *s, csi *w)
Definition: cs.c:811
static void cs_unmatched(csi m, const csi *wi, csi *p, csi *rr, csi set)
Definition: cs.c:657
Definition: cs.c:971
Definition: cs.h:20
csn * cs_chol(const cs *A, const css *S)
Definition: cs.c:393
double cs_norm(const cs *A)
Definition: cs.c:1288
#define CS_UNFLIP(i)
Definition: cs.h:154
#define CS_CSC(A)
Definition: cs.h:157
csd * cs_dmperm(const cs *A, csi seed)
Definition: cs.c:672
csi cc[5]
Definition: cs.h:95
csi cs_dropzeros(cs *A)
Definition: cs.c:761
csi * pinv
Definition: cs.h:69
csi cs_dfs(csi j, cs *G, csi top, csi *xi, csi *pstack, const csi *pinv)
Definition: cs.c:572
csi cs_usolve(const cs *U, double *x)
Definition: cs.c:1890
csi nzmax
Definition: cs.h:35
csd * cs_dfree(csd *D)
Definition: cs.c:1982
csi * p
Definition: cs.h:89
Definition: cs.h:33
csi * cs_post(const csi *parent, csi n)
Definition: cs.c:1336
csi * cs_maxtrans(const cs *A, csi seed)
Definition: cs.c:1204
csi * q
Definition: cs.h:70
csi * cs_etree(const cs *A, csi ata)
Definition: cs.c:833
cs * cs_add(const cs *A, const cs *B, double alpha, double beta)
Definition: cs.c:3
csi cs_spsolve(cs *G, const cs *B, csi k, csi *xi, double *x, const csi *pinv, csi lo)
Definition: cs.c:1656
#define csi
Definition: cs.h:27
cs * cs_multiply(const cs *A, const cs *B)
Definition: cs.c:1254
csi nz
Definition: cs.h:41
csi cs_tdfs(csi j, csi k, csi *head, const csi *next, csi *post, csi *stack)
Definition: cs.c:1807
#define NEXT(J)
Definition: cs.c:497
csi cs_ipvec(const csi *p, const double *b, double *x, csi n)
Definition: cs.c:942
static csi cs_rprune(csi i, csi j, double aij, void *other)
Definition: cs.c:665
csn * cs_nfree(csn *N)
Definition: cs.c:1946
cs * cs_spfree(cs *A)
Definition: cs.c:1936
static csi cs_tol(csi i, csi j, double aij, void *tol)
Definition: cs.c:749
static csi cs_wclear(csi mark, csi lemax, csi *w, csi n)
Definition: cs.c:30
#define CS_TRIPLET(A)
Definition: cs.h:158
cholmod_common c
Definition: chm_common.c:15
csd * cs_scc(cs *A)
Definition: cs.c:1591
#define CS_MIN(a, b)
Definition: cs.h:152
double * x
Definition: cs.h:40
csi cs_lsolve(const cs *L, double *x)
Definition: cs.c:985
csi * r
Definition: cs.h:91
csi cs_dupl(cs *A)
Definition: cs.c:766
static csi cs_nonzero(csi i, csi j, double aij, void *other)
Definition: cs.c:757
static csi cs_vcount(const cs *A, css *S)
Definition: cs.c:1683
csi * s
Definition: cs.h:92
csi cs_sprealloc(cs *A, csi nzmax)
Definition: cs.c:1922
void * cs_malloc(csi n, size_t size)
Definition: cs.c:1136
css * cs_sqr(csi order, const cs *A, csi qr)
Definition: cs.c:1740
cs * cs_compress(const cs *T)
Definition: cs.c:475
csi * cs_randperm(csi n, csi seed)
Definition: cs.c:1531
cs * L
Definition: cs.h:81
Definition: cs.h:67
cs * cs_transpose(const cs *A, csi values)
Definition: cs.c:1830
csi cs_fkeep(cs *A, csi(*fkeep)(csi, csi, double, void *), void *other)
Definition: cs.c:862
csi cs_updown(cs *L, csi sigma, const cs *C, const csi *parent)
Definition: cs.c:1854
csi cs_happly(const cs *V, csi i, double beta, double *x)
Definition: cs.c:902
#define CS_DATE
Definition: cs.h:19
cs * cs_spalloc(csi m, csi n, csi nzmax, csi values, csi triplet)
Definition: cs.c:1907
csn * cs_ndone(csn *N, cs *C, void *w, void *x, csi ok)
Definition: cs.c:2009