← Back to branch summary

~ubuntu-branches/ubuntu/oneiric/suitesparse/oneiric

« back to all changes in this revision

Viewing changes to CXSparse/Source/cs_dmperm.c

  • Committer: Bazaar Package Importer
  • Author(s): Nick Ellery
  • Date: 2009-06-14 19:15:52 UTC
  • mfrom: (7.2.2 sid)
  • Revision ID: james.westby@ubuntu.com-20090614191552-2hliya5q8n1quseu
Tags: 1:3.4.0-1ubuntu1
https://launchpad.net/bugs/387137
* Merge from debian unstable, remaining changes (LP: #387137):
  - debian/control:
    - demote libatlas-doc from recommends to suggests as it is not in main
    - drop recommends on doc-central as it is not in main

Show diffs side-by-side

added added

removed removed

Lines of Context:
CXSparse/Source/cs_dmperm.c
5
5
{
6
6
    CS_INT *Ap, *Ai, head = 0, tail = 0, j, i, p, j2 ;
7
7
    cs *C ;
8
 
    for (j = 0 ; j < n ; j++)           /* place all unmatched nodes in queue */
 
8
    for (j = 0 ; j < n ; j++)           /* place all unmatched nodes in queue */
9
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 */
 
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
13
    }
14
 
    if (tail == 0) return (1) ;         /* quick return if no unmatched nodes */
 
14
    if (tail == 0) return (1) ;         /* quick return if no unmatched nodes */
15
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 */
 
16
    if (!C) return (0) ;                /* bfs of C=A' to find R3,C3 from R0 */
17
17
    Ap = C->p ; Ai = C->i ;
18
 
    while (head < tail)                 /* while queue is not empty */
 
18
    while (head < tail)                 /* while queue is not empty */
19
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
 
        }
 
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
31
    }
32
 
    if (mark != 1) cs_spfree (C) ;      /* free A' if it was created */
 
32
    if (mark != 1) cs_spfree (C) ;      /* free A' if it was created */
33
33
    return (1) ;
34
34
}
35
35
 
41
41
    CS_INT kr = rr [set-1] ;
42
42
    for (j = 0 ; j < n ; j++)
43
43
    {
44
 
        if (wj [j] != mark) continue ;      /* skip if j is not in C set */
45
 
        p [kr++] = imatch [j] ;
46
 
        q [kc++] = j ;
 
44
        if (wj [j] != mark) continue ;      /* skip if j is not in C set */
 
45
        p [kr++] = imatch [j] ;
 
46
        q [kc++] = j ;
47
47
    }
48
48
    cc [set+1] = kc ;
49
49
    rr [set] = kr ;
68
68
csd *cs_dmperm (const cs *A, CS_INT seed)
69
69
{
70
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 ;
 
71
        *ps, *rs, nb1, nb2, *p, *q, *cc, *rr, *r, *s, ok ;
72
72
    cs *C ;
73
73
    csd *D, *scc ;
74
74
    /* --- Maximum matching ------------------------------------------------- */
75
 
    if (!CS_CSC (A)) return (NULL) ;            /* check inputs */
 
75
    if (!CS_CSC (A)) return (NULL) ;            /* check inputs */
76
76
    m = A->m ; n = A->n ;
77
 
    D = cs_dalloc (m, n) ;                      /* allocate result */
 
77
    D = cs_dalloc (m, n) ;                      /* allocate result */
78
78
    if (!D) return (NULL) ;
79
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 */
 
80
    jmatch = cs_maxtrans (A, seed) ;            /* max transversal */
 
81
    imatch = jmatch + m ;                       /* imatch = inverse of jmatch */
82
82
    if (!jmatch) return (cs_ddone (D, NULL, jmatch, 0)) ;
83
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*/
 
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
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 */
 
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
95
    cs_free (jmatch) ;
96
96
    /* --- Fine decomposition ----------------------------------------------- */
97
 
    pinv = cs_pinv (p, m) ;         /* pinv=p' */
 
97
    pinv = cs_pinv (p, m) ;         /* pinv=p' */
98
98
    if (!pinv) return (cs_ddone (D, NULL, NULL, 0)) ;
99
99
    C = cs_permute (A, pinv, q, 0) ;/* C=A(p,q) (it will hold A(R2,C2)) */
100
100
    cs_free (pinv) ;
101
101
    if (!C) return (cs_ddone (D, NULL, NULL, 0)) ;
102
102
    Cp = C->p ;
103
 
    nc = cc [3] - cc [2] ;          /* delete cols C0, C1, and C3 from C */
 
103
    nc = cc [3] - cc [2] ;          /* delete cols C0, C1, and C3 from C */
104
104
    if (cc [2] > 0) for (j = cc [2] ; j <= cc [3] ; j++) Cp [j-cc[2]] = Cp [j] ;
105
105
    C->n = nc ;
106
 
    if (rr [2] - rr [1] < m)        /* delete rows R0, R1, and R3 from C */
 
106
    if (rr [2] - rr [1] < m)        /* delete rows R0, R1, and R3 from C */
107
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] ;
 
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
112
    }
113
113
    C->m = nc ;
114
 
    scc = cs_scc (C) ;              /* find strongly connected components of C*/
 
114
    scc = cs_scc (C) ;              /* find strongly connected components of C*/
115
115
    if (!scc) return (cs_ddone (D, C, NULL, 0)) ;
116
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) */
 
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
120
    for (k = 0 ; k < nc ; k++) wj [k] = q [ps [k] + cc [2]] ;
121
121
    for (k = 0 ; k < nc ; k++) q [k + cc [2]] = wj [k] ;
122
122
    for (k = 0 ; k < nc ; k++) wi [k] = p [ps [k] + rr [1]] ;
123
123
    for (k = 0 ; k < nc ; k++) p [k + rr [1]] = wi [k] ;
124
 
    nb2 = 0 ;                       /* create the fine block partitions */
 
124
    nb2 = 0 ;                       /* create the fine block partitions */
125
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) */
 
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
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++ ;
 
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
132
    }
133
133
    if (rr [2] < m)
134
134
    {
135
 
        r [nb2] = rr [2] ;          /* trailing coarse block A ([R3 R0], C3) */
136
 
        s [nb2] = cc [3] ;
137
 
        nb2++ ;
 
135
        r [nb2] = rr [2] ;          /* trailing coarse block A ([R3 R0], C3) */
 
136
        s [nb2] = cc [3] ;
 
137
        nb2++ ;
138
138
    }
139
139
    r [nb2] = m ;
140
140
    s [nb2] = n ;