lemon-project-template-glpk

view deps/glpk/src/amd/amd_1.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
line source
1 /* ========================================================================= */
2 /* === AMD_1 =============================================================== */
3 /* ========================================================================= */
4
5 /* ------------------------------------------------------------------------- */
6 /* AMD, Copyright (c) Timothy A. Davis, */
7 /* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
8 /* email: davis at cise.ufl.edu CISE Department, Univ. of Florida. */
9 /* web: http://www.cise.ufl.edu/research/sparse/amd */
10 /* ------------------------------------------------------------------------- */
11
12 /* AMD_1: Construct A+A' for a sparse matrix A and perform the AMD ordering.
13 *
14 * The n-by-n sparse matrix A can be unsymmetric. It is stored in MATLAB-style
15 * compressed-column form, with sorted row indices in each column, and no
16 * duplicate entries. Diagonal entries may be present, but they are ignored.
17 * Row indices of column j of A are stored in Ai [Ap [j] ... Ap [j+1]-1].
18 * Ap [0] must be zero, and nz = Ap [n] is the number of entries in A. The
19 * size of the matrix, n, must be greater than or equal to zero.
20 *
21 * This routine must be preceded by a call to AMD_aat, which computes the
22 * number of entries in each row/column in A+A', excluding the diagonal.
23 * Len [j], on input, is the number of entries in row/column j of A+A'. This
24 * routine constructs the matrix A+A' and then calls AMD_2. No error checking
25 * is performed (this was done in AMD_valid).
26 */
27
28 #include "amd_internal.h"
29
30 GLOBAL void AMD_1
31 (
32 Int n, /* n > 0 */
33 const Int Ap [ ], /* input of size n+1, not modified */
34 const Int Ai [ ], /* input of size nz = Ap [n], not modified */
35 Int P [ ], /* size n output permutation */
36 Int Pinv [ ], /* size n output inverse permutation */
37 Int Len [ ], /* size n input, undefined on output */
38 Int slen, /* slen >= sum (Len [0..n-1]) + 7n,
39 * ideally slen = 1.2 * sum (Len) + 8n */
40 Int S [ ], /* size slen workspace */
41 double Control [ ], /* input array of size AMD_CONTROL */
42 double Info [ ] /* output array of size AMD_INFO */
43 )
44 {
45 Int i, j, k, p, pfree, iwlen, pj, p1, p2, pj2, *Iw, *Pe, *Nv, *Head,
46 *Elen, *Degree, *s, *W, *Sp, *Tp ;
47
48 /* --------------------------------------------------------------------- */
49 /* construct the matrix for AMD_2 */
50 /* --------------------------------------------------------------------- */
51
52 ASSERT (n > 0) ;
53
54 iwlen = slen - 6*n ;
55 s = S ;
56 Pe = s ; s += n ;
57 Nv = s ; s += n ;
58 Head = s ; s += n ;
59 Elen = s ; s += n ;
60 Degree = s ; s += n ;
61 W = s ; s += n ;
62 Iw = s ; s += iwlen ;
63
64 ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
65
66 /* construct the pointers for A+A' */
67 Sp = Nv ; /* use Nv and W as workspace for Sp and Tp [ */
68 Tp = W ;
69 pfree = 0 ;
70 for (j = 0 ; j < n ; j++)
71 {
72 Pe [j] = pfree ;
73 Sp [j] = pfree ;
74 pfree += Len [j] ;
75 }
76
77 /* Note that this restriction on iwlen is slightly more restrictive than
78 * what is strictly required in AMD_2. AMD_2 can operate with no elbow
79 * room at all, but it will be very slow. For better performance, at
80 * least size-n elbow room is enforced. */
81 ASSERT (iwlen >= pfree + n) ;
82
83 #ifndef NDEBUG
84 for (p = 0 ; p < iwlen ; p++) Iw [p] = EMPTY ;
85 #endif
86
87 for (k = 0 ; k < n ; k++)
88 {
89 AMD_DEBUG1 (("Construct row/column k= "ID" of A+A'\n", k)) ;
90 p1 = Ap [k] ;
91 p2 = Ap [k+1] ;
92
93 /* construct A+A' */
94 for (p = p1 ; p < p2 ; )
95 {
96 /* scan the upper triangular part of A */
97 j = Ai [p] ;
98 ASSERT (j >= 0 && j < n) ;
99 if (j < k)
100 {
101 /* entry A (j,k) in the strictly upper triangular part */
102 ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
103 ASSERT (Sp [k] < (k == n-1 ? pfree : Pe [k+1])) ;
104 Iw [Sp [j]++] = k ;
105 Iw [Sp [k]++] = j ;
106 p++ ;
107 }
108 else if (j == k)
109 {
110 /* skip the diagonal */
111 p++ ;
112 break ;
113 }
114 else /* j > k */
115 {
116 /* first entry below the diagonal */
117 break ;
118 }
119 /* scan lower triangular part of A, in column j until reaching
120 * row k. Start where last scan left off. */
121 ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
122 pj2 = Ap [j+1] ;
123 for (pj = Tp [j] ; pj < pj2 ; )
124 {
125 i = Ai [pj] ;
126 ASSERT (i >= 0 && i < n) ;
127 if (i < k)
128 {
129 /* A (i,j) is only in the lower part, not in upper */
130 ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;
131 ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
132 Iw [Sp [i]++] = j ;
133 Iw [Sp [j]++] = i ;
134 pj++ ;
135 }
136 else if (i == k)
137 {
138 /* entry A (k,j) in lower part and A (j,k) in upper */
139 pj++ ;
140 break ;
141 }
142 else /* i > k */
143 {
144 /* consider this entry later, when k advances to i */
145 break ;
146 }
147 }
148 Tp [j] = pj ;
149 }
150 Tp [k] = p ;
151 }
152
153 /* clean up, for remaining mismatched entries */
154 for (j = 0 ; j < n ; j++)
155 {
156 for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
157 {
158 i = Ai [pj] ;
159 ASSERT (i >= 0 && i < n) ;
160 /* A (i,j) is only in the lower part, not in upper */
161 ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;
162 ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
163 Iw [Sp [i]++] = j ;
164 Iw [Sp [j]++] = i ;
165 }
166 }
167
168 #ifndef NDEBUG
169 for (j = 0 ; j < n-1 ; j++) ASSERT (Sp [j] == Pe [j+1]) ;
170 ASSERT (Sp [n-1] == pfree) ;
171 #endif
172
173 /* Tp and Sp no longer needed ] */
174
175 /* --------------------------------------------------------------------- */
176 /* order the matrix */
177 /* --------------------------------------------------------------------- */
178
179 AMD_2 (n, Pe, Iw, Len, iwlen, pfree,
180 Nv, Pinv, P, Head, Elen, Degree, W, Control, Info) ;
181 }