cs_lu.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* [L,U,pinv]=lu(A, [q lnz unz]). lnz and unz can be guess */
3 csn *cs_lu (const cs *A, const css *S, double tol)
4 {
5  cs *L, *U ;
6  csn *N ;
7  double pivot, *Lx, *Ux, *x, a, t ;
8  int *Lp, *Li, *Up, *Ui, *pinv, *xi, *q, n, ipiv, k, top, p, i, col, lnz,unz;
9  if (!CS_CSC (A) || !S) return (NULL) ; /* check inputs */
10  n = A->n ;
11  q = S->q ; lnz = (int)S->lnz ; unz = (int)S->unz ;
12  x = (double*) cs_malloc (n, sizeof (double)) ; /* get double workspace */
13  xi = (int*) cs_malloc (2*n, sizeof (int)) ; /* get int workspace */
14  N = (csn*) cs_calloc (1, sizeof (csn)) ; /* allocate result */
15  if (!x || !xi || !N) return (cs_ndone (N, NULL, xi, x, 0)) ;
16  N->L = L = cs_spalloc (n, n, lnz, 1, 0) ; /* allocate result L */
17  N->U = U = cs_spalloc (n, n, unz, 1, 0) ; /* allocate result U */
18  N->pinv = pinv = (int*) cs_malloc (n, sizeof (int)) ; /* allocate result pinv */
19  if (!L || !U || !pinv) return (cs_ndone (N, NULL, xi, x, 0)) ;
20  Lp = L->p ; Up = U->p ;
21  for (i = 0 ; i < n ; i++) x [i] = 0 ; /* clear workspace */
22  for (i = 0 ; i < n ; i++) pinv [i] = -1 ; /* no rows pivotal yet */
23  for (k = 0 ; k <= n ; k++) Lp [k] = 0 ; /* no cols of L yet */
24  lnz = unz = 0 ;
25  for (k = 0 ; k < n ; k++) /* compute L(:,k) and U(:,k) */
26  {
27  /* --- Triangular solve --------------------------------------------- */
28  Lp [k] = lnz ; /* L(:,k) starts here */
29  Up [k] = unz ; /* U(:,k) starts here */
30  if ((lnz + n > L->nzmax && !cs_sprealloc (L, 2*L->nzmax + n)) ||
31  (unz + n > U->nzmax && !cs_sprealloc (U, 2*U->nzmax + n)))
32  {
33  return (cs_ndone (N, NULL, xi, x, 0)) ;
34  }
35  Li = L->i ; Lx = L->x ; Ui = U->i ; Ux = U->x ;
36  col = q ? (q [k]) : k ;
37  top = cs_spsolve (L, A, col, xi, x, pinv, 1) ; /* x = L\A(:,col) */
38  /* --- Find pivot --------------------------------------------------- */
39  ipiv = -1 ;
40  a = -1 ;
41  for (p = top ; p < n ; p++)
42  {
43  i = xi [p] ; /* x(i) is nonzero */
44  if (pinv [i] < 0) /* row i is not yet pivotal */
45  {
46  if ((t = fabs (x [i])) > a)
47  {
48  a = t ; /* largest pivot candidate so far */
49  ipiv = i ;
50  }
51  }
52  else /* x(i) is the entry U(pinv[i],k) */
53  {
54  Ui [unz] = pinv [i] ;
55  Ux [unz++] = x [i] ;
56  }
57  }
58  if (ipiv == -1 || a <= 0) return (cs_ndone (N, NULL, xi, x, 0)) ;
59  if (pinv [col] < 0 && fabs (x [col]) >= a*tol) ipiv = col ;
60  /* --- Divide by pivot ---------------------------------------------- */
61  pivot = x [ipiv] ; /* the chosen pivot */
62  Ui [unz] = k ; /* last entry in U(:,k) is U(k,k) */
63  Ux [unz++] = pivot ;
64  pinv [ipiv] = k ; /* ipiv is the kth pivot row */
65  Li [lnz] = ipiv ; /* first entry in L(:,k) is L(k,k) = 1 */
66  Lx [lnz++] = 1 ;
67  for (p = top ; p < n ; p++) /* L(k+1:n,k) = x / pivot */
68  {
69  i = xi [p] ;
70  if (pinv [i] < 0) /* x(i) is an entry in L(:,k) */
71  {
72  Li [lnz] = i ; /* save unpermuted row in L */
73  Lx [lnz++] = x [i] / pivot ; /* scale pivot column */
74  }
75  x [i] = 0 ; /* x [0..n-1] = 0 for next k */
76  }
77  }
78  /* --- Finalize L and U ------------------------------------------------- */
79  Lp [n] = lnz ;
80  Up [n] = unz ;
81  Li = L->i ; /* fix row indices of L for final pinv */
82  for (p = 0 ; p < lnz ; p++) Li [p] = pinv [Li [p]] ;
83  cs_sprealloc (L, 0) ; /* remove extra space from L and U */
84  cs_sprealloc (U, 0) ;
85  return (cs_ndone (N, NULL, xi, x, 1)) ; /* success */
86 }
#define N
int * pinv
Definition: cs.h:66
cs * cs_spalloc(int m, int n, int nzmax, int values, int triplet)
Definition: cs_util.c:3
int n
Definition: cs.h:20
cs * U
Definition: cs.h:65
double lnz
Definition: cs.h:58
Definition: cs.h:62
int nzmax
Definition: cs.h:18
int cs_sprealloc(cs *A, int nzmax)
Definition: cs_util.c:18
double unz
Definition: cs.h:59
#define CS_CSC(A)
Definition: cs.h:140
int * p
Definition: cs.h:21
Definition: cs.h:16
int * i
Definition: cs.h:22
void * cs_calloc(int n, size_t size)
Definition: cs_malloc.c:16
int * q
Definition: cs.h:53
void * cs_malloc(int n, size_t size)
Definition: cs_malloc.c:10
#define L
double * x
Definition: cs.h:23
csn * cs_ndone(csn *N, cs *C, void *w, void *x, int ok)
Definition: cs_util.c:105
ColXpr col(Index i)
Definition: BlockMethods.h:708
int cs_spsolve(cs *G, const cs *B, int k, int *xi, double *x, const int *pinv, int lo)
Definition: cs_spsolve.c:3
csn * cs_lu(const cs *A, const css *S, double tol)
Definition: cs_lu.c:3
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