EVSL  1.1.0
EigenValues Slicing Library
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Macros
cs_dmperm.c
Go to the documentation of this file.
1 #include "cs.h"
2 /* breadth-first search for coarse decomposition (C0,C1,R1 or R0,R3,C3) */
3 static CS_INT cs_bfs (const cs *A, CS_INT n, CS_INT *wi, CS_INT *wj, CS_INT *queue,
4  const CS_INT *imatch, const CS_INT *jmatch, CS_INT mark)
5 {
6  CS_INT *Ap, *Ai, head = 0, tail = 0, j, i, p, j2 ;
7  cs *C ;
8  for (j = 0 ; j < n ; j++) /* place all unmatched nodes in queue */
9  {
10  if (imatch [j] >= 0) continue ; /* skip j if matched */
11  wj [j] = 0 ; /* j in set C0 (R0 if transpose) */
12  queue [tail++] = j ; /* place unmatched col j in queue */
13  }
14  if (tail == 0) return (1) ; /* quick return if no unmatched nodes */
15  C = (mark == 1) ? ((cs *) A) : cs_transpose (A, 0) ;
16  if (!C) return (0) ; /* bfs of C=A' to find R3,C3 from R0 */
17  Ap = C->p ; Ai = C->i ;
18  while (head < tail) /* while queue is not empty */
19  {
20  j = queue [head++] ; /* get the head of the queue */
21  for (p = Ap [j] ; p < Ap [j+1] ; p++)
22  {
23  i = Ai [p] ;
24  if (wi [i] >= 0) continue ; /* skip if i is marked */
25  wi [i] = mark ; /* i in set R1 (C3 if transpose) */
26  j2 = jmatch [i] ; /* traverse alternating path to j2 */
27  if (wj [j2] >= 0) continue ;/* skip j2 if it is marked */
28  wj [j2] = mark ; /* j2 in set C1 (R3 if transpose) */
29  queue [tail++] = j2 ; /* add j2 to queue */
30  }
31  }
32  if (mark != 1) cs_spfree (C) ; /* free A' if it was created */
33  return (1) ;
34 }
35 
36 /* collect matched rows and columns into p and q */
37 static void cs_matched (CS_INT n, const CS_INT *wj, const CS_INT *imatch, CS_INT *p, CS_INT *q,
38  CS_INT *cc, CS_INT *rr, CS_INT set, CS_INT mark)
39 {
40  CS_INT kc = cc [set], j ;
41  CS_INT kr = rr [set-1] ;
42  for (j = 0 ; j < n ; j++)
43  {
44  if (wj [j] != mark) continue ; /* skip if j is not in C set */
45  p [kr++] = imatch [j] ;
46  q [kc++] = j ;
47  }
48  cc [set+1] = kc ;
49  rr [set] = kr ;
50 }
51 
52 /* collect unmatched rows into the permutation vector p */
53 static void cs_unmatched (CS_INT m, const CS_INT *wi, CS_INT *p, CS_INT *rr, CS_INT set)
54 {
55  CS_INT i, kr = rr [set] ;
56  for (i = 0 ; i < m ; i++) if (wi [i] == 0) p [kr++] = i ;
57  rr [set+1] = kr ;
58 }
59 
60 /* return 1 if row i is in R2 */
61 static CS_INT cs_rprune (CS_INT i, CS_INT j, CS_ENTRY aij, void *other)
62 {
63  CS_INT *rr = (CS_INT *) other ;
64  return (i >= rr [1] && i < rr [2]) ;
65 }
66 
67 /* Given A, compute coarse and then fine dmperm */
68 csd *cs_dmperm (const cs *A, CS_INT seed)
69 {
70  CS_INT m, n, i, j, k, cnz, nc, *jmatch, *imatch, *wi, *wj, *pinv, *Cp, *Ci,
71  *ps, *rs, nb1, nb2, *p, *q, *cc, *rr, *r, *s, ok ;
72  cs *C ;
73  csd *D, *scc ;
74  /* --- Maximum matching ------------------------------------------------- */
75  if (!CS_CSC (A)) return (NULL) ; /* check inputs */
76  m = A->m ; n = A->n ;
77  D = cs_dalloc (m, n) ; /* allocate result */
78  if (!D) return (NULL) ;
79  p = D->p ; q = D->q ; r = D->r ; s = D->s ; cc = D->cc ; rr = D->rr ;
80  jmatch = cs_maxtrans (A, seed) ; /* max transversal */
81  imatch = jmatch + m ; /* imatch = inverse of jmatch */
82  if (!jmatch) return (cs_ddone (D, NULL, jmatch, 0)) ;
83  /* --- Coarse decomposition --------------------------------------------- */
84  wi = r ; wj = s ; /* use r and s as workspace */
85  for (j = 0 ; j < n ; j++) wj [j] = -1 ; /* unmark all cols for bfs */
86  for (i = 0 ; i < m ; i++) wi [i] = -1 ; /* unmark all rows for bfs */
87  cs_bfs (A, n, wi, wj, q, imatch, jmatch, 1) ; /* find C1, R1 from C0*/
88  ok = cs_bfs (A, m, wj, wi, p, jmatch, imatch, 3) ; /* find R3, C3 from R0*/
89  if (!ok) return (cs_ddone (D, NULL, jmatch, 0)) ;
90  cs_unmatched (n, wj, q, cc, 0) ; /* unmatched set C0 */
91  cs_matched (n, wj, imatch, p, q, cc, rr, 1, 1) ; /* set R1 and C1 */
92  cs_matched (n, wj, imatch, p, q, cc, rr, 2, -1) ; /* set R2 and C2 */
93  cs_matched (n, wj, imatch, p, q, cc, rr, 3, 3) ; /* set R3 and C3 */
94  cs_unmatched (m, wi, p, rr, 3) ; /* unmatched set R0 */
95  cs_free (jmatch) ;
96  /* --- Fine decomposition ----------------------------------------------- */
97  pinv = cs_pinv (p, m) ; /* pinv=p' */
98  if (!pinv) return (cs_ddone (D, NULL, NULL, 0)) ;
99  C = cs_permute (A, pinv, q, 0) ;/* C=A(p,q) (it will hold A(R2,C2)) */
100  cs_free (pinv) ;
101  if (!C) return (cs_ddone (D, NULL, NULL, 0)) ;
102  Cp = C->p ;
103  nc = cc [3] - cc [2] ; /* delete cols C0, C1, and C3 from C */
104  if (cc [2] > 0) for (j = cc [2] ; j <= cc [3] ; j++) Cp [j-cc[2]] = Cp [j] ;
105  C->n = nc ;
106  if (rr [2] - rr [1] < m) /* delete rows R0, R1, and R3 from C */
107  {
108  cs_fkeep (C, cs_rprune, rr) ;
109  cnz = Cp [nc] ;
110  Ci = C->i ;
111  if (rr [1] > 0) for (k = 0 ; k < cnz ; k++) Ci [k] -= rr [1] ;
112  }
113  C->m = nc ;
114  scc = cs_scc (C) ; /* find strongly connected components of C*/
115  if (!scc) return (cs_ddone (D, C, NULL, 0)) ;
116  /* --- Combine coarse and fine decompositions --------------------------- */
117  ps = scc->p ; /* C(ps,ps) is the permuted matrix */
118  rs = scc->r ; /* kth block is rs[k]..rs[k+1]-1 */
119  nb1 = scc->nb ; /* # of blocks of A(R2,C2) */
120  for (k = 0 ; k < nc ; k++) wj [k] = q [ps [k] + cc [2]] ;
121  for (k = 0 ; k < nc ; k++) q [k + cc [2]] = wj [k] ;
122  for (k = 0 ; k < nc ; k++) wi [k] = p [ps [k] + rr [1]] ;
123  for (k = 0 ; k < nc ; k++) p [k + rr [1]] = wi [k] ;
124  nb2 = 0 ; /* create the fine block partitions */
125  r [0] = s [0] = 0 ;
126  if (cc [2] > 0) nb2++ ; /* leading coarse block A (R1, [C0 C1]) */
127  for (k = 0 ; k < nb1 ; k++) /* coarse block A (R2,C2) */
128  {
129  r [nb2] = rs [k] + rr [1] ; /* A (R2,C2) splits into nb1 fine blocks */
130  s [nb2] = rs [k] + cc [2] ;
131  nb2++ ;
132  }
133  if (rr [2] < m)
134  {
135  r [nb2] = rr [2] ; /* trailing coarse block A ([R3 R0], C3) */
136  s [nb2] = cc [3] ;
137  nb2++ ;
138  }
139  r [nb2] = m ;
140  s [nb2] = n ;
141  D->nb = nb2 ;
142  cs_dfree (scc) ;
143  return (cs_ddone (D, C, NULL, 1)) ;
144 }
csd * cs_scc(cs *A)
Definition: cs_scc.c:3
void * cs_free(void *p)
Definition: cs_malloc.c:22
CS_INT * cs_pinv(CS_INT const *p, CS_INT n)
Definition: cs_pinv.c:3
#define cs
Definition: cs.h:637
csd * cs_dmperm(const cs *A, CS_INT seed)
Definition: cs_dmperm.c:68
#define CS_ENTRY
Definition: cs.h:635
cs * cs_spfree(cs *A)
Definition: cs_util.c:33
#define csd
Definition: cs.h:690
#define CS_CSC(A)
Definition: cs.h:659
CS_INT cs_fkeep(cs *A, CS_INT(*fkeep)(CS_INT, CS_INT, CS_ENTRY, void *), void *other)
Definition: cs_fkeep.c:3
cs * cs_transpose(const cs *A, CS_INT values)
Definition: cs_transpose.c:3
#define CS_INT
Definition: cs.h:627
cs * cs_permute(const cs *A, const CS_INT *pinv, const CS_INT *q, CS_INT values)
Definition: cs_permute.c:3
csd * cs_dalloc(CS_INT m, CS_INT n)
Definition: cs_util.c:66
csd * cs_ddone(csd *D, cs *C, void *w, CS_INT ok)
Definition: cs_util.c:115
CS_INT * cs_maxtrans(const cs *A, CS_INT seed)
Definition: cs_maxtrans.c:44
csd * cs_dfree(csd *D)
Definition: cs_util.c:79