cs_chol.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* L = chol (A, [pinv parent cp]), pinv is optional */
3 csn *cs_chol (const cs *A, const css *S)
4 {
5  double d, lki, *Lx, *x, *Cx ;
6  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 = (csn*) cs_calloc (1, sizeof (csn)) ; /* allocate result */
12  c = (int*) cs_malloc (2*n, sizeof (int)) ; /* get int workspace */
13  x = (double*) cs_malloc (n, sizeof (double)) ; /* get double 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 * 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] = lki ;
49  }
50  /* --- Compute L(k,k) ----------------------------------------------- */
51  if (d <= 0) return (cs_ndone (N, E, c, x, 0)) ; /* not pos def */
52  p = c [k]++ ;
53  Li [p] = k ; /* store L(k,k) = sqrt (d) in column k */
54  Lx [p] = sqrt (d) ;
55  }
56  Lp [n] = cp [n] ; /* finalize L */
57  return (cs_ndone (N, E, c, x, 1)) ; /* success: free E,s,x; return N */
58 }
#define N
IntermediateState sqrt(const Expression &arg)
cs * cs_spalloc(int m, int n, int nzmax, int values, int triplet)
Definition: cs_util.c:3
cs * cs_symperm(const cs *A, const int *pinv, int values)
Definition: cs_symperm.c:3
csn * cs_chol(const cs *A, const css *S)
Definition: cs_chol.c:3
int n
Definition: cs.h:20
Definition: cs.h:62
int cs_ereach(const cs *A, int k, const int *parent, int *s, int *w)
Definition: cs_ereach.c:3
#define CS_CSC(A)
Definition: cs.h:140
int * p
Definition: cs.h:21
#define E
Definition: cs.h:16
int * i
Definition: cs.h:22
int * cp
Definition: cs.h:55
void * cs_calloc(int n, size_t size)
Definition: cs_malloc.c:16
int * pinv
Definition: cs.h:52
void * cs_malloc(int n, size_t size)
Definition: cs_malloc.c:10
#define L
double * x
Definition: cs.h:23
int * parent
Definition: cs.h:54
csn * cs_ndone(csn *N, cs *C, void *w, void *x, int ok)
Definition: cs_util.c:105
cs * L
Definition: cs.h:64
Definition: cs.h:50


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:34:31