EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
cs_maxtrans.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* find an augmenting path starting at column k and extend the match if found */
3 static void cs_augment (CS_INT k, const cs *A, CS_INT *jmatch, CS_INT *cheap, CS_INT *w,
4  CS_INT *js, CS_INT *is, CS_INT *ps)
5 {
6  CS_INT found = 0, p, i = -1, *Ap = A->p, *Ai = A->i, head = 0, j ;
7  js [0] = k ; /* start with just node k in jstack */
8  while (head >= 0)
9  {
10  /* --- Start (or continue) depth-first-search at node j ------------- */
11  j = js [head] ; /* get j from top of jstack */
12  if (w [j] != k) /* 1st time j visited for kth path */
13  {
14  w [j] = k ; /* mark j as visited for kth path */
15  for (p = cheap [j] ; p < Ap [j+1] && !found ; p++)
16  {
17  i = Ai [p] ; /* try a cheap assignment (i,j) */
18  found = (jmatch [i] == -1) ;
19  }
20  cheap [j] = p ; /* start here next time j is traversed*/
21  if (found)
22  {
23  is [head] = i ; /* column j matched with row i */
24  break ; /* end of augmenting path */
25  }
26  ps [head] = Ap [j] ; /* no cheap match: start dfs for j */
27  }
28  /* --- Depth-first-search of neighbors of j ------------------------- */
29  for (p = ps [head] ; p < Ap [j+1] ; p++)
30  {
31  i = Ai [p] ; /* consider row i */
32  if (w [jmatch [i]] == k) continue ; /* skip jmatch [i] if marked */
33  ps [head] = p + 1 ; /* pause dfs of node j */
34  is [head] = i ; /* i will be matched with j if found */
35  js [++head] = jmatch [i] ; /* start dfs at column jmatch [i] */
36  break ;
37  }
38  if (p == Ap [j+1]) head-- ; /* node j is done; pop from stack */
39  } /* augment the match if path found: */
40  if (found) for (p = head ; p >= 0 ; p--) jmatch [is [p]] = js [p] ;
41 }
42 
43 /* find a maximum transveral */
44 CS_INT *cs_maxtrans (const cs *A, CS_INT seed) /*[jmatch [0..m-1]; imatch [0..n-1]]*/
45 {
46  CS_INT i, j, k, n, m, p, n2 = 0, m2 = 0, *Ap, *jimatch, *w, *cheap, *js, *is,
47  *ps, *Ai, *Cp, *jmatch, *imatch, *q ;
48  cs *C ;
49  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
50  n = A->n ; m = A->m ; Ap = A->p ; Ai = A->i ;
51  w = jimatch = cs_calloc (m+n, sizeof (CS_INT)) ; /* allocate result */
52  if (!jimatch) return (NULL) ;
53  for (k = 0, j = 0 ; j < n ; j++) /* count nonempty rows and columns */
54  {
55  n2 += (Ap [j] < Ap [j+1]) ;
56  for (p = Ap [j] ; p < Ap [j+1] ; p++)
57  {
58  w [Ai [p]] = 1 ;
59  k += (j == Ai [p]) ; /* count entries already on diagonal */
60  }
61  }
62  if (k == CS_MIN (m,n)) /* quick return if diagonal zero-free */
63  {
64  jmatch = jimatch ; imatch = jimatch + m ;
65  for (i = 0 ; i < k ; i++) jmatch [i] = i ;
66  for ( ; i < m ; i++) jmatch [i] = -1 ;
67  for (j = 0 ; j < k ; j++) imatch [j] = j ;
68  for ( ; j < n ; j++) imatch [j] = -1 ;
69  return (cs_idone (jimatch, NULL, NULL, 1)) ;
70  }
71  for (i = 0 ; i < m ; i++) m2 += w [i] ;
72  C = (m2 < n2) ? cs_transpose (A,0) : ((cs *) A) ; /* transpose if needed */
73  if (!C) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, NULL, 0)) ;
74  n = C->n ; m = C->m ; Cp = C->p ;
75  jmatch = (m2 < n2) ? jimatch + n : jimatch ;
76  imatch = (m2 < n2) ? jimatch : jimatch + m ;
77  w = cs_malloc (5*n, sizeof (CS_INT)) ; /* get workspace */
78  if (!w) return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 0)) ;
79  cheap = w + n ; js = w + 2*n ; is = w + 3*n ; ps = w + 4*n ;
80  for (j = 0 ; j < n ; j++) cheap [j] = Cp [j] ; /* for cheap assignment */
81  for (j = 0 ; j < n ; j++) w [j] = -1 ; /* all columns unflagged */
82  for (i = 0 ; i < m ; i++) jmatch [i] = -1 ; /* nothing matched yet */
83  q = cs_randperm (n, seed) ; /* q = random permutation */
84  for (k = 0 ; k < n ; k++) /* augment, starting at column q[k] */
85  {
86  cs_augment (q ? q [k]: k, C, jmatch, cheap, w, js, is, ps) ;
87  }
88  cs_free (q) ;
89  for (j = 0 ; j < n ; j++) imatch [j] = -1 ; /* find row match */
90  for (i = 0 ; i < m ; i++) if (jmatch [i] >= 0) imatch [jmatch [i]] = i ;
91  return (cs_idone (jimatch, (m2 < n2) ? C : NULL, w, 1)) ;
92 }
void * cs_free(void *p)
Definition: cs_malloc.c:22
#define cs
Definition: cs.h:637
#define CS_CSC(A)
Definition: cs.h:659
void * cs_calloc(CS_INT n, size_t size)
Definition: cs_malloc.c:16
CS_INT * cs_randperm(CS_INT n, CS_INT seed)
Definition: cs_randperm.c:5
cs * cs_transpose(const cs *A, CS_INT values)
Definition: cs_transpose.c:3
#define CS_INT
Definition: cs.h:627
#define CS_MIN(a, b)
Definition: cs.h:654
void * cs_malloc(CS_INT n, size_t size)
Definition: cs_malloc.c:10
CS_INT * cs_maxtrans(const cs *A, CS_INT seed)
Definition: cs_maxtrans.c:44
CS_INT * cs_idone(CS_INT *p, cs *C, void *w, CS_INT ok)
Definition: cs_util.c:98