00001 #include "cs.h"
00002
00003 csn *cs_chol (const cs *A, const css *S)
00004 {
00005 double d, lki, *Lx, *x, *Cx ;
00006 int top, i, p, k, n, *Li, *Lp, *cp, *pinv, *s, *c, *parent, *Cp, *Ci ;
00007 cs *L, *C, *E ;
00008 csn *N ;
00009 if (!CS_CSC (A) || !S || !S->cp || !S->parent) return (NULL) ;
00010 n = A->n ;
00011 N = cs_calloc (1, sizeof (csn)) ;
00012 c = cs_malloc (2*n, sizeof (int)) ;
00013 x = cs_malloc (n, sizeof (double)) ;
00014 cp = S->cp ; pinv = S->pinv ; parent = S->parent ;
00015 C = pinv ? cs_symperm (A, pinv, 1) : ((cs *) A) ;
00016 E = pinv ? C : NULL ;
00017 if (!N || !c || !x || !C) return (cs_ndone (N, E, c, x, 0)) ;
00018 s = c + n ;
00019 Cp = C->p ; Ci = C->i ; Cx = C->x ;
00020 N->L = L = cs_spalloc (n, n, cp [n], 1, 0) ;
00021 if (!L) return (cs_ndone (N, E, c, x, 0)) ;
00022 Lp = L->p ; Li = L->i ; Lx = L->x ;
00023 for (k = 0 ; k < n ; k++) Lp [k] = c [k] = cp [k] ;
00024 for (k = 0 ; k < n ; k++)
00025 {
00026
00027 top = cs_ereach (C, k, parent, s, c) ;
00028 x [k] = 0 ;
00029 for (p = Cp [k] ; p < Cp [k+1] ; p++)
00030 {
00031 if (Ci [p] <= k) x [Ci [p]] = Cx [p] ;
00032 }
00033 d = x [k] ;
00034 x [k] = 0 ;
00035
00036 for ( ; top < n ; top++)
00037 {
00038 i = s [top] ;
00039 lki = x [i] / Lx [Lp [i]] ;
00040 x [i] = 0 ;
00041 for (p = Lp [i] + 1 ; p < c [i] ; p++)
00042 {
00043 x [Li [p]] -= Lx [p] * lki ;
00044 }
00045 d -= lki * lki ;
00046 p = c [i]++ ;
00047 Li [p] = k ;
00048 Lx [p] = lki ;
00049 }
00050
00051 if (d <= 0) return (cs_ndone (N, E, c, x, 0)) ;
00052 p = c [k]++ ;
00053 Li [p] = k ;
00054 Lx [p] = sqrt (d) ;
00055 }
00056 Lp [n] = cp [n] ;
00057 return (cs_ndone (N, E, c, x, 1)) ;
00058 }