1 #include "cs.h"
2 /* L = chol (A, [pinv parent cp]), pinv is optional */
cs_chol(const cs * A,const css * S)3 csn *cs_chol (const cs *A, const css *S)
4 {
5 CS_ENTRY d, lki, *Lx, *x, *Cx ;
6 CS_INT top, i, p, k, n, *Li, *Lp, *cp, *pinv, *s, *c, *parent, *Cp, *Ci ;
7 cs *L, *C, *E ;
8 csn *N ;
9 if (!CS_CSC (A) || !S || !S->cp || !S->parent) return (NULL) ;
10 n = A->n ;
11 N = cs_calloc (1, sizeof (csn)) ; /* allocate result */
12 c = cs_malloc (2*n, sizeof (CS_INT)) ; /* get CS_INT workspace */
13 x = cs_malloc (n, sizeof (CS_ENTRY)) ; /* get CS_ENTRY workspace */
14 cp = S->cp ; pinv = S->pinv ; parent = S->parent ;
15 C = pinv ? cs_symperm (A, pinv, 1) : ((cs *) A) ;
16 E = pinv ? C : NULL ; /* E is alias for A, or a copy E=A(p,p) */
17 if (!N || !c || !x || !C) return (cs_ndone (N, E, c, x, 0)) ;
18 s = c + n ;
19 Cp = C->p ; Ci = C->i ; Cx = C->x ;
20 N->L = L = cs_spalloc (n, n, cp [n], 1, 0) ; /* allocate result */
21 if (!L) return (cs_ndone (N, E, c, x, 0)) ;
22 Lp = L->p ; Li = L->i ; Lx = L->x ;
23 for (k = 0 ; k < n ; k++) Lp [k] = c [k] = cp [k] ;
24 for (k = 0 ; k < n ; k++) /* compute L(k,:) for L*L' = C */
25 {
26 /* --- Nonzero pattern of L(k,:) ------------------------------------ */
27 top = cs_ereach (C, k, parent, s, c) ; /* find pattern of L(k,:) */
28 x [k] = 0 ; /* x (0:k) is now zero */
29 for (p = Cp [k] ; p < Cp [k+1] ; p++) /* x = full(triu(C(:,k))) */
30 {
31 if (Ci [p] <= k) x [Ci [p]] = Cx [p] ;
32 }
33 d = x [k] ; /* d = C(k,k) */
34 x [k] = 0 ; /* clear x for k+1st iteration */
35 /* --- Triangular solve --------------------------------------------- */
36 for ( ; top < n ; top++) /* solve L(0:k-1,0:k-1) * x = C(:,k) */
37 {
38 i = s [top] ; /* s [top..n-1] is pattern of L(k,:) */
39 lki = x [i] / Lx [Lp [i]] ; /* L(k,i) = x (i) / L(i,i) */
40 x [i] = 0 ; /* clear x for k+1st iteration */
41 for (p = Lp [i] + 1 ; p < c [i] ; p++)
42 {
43 x [Li [p]] -= Lx [p] * lki ;
44 }
45 d -= lki * CS_CONJ (lki) ; /* d = d - L(k,i)*L(k,i) */
46 p = c [i]++ ;
47 Li [p] = k ; /* store L(k,i) in column i */
48 Lx [p] = CS_CONJ (lki) ;
49 }
50 /* --- Compute L(k,k) ----------------------------------------------- */
51 if (CS_REAL (d) <= 0 || CS_IMAG (d) != 0)
52 return (cs_ndone (N, E, c, x, 0)) ; /* not pos def */
53 p = c [k]++ ;
54 Li [p] = k ; /* store L(k,k) = sqrt (d) in column k */
55 Lx [p] = sqrt (d) ;
56 }
57 Lp [n] = cp [n] ; /* finalize L */
58 return (cs_ndone (N, E, c, x, 1)) ; /* success: free E,s,x; return N */
59 }
60