|
1 /* ========================================================================= */ |
|
2 /* === AMD_order =========================================================== */ |
|
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 /* User-callable AMD minimum degree ordering routine. See amd.h for |
|
13 * documentation. |
|
14 */ |
|
15 |
|
16 #include "amd_internal.h" |
|
17 |
|
18 /* ========================================================================= */ |
|
19 /* === AMD_order =========================================================== */ |
|
20 /* ========================================================================= */ |
|
21 |
|
22 GLOBAL Int AMD_order |
|
23 ( |
|
24 Int n, |
|
25 const Int Ap [ ], |
|
26 const Int Ai [ ], |
|
27 Int P [ ], |
|
28 double Control [ ], |
|
29 double Info [ ] |
|
30 ) |
|
31 { |
|
32 Int *Len, *S, nz, i, *Pinv, info, status, *Rp, *Ri, *Cp, *Ci, ok ; |
|
33 size_t nzaat, slen ; |
|
34 double mem = 0 ; |
|
35 |
|
36 #ifndef NDEBUG |
|
37 AMD_debug_init ("amd") ; |
|
38 #endif |
|
39 |
|
40 /* clear the Info array, if it exists */ |
|
41 info = Info != (double *) NULL ; |
|
42 if (info) |
|
43 { |
|
44 for (i = 0 ; i < AMD_INFO ; i++) |
|
45 { |
|
46 Info [i] = EMPTY ; |
|
47 } |
|
48 Info [AMD_N] = n ; |
|
49 Info [AMD_STATUS] = AMD_OK ; |
|
50 } |
|
51 |
|
52 /* make sure inputs exist and n is >= 0 */ |
|
53 if (Ai == (Int *) NULL || Ap == (Int *) NULL || P == (Int *) NULL || n < 0) |
|
54 { |
|
55 if (info) Info [AMD_STATUS] = AMD_INVALID ; |
|
56 return (AMD_INVALID) ; /* arguments are invalid */ |
|
57 } |
|
58 |
|
59 if (n == 0) |
|
60 { |
|
61 return (AMD_OK) ; /* n is 0 so there's nothing to do */ |
|
62 } |
|
63 |
|
64 nz = Ap [n] ; |
|
65 if (info) |
|
66 { |
|
67 Info [AMD_NZ] = nz ; |
|
68 } |
|
69 if (nz < 0) |
|
70 { |
|
71 if (info) Info [AMD_STATUS] = AMD_INVALID ; |
|
72 return (AMD_INVALID) ; |
|
73 } |
|
74 |
|
75 /* check if n or nz will cause size_t overflow */ |
|
76 if (((size_t) n) >= SIZE_T_MAX / sizeof (Int) |
|
77 || ((size_t) nz) >= SIZE_T_MAX / sizeof (Int)) |
|
78 { |
|
79 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; |
|
80 return (AMD_OUT_OF_MEMORY) ; /* problem too large */ |
|
81 } |
|
82 |
|
83 /* check the input matrix: AMD_OK, AMD_INVALID, or AMD_OK_BUT_JUMBLED */ |
|
84 status = AMD_valid (n, n, Ap, Ai) ; |
|
85 |
|
86 if (status == AMD_INVALID) |
|
87 { |
|
88 if (info) Info [AMD_STATUS] = AMD_INVALID ; |
|
89 return (AMD_INVALID) ; /* matrix is invalid */ |
|
90 } |
|
91 |
|
92 /* allocate two size-n integer workspaces */ |
|
93 Len = amd_malloc (n * sizeof (Int)) ; |
|
94 Pinv = amd_malloc (n * sizeof (Int)) ; |
|
95 mem += n ; |
|
96 mem += n ; |
|
97 if (!Len || !Pinv) |
|
98 { |
|
99 /* :: out of memory :: */ |
|
100 amd_free (Len) ; |
|
101 amd_free (Pinv) ; |
|
102 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; |
|
103 return (AMD_OUT_OF_MEMORY) ; |
|
104 } |
|
105 |
|
106 if (status == AMD_OK_BUT_JUMBLED) |
|
107 { |
|
108 /* sort the input matrix and remove duplicate entries */ |
|
109 AMD_DEBUG1 (("Matrix is jumbled\n")) ; |
|
110 Rp = amd_malloc ((n+1) * sizeof (Int)) ; |
|
111 Ri = amd_malloc (MAX (nz,1) * sizeof (Int)) ; |
|
112 mem += (n+1) ; |
|
113 mem += MAX (nz,1) ; |
|
114 if (!Rp || !Ri) |
|
115 { |
|
116 /* :: out of memory :: */ |
|
117 amd_free (Rp) ; |
|
118 amd_free (Ri) ; |
|
119 amd_free (Len) ; |
|
120 amd_free (Pinv) ; |
|
121 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; |
|
122 return (AMD_OUT_OF_MEMORY) ; |
|
123 } |
|
124 /* use Len and Pinv as workspace to create R = A' */ |
|
125 AMD_preprocess (n, Ap, Ai, Rp, Ri, Len, Pinv) ; |
|
126 Cp = Rp ; |
|
127 Ci = Ri ; |
|
128 } |
|
129 else |
|
130 { |
|
131 /* order the input matrix as-is. No need to compute R = A' first */ |
|
132 Rp = NULL ; |
|
133 Ri = NULL ; |
|
134 Cp = (Int *) Ap ; |
|
135 Ci = (Int *) Ai ; |
|
136 } |
|
137 |
|
138 /* --------------------------------------------------------------------- */ |
|
139 /* determine the symmetry and count off-diagonal nonzeros in A+A' */ |
|
140 /* --------------------------------------------------------------------- */ |
|
141 |
|
142 nzaat = AMD_aat (n, Cp, Ci, Len, P, Info) ; |
|
143 AMD_DEBUG1 (("nzaat: %g\n", (double) nzaat)) ; |
|
144 ASSERT ((MAX (nz-n, 0) <= nzaat) && (nzaat <= 2 * (size_t) nz)) ; |
|
145 |
|
146 /* --------------------------------------------------------------------- */ |
|
147 /* allocate workspace for matrix, elbow room, and 6 size-n vectors */ |
|
148 /* --------------------------------------------------------------------- */ |
|
149 |
|
150 S = NULL ; |
|
151 slen = nzaat ; /* space for matrix */ |
|
152 ok = ((slen + nzaat/5) >= slen) ; /* check for size_t overflow */ |
|
153 slen += nzaat/5 ; /* add elbow room */ |
|
154 for (i = 0 ; ok && i < 7 ; i++) |
|
155 { |
|
156 ok = ((slen + n) > slen) ; /* check for size_t overflow */ |
|
157 slen += n ; /* size-n elbow room, 6 size-n work */ |
|
158 } |
|
159 mem += slen ; |
|
160 ok = ok && (slen < SIZE_T_MAX / sizeof (Int)) ; /* check for overflow */ |
|
161 ok = ok && (slen < Int_MAX) ; /* S[i] for Int i must be OK */ |
|
162 if (ok) |
|
163 { |
|
164 S = amd_malloc (slen * sizeof (Int)) ; |
|
165 } |
|
166 AMD_DEBUG1 (("slen %g\n", (double) slen)) ; |
|
167 if (!S) |
|
168 { |
|
169 /* :: out of memory :: (or problem too large) */ |
|
170 amd_free (Rp) ; |
|
171 amd_free (Ri) ; |
|
172 amd_free (Len) ; |
|
173 amd_free (Pinv) ; |
|
174 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; |
|
175 return (AMD_OUT_OF_MEMORY) ; |
|
176 } |
|
177 if (info) |
|
178 { |
|
179 /* memory usage, in bytes. */ |
|
180 Info [AMD_MEMORY] = mem * sizeof (Int) ; |
|
181 } |
|
182 |
|
183 /* --------------------------------------------------------------------- */ |
|
184 /* order the matrix */ |
|
185 /* --------------------------------------------------------------------- */ |
|
186 |
|
187 AMD_1 (n, Cp, Ci, P, Pinv, Len, slen, S, Control, Info) ; |
|
188 |
|
189 /* --------------------------------------------------------------------- */ |
|
190 /* free the workspace */ |
|
191 /* --------------------------------------------------------------------- */ |
|
192 |
|
193 amd_free (Rp) ; |
|
194 amd_free (Ri) ; |
|
195 amd_free (Len) ; |
|
196 amd_free (Pinv) ; |
|
197 amd_free (S) ; |
|
198 if (info) Info [AMD_STATUS] = status ; |
|
199 return (status) ; /* successful ordering */ |
|
200 } |