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