|
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 } |