cs_chol.c
Go to the documentation of this file.
00001 #include "cs.h"
00002 /* L = chol (A, [pinv parent cp]), pinv is optional */
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)) ;       /* allocate result */
00012     c = cs_malloc (2*n, sizeof (int)) ;     /* get int workspace */
00013     x = cs_malloc (n, sizeof (double)) ;    /* get double workspace */
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 ;           /* E is alias for A, or a copy E=A(p,p) */
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) ;    /* allocate result */
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++)       /* compute L(k,:) for L*L' = C */
00025     {
00026         /* --- Nonzero pattern of L(k,:) ------------------------------------ */
00027         top = cs_ereach (C, k, parent, s, c) ;      /* find pattern of L(k,:) */
00028         x [k] = 0 ;                                 /* x (0:k) is now zero */
00029         for (p = Cp [k] ; p < Cp [k+1] ; p++)       /* x = full(triu(C(:,k))) */
00030         {
00031             if (Ci [p] <= k) x [Ci [p]] = Cx [p] ;
00032         }
00033         d = x [k] ;                     /* d = C(k,k) */
00034         x [k] = 0 ;                     /* clear x for k+1st iteration */
00035         /* --- Triangular solve --------------------------------------------- */
00036         for ( ; top < n ; top++)    /* solve L(0:k-1,0:k-1) * x = C(:,k) */
00037         {
00038             i = s [top] ;               /* s [top..n-1] is pattern of L(k,:) */
00039             lki = x [i] / Lx [Lp [i]] ; /* L(k,i) = x (i) / L(i,i) */
00040             x [i] = 0 ;                 /* clear x for k+1st iteration */
00041             for (p = Lp [i] + 1 ; p < c [i] ; p++)
00042             {
00043                 x [Li [p]] -= Lx [p] * lki ;
00044             }
00045             d -= lki * lki ;            /* d = d - L(k,i)*L(k,i) */
00046             p = c [i]++ ;
00047             Li [p] = k ;                /* store L(k,i) in column i */
00048             Lx [p] = lki ;
00049         }
00050         /* --- Compute L(k,k) ----------------------------------------------- */
00051         if (d <= 0) return (cs_ndone (N, E, c, x, 0)) ; /* not pos def */
00052         p = c [k]++ ;
00053         Li [p] = k ;                /* store L(k,k) = sqrt (d) in column k */
00054         Lx [p] = sqrt (d) ;
00055     }
00056     Lp [n] = cp [n] ;               /* finalize L */
00057     return (cs_ndone (N, E, c, x, 1)) ; /* success: free E,s,x; return N */
00058 }


re_vision
Author(s): Dorian Galvez-Lopez
autogenerated on Sun Jan 5 2014 11:30:59