cs_sqr.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* compute nnz(V) = S->lnz, S->pinv, S->leftmost, S->m2 from A and S->parent */
3 static int cs_vcount (const cs *A, css *S)
4 {
5  int i, k, p, pa, n = A->n, m = A->m, *Ap = A->p, *Ai = A->i, *next, *head,
6  *tail, *nque, *pinv, *leftmost, *w, *parent = S->parent ;
7  S->pinv = pinv = (int*)cs_malloc (m+n, sizeof (int)) ; /* allocate pinv, */
8  S->leftmost = leftmost = (int*)cs_malloc (m, sizeof (int)) ; /* and leftmost */
9  w = (int*)cs_malloc (m+3*n, sizeof (int)) ; /* get workspace */
10  if (!pinv || !w || !leftmost)
11  {
12  cs_free (w) ; /* pinv and leftmost freed later */
13  return (0) ; /* out of memory */
14  }
15  next = w ; head = w + m ; tail = w + m + n ; nque = w + m + 2*n ;
16  for (k = 0 ; k < n ; k++) head [k] = -1 ; /* queue k is empty */
17  for (k = 0 ; k < n ; k++) tail [k] = -1 ;
18  for (k = 0 ; k < n ; k++) nque [k] = 0 ;
19  for (i = 0 ; i < m ; i++) leftmost [i] = -1 ;
20  for (k = n-1 ; k >= 0 ; k--)
21  {
22  for (p = Ap [k] ; p < Ap [k+1] ; p++)
23  {
24  leftmost [Ai [p]] = k ; /* leftmost[i] = min(find(A(i,:)))*/
25  }
26  }
27  for (i = m-1 ; i >= 0 ; i--) /* scan rows in reverse order */
28  {
29  pinv [i] = -1 ; /* row i is not yet ordered */
30  k = leftmost [i] ;
31  if (k == -1) continue ; /* row i is empty */
32  if (nque [k]++ == 0) tail [k] = i ; /* first row in queue k */
33  next [i] = head [k] ; /* put i at head of queue k */
34  head [k] = i ;
35  }
36  S->lnz = 0 ;
37  S->m2 = m ;
38  for (k = 0 ; k < n ; k++) /* find row permutation and nnz(V)*/
39  {
40  i = head [k] ; /* remove row i from queue k */
41  S->lnz++ ; /* count V(k,k) as nonzero */
42  if (i < 0) i = S->m2++ ; /* add a fictitious row */
43  pinv [i] = k ; /* associate row i with V(:,k) */
44  if (--nque [k] <= 0) continue ; /* skip if V(k+1:m,k) is empty */
45  S->lnz += nque [k] ; /* nque [k] is nnz (V(k+1:m,k)) */
46  if ((pa = parent [k]) != -1) /* move all rows to parent of k */
47  {
48  if (nque [pa] == 0) tail [pa] = tail [k] ;
49  next [tail [k]] = head [pa] ;
50  head [pa] = next [i] ;
51  nque [pa] += nque [k] ;
52  }
53  }
54  for (i = 0 ; i < m ; i++) if (pinv [i] < 0) pinv [i] = k++ ;
55  cs_free (w) ;
56  return (1) ;
57 }
58 
59 /* symbolic ordering and analysis for QR or LU */
60 css *cs_sqr (int order, const cs *A, int qr)
61 {
62  int n, k, ok = 1, *post ;
63  css *S ;
64  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
65  n = A->n ;
66  S = (css*)cs_calloc (1, sizeof (css)) ; /* allocate result S */
67  if (!S) return (NULL) ; /* out of memory */
68  S->q = cs_amd (order, A) ; /* fill-reducing ordering */
69  if (order && !S->q) return (cs_sfree (S)) ;
70  if (qr) /* QR symbolic analysis */
71  {
72  cs *C = order ? cs_permute (A, NULL, S->q, 0) : ((cs *) A) ;
73  S->parent = cs_etree (C, 1) ; /* etree of C'*C, where C=A(:,q) */
74  post = cs_post (S->parent, n) ;
75  S->cp = cs_counts (C, S->parent, post, 1) ; /* col counts chol(C'*C) */
76  cs_free (post) ;
77  ok = C && S->parent && S->cp && cs_vcount (C, S) ;
78  if (ok) for (S->unz = 0, k = 0 ; k < n ; k++) S->unz += S->cp [k] ;
79  ok = ok && S->lnz >= 0 && S->unz >= 0 ; /* int overflow guard */
80  if (order) cs_spfree (C) ;
81  }
82  else
83  {
84  S->unz = 4*(A->p [n]) + n ; /* for LU factorization only, */
85  S->lnz = S->unz ; /* guess nnz(L) and nnz(U) */
86  }
87  return (ok ? S : cs_sfree (S)) ; /* return result S */
88 }
css * cs_sqr(int order, const cs *A, int qr)
Definition: cs_sqr.c:60
cs * cs_spfree(cs *A)
Definition: cs_util.c:32
cs * cs_permute(const cs *A, const int *pinv, const int *q, int values)
Definition: cs_permute.c:3
int n
Definition: cs.h:20
void * cs_free(void *p)
Definition: cs_malloc.c:22
int * cs_post(const int *parent, int n)
Definition: cs_post.c:3
double lnz
Definition: cs.h:58
double unz
Definition: cs.h:59
static int cs_vcount(const cs *A, css *S)
Definition: cs_sqr.c:3
int * cs_etree(const cs *A, int ata)
Definition: cs_etree.c:3
int * leftmost
Definition: cs.h:56
#define CS_CSC(A)
Definition: cs.h:140
SegmentReturnType tail(Index vecSize)
Definition: BlockMethods.h:810
int * p
Definition: cs.h:21
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
int * q
Definition: cs.h:53
void * cs_malloc(int n, size_t size)
Definition: cs_malloc.c:10
SegmentReturnType head(Index vecSize)
Definition: BlockMethods.h:781
int * cs_amd(int order, const cs *A)
Definition: cs_amd.c:18
int m2
Definition: cs.h:57
Expression next(const Expression &arg)
css * cs_sfree(css *S)
Definition: cs_util.c:53
int * parent
Definition: cs.h:54
int m
Definition: cs.h:19
int * cs_counts(const cs *A, const int *parent, const int *post, int ata)
Definition: cs_counts.c:17
Definition: cs.h:50


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