cs_qr.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* sparse QR factorization [V,beta,pinv,R] = qr (A) */
3 csn *cs_qr (const cs *A, const css *S)
4 {
5  double *Rx, *Vx, *Ax, *x, *Beta ;
6  int i, k, p, n, vnz, p1, top, m2, len, col, rnz, *s, *leftmost, *Ap, *Ai,
7  *parent, *Rp, *Ri, *Vp, *Vi, *w, *pinv, *q ;
8  cs *R, *V ;
9  csn *N ;
10  if (!CS_CSC (A) || !S) return (NULL) ;
11  n = A->n ; Ap = A->p ; Ai = A->i ; Ax = A->x ;
12  q = S->q ; parent = S->parent ; pinv = S->pinv ; m2 = S->m2 ;
13  vnz = (int)S->lnz ; rnz = (int)S->unz ; leftmost = S->leftmost ;
14  w = (int*) cs_malloc (m2+n, sizeof (int)) ; /* get int workspace */
15  x = (double*) cs_malloc (m2, sizeof (double)) ; /* get double workspace */
16  N = (csn*) cs_calloc (1, sizeof (csn)) ; /* allocate result */
17  if (!w || !x || !N) return (cs_ndone (N, NULL, w, x, 0)) ;
18  s = w + m2 ; /* s is size n */
19  for (k = 0 ; k < m2 ; k++) x [k] = 0 ; /* clear workspace x */
20  N->L = V = cs_spalloc (m2, n, vnz, 1, 0) ; /* allocate result V */
21  N->U = R = cs_spalloc (m2, n, rnz, 1, 0) ; /* allocate result R */
22  N->B = Beta = (double*) cs_malloc (n, sizeof (double)) ; /* allocate result Beta */
23  if (!R || !V || !Beta) return (cs_ndone (N, NULL, w, x, 0)) ;
24  Rp = R->p ; Ri = R->i ; Rx = R->x ;
25  Vp = V->p ; Vi = V->i ; Vx = V->x ;
26  for (i = 0 ; i < m2 ; i++) w [i] = -1 ; /* clear w, to mark nodes */
27  rnz = 0 ; vnz = 0 ;
28  for (k = 0 ; k < n ; k++) /* compute V and R */
29  {
30  Rp [k] = rnz ; /* R(:,k) starts here */
31  Vp [k] = p1 = vnz ; /* V(:,k) starts here */
32  w [k] = k ; /* add V(k,k) to pattern of V */
33  Vi [vnz++] = k ;
34  top = n ;
35  col = q ? q [k] : k ;
36  for (p = Ap [col] ; p < Ap [col+1] ; p++) /* find R(:,k) pattern */
37  {
38  i = leftmost [Ai [p]] ; /* i = min(find(A(i,q))) */
39  for (len = 0 ; w [i] != k ; i = parent [i]) /* traverse up to k */
40  {
41  s [len++] = i ;
42  w [i] = k ;
43  }
44  while (len > 0) s [--top] = s [--len] ; /* push path on stack */
45  i = pinv [Ai [p]] ; /* i = permuted row of A(:,col) */
46  x [i] = Ax [p] ; /* x (i) = A(:,col) */
47  if (i > k && w [i] < k) /* pattern of V(:,k) = x (k+1:m) */
48  {
49  Vi [vnz++] = i ; /* add i to pattern of V(:,k) */
50  w [i] = k ;
51  }
52  }
53  for (p = top ; p < n ; p++) /* for each i in pattern of R(:,k) */
54  {
55  i = s [p] ; /* R(i,k) is nonzero */
56  cs_happly (V, i, Beta [i], x) ; /* apply (V(i),Beta(i)) to x */
57  Ri [rnz] = i ; /* R(i,k) = x(i) */
58  Rx [rnz++] = x [i] ;
59  x [i] = 0 ;
60  if (parent [i] == k) vnz = cs_scatter (V, i, 0, w, NULL, k, V, vnz);
61  }
62  for (p = p1 ; p < vnz ; p++) /* gather V(:,k) = x */
63  {
64  Vx [p] = x [Vi [p]] ;
65  x [Vi [p]] = 0 ;
66  }
67  Ri [rnz] = k ; /* R(k,k) = norm (x) */
68  Rx [rnz++] = cs_house (Vx+p1, Beta+k, vnz-p1) ; /* [v,beta]=house(x) */
69  }
70  Rp [n] = rnz ; /* finalize R */
71  Vp [n] = vnz ; /* finalize V */
72  return (cs_ndone (N, NULL, w, x, 1)) ; /* success */
73 }
#define N
double * B
Definition: cs.h:67
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 cs_house(double *x, double *beta, int n)
Definition: cs_house.c:4
double lnz
Definition: cs.h:58
Definition: cs.h:62
#define Beta
int cs_happly(const cs *V, int i, double beta, double *x)
Definition: cs_happly.c:3
double unz
Definition: cs.h:59
int * leftmost
Definition: cs.h:56
#define CS_CSC(A)
Definition: cs.h:140
int * p
Definition: cs.h:21
int cs_scatter(const cs *A, int j, double beta, int *w, double *x, int mark, cs *C, int nz)
Definition: cs_scatter.c:3
Definition: cs.h:16
int * i
Definition: cs.h:22
void * cs_calloc(int n, size_t size)
Definition: cs_malloc.c:16
int * pinv
Definition: cs.h:52
int * q
Definition: cs.h:53
void * cs_malloc(int n, size_t size)
Definition: cs_malloc.c:10
int m2
Definition: cs.h:57
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
ColXpr col(Index i)
Definition: BlockMethods.h:708
cs * L
Definition: cs.h:64
Definition: cs.h:50
#define R
csn * cs_qr(const cs *A, const css *S)
Definition: cs_qr.c:3


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