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 module amd_order; 16 17 nothrow @nogc extern(C): 18 19 import glob_opts; 20 import amd_internal; 21 import amd; 22 import amd_1; 23 import amd_valid; 24 import amd_preprocess; 25 import amd_aat; 26 import SuiteSparse_config; 27 28 /* ========================================================================= */ 29 /* === AMD_order =========================================================== */ 30 /* ========================================================================= */ 31 32 version (DLONG){ 33 SuiteSparse_long AMD_order 34 ( 35 SuiteSparse_long n, 36 const SuiteSparse_long * Ap, 37 const SuiteSparse_long * Ai, 38 SuiteSparse_long * P, 39 c_float * Control, 40 c_float * Info 41 ) 42 { 43 Int *Len; 44 Int *S; 45 Int nz; 46 Int i; 47 Int *Pinv; 48 Int info; 49 Int status; 50 Int *Rp; 51 Int *Ri; 52 Int *Cp; 53 Int *Ci; 54 Int ok; 55 size_t nzaat; 56 size_t slen; 57 c_float mem = 0 ; 58 59 version(NDEBUG){} 60 else { 61 AMD_debug_init ("amd") ; 62 } 63 64 /* clear the Info array, if it exists */ 65 info = Info != cast(c_float *) NULL ; 66 if (info) 67 { 68 for (i = 0 ; i < AMD_INFO ; i++) 69 { 70 Info [i] = EMPTY ; 71 } 72 Info [AMD_N] = n ; 73 Info [AMD_STATUS] = AMD_OK ; 74 } 75 76 /* make sure inputs exist and n is >= 0 */ 77 if (Ai == cast(Int *) NULL || Ap == cast(Int *) NULL || P == cast(Int *) NULL || n < 0) 78 { 79 if (info) Info [AMD_STATUS] = AMD_INVALID ; 80 return (AMD_INVALID) ; /* arguments are invalid */ 81 } 82 83 if (n == 0) 84 { 85 return (AMD_OK) ; /* n is 0 so there's nothing to do */ 86 } 87 88 nz = Ap [n] ; 89 if (info) 90 { 91 Info [AMD_NZ] = nz ; 92 } 93 if (nz < 0) 94 { 95 if (info) Info [AMD_STATUS] = AMD_INVALID ; 96 return (AMD_INVALID) ; 97 } 98 99 /* check if n or nz will cause size_t overflow */ 100 if ((cast(size_t) n) >= SIZE_T_MAX / (Int.sizeof) 101 || (cast(size_t) nz) >= SIZE_T_MAX / (Int.sizeof)) 102 { 103 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; 104 return (AMD_OUT_OF_MEMORY) ; /* problem too large */ 105 } 106 107 /* check the input matrix: AMD_OK, AMD_INVALID, or AMD_OK_BUT_JUMBLED */ 108 status = AMD_valid (n, n, Ap, Ai) ; 109 110 if (status == AMD_INVALID) 111 { 112 if (info) Info [AMD_STATUS] = AMD_INVALID ; 113 return (AMD_INVALID) ; /* matrix is invalid */ 114 } 115 116 /* allocate two size-n integer workspaces */ 117 Len = cast(Int*)SuiteSparse_malloc (n, (Int.sizeof)) ; 118 Pinv = cast(Int*)SuiteSparse_malloc (n, (Int.sizeof)) ; 119 mem += n ; 120 mem += n ; 121 if (!Len || !Pinv) 122 { 123 /* :: out of memory :: */ 124 SuiteSparse_free (Len) ; 125 SuiteSparse_free (Pinv) ; 126 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; 127 return (AMD_OUT_OF_MEMORY) ; 128 } 129 130 if (status == AMD_OK_BUT_JUMBLED) 131 { 132 /* sort the input matrix and remove duplicate entries */ 133 AMD_DEBUG1 (("Matrix is jumbled\n")) ; 134 Rp = cast(Int*)SuiteSparse_malloc (n+1, (Int.sizeof)) ; 135 Ri = cast(Int*)SuiteSparse_malloc (nz, (Int.sizeof)) ; 136 mem += (n+1) ; 137 mem += MAX (nz,1) ; 138 if (!Rp || !Ri) 139 { 140 /* :: out of memory :: */ 141 SuiteSparse_free (Rp) ; 142 SuiteSparse_free (Ri) ; 143 SuiteSparse_free (Len) ; 144 SuiteSparse_free (Pinv) ; 145 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; 146 return (AMD_OUT_OF_MEMORY) ; 147 } 148 /* use Len and Pinv as workspace to create R = A' */ 149 AMD_preprocess (n, Ap, Ai, Rp, Ri, Len, Pinv) ; 150 Cp = Rp ; 151 Ci = Ri ; 152 } 153 else 154 { 155 /* order the input matrix as-is. No need to compute R = A' first */ 156 Rp = NULL ; 157 Ri = NULL ; 158 Cp = cast(Int *) Ap ; 159 Ci = cast(Int *) Ai ; 160 } 161 162 /* --------------------------------------------------------------------- */ 163 /* determine the symmetry and count off-diagonal nonzeros in A+A' */ 164 /* --------------------------------------------------------------------- */ 165 166 nzaat = AMD_aat (n, Cp, Ci, Len, P, Info) ; 167 AMD_DEBUG1 ("nzaat: %g\n", cast(c_float) nzaat) ; 168 ASSERT ((MAX (nz-n, 0) <= nzaat) && (nzaat <= 2 * cast(size_t) nz)) ; 169 170 /* --------------------------------------------------------------------- */ 171 /* allocate workspace for matrix, elbow room, and 6 size-n vectors */ 172 /* --------------------------------------------------------------------- */ 173 174 S = NULL ; 175 slen = nzaat ; /* space for matrix */ 176 ok = ((slen + nzaat/5) >= slen) ; /* check for size_t overflow */ 177 slen += nzaat/5 ; /* add elbow room */ 178 for (i = 0 ; ok && i < 7 ; i++) 179 { 180 ok = ((slen + n) > slen) ; /* check for size_t overflow */ 181 slen += n ; /* size-n elbow room, 6 size-n work */ 182 } 183 mem += slen ; 184 ok = ok && (slen < SIZE_T_MAX / (Int.sizeof)) ; /* check for overflow */ 185 ok = ok && (slen < Int_MAX) ; /* S[i] for Int i must be OK */ 186 if (ok) 187 { 188 S = cast(Int*)SuiteSparse_malloc (slen, (Int.sizeof)) ; 189 } 190 AMD_DEBUG1 ("slen %g\n", cast(c_float) slen) ; 191 if (!S) 192 { 193 /* :: out of memory :: (or problem too large) */ 194 SuiteSparse_free (Rp) ; 195 SuiteSparse_free (Ri) ; 196 SuiteSparse_free (Len) ; 197 SuiteSparse_free (Pinv) ; 198 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; 199 return (AMD_OUT_OF_MEMORY) ; 200 } 201 if (info) 202 { 203 /* memory usage, in bytes. */ 204 Info [AMD_MEMORY] = mem * (Int.sizeof) ; 205 } 206 207 /* --------------------------------------------------------------------- */ 208 /* order the matrix */ 209 /* --------------------------------------------------------------------- */ 210 211 AMD_1 (n, Cp, Ci, P, Pinv, Len, cast(Int)slen, S, Control, Info) ; 212 213 /* --------------------------------------------------------------------- */ 214 /* free the workspace */ 215 /* --------------------------------------------------------------------- */ 216 217 SuiteSparse_free (Rp) ; 218 SuiteSparse_free (Ri) ; 219 SuiteSparse_free (Len) ; 220 SuiteSparse_free (Pinv) ; 221 SuiteSparse_free (S) ; 222 if (info) Info [AMD_STATUS] = status ; 223 return (status) ; /* successful ordering */ 224 } 225 } 226 else { 227 int AMD_order /* returns AMD_OK, AMD_OK_BUT_JUMBLED, 228 * AMD_INVALID, or AMD_OUT_OF_MEMORY */ 229 ( 230 int n, /* A is n-by-n. n must be >= 0. */ 231 const int * Ap, /* column pointers for A, of size n+1 */ 232 const int * Ai, /* row indices of A, of size nz = Ap [n] */ 233 int * P, /* output permutation, of size n */ 234 c_float * Control, /* input Control settings, of size AMD_CONTROL */ 235 c_float * Info /* output Info statistics, of size AMD_INFO */ 236 ) 237 { 238 Int *Len; 239 Int *S; 240 Int nz; 241 Int i; 242 Int *Pinv; 243 Int info; 244 Int status; 245 Int *Rp; 246 Int *Ri; 247 Int *Cp; 248 Int *Ci; 249 Int ok; 250 size_t nzaat; 251 size_t slen; 252 c_float mem = 0 ; 253 254 version(NDEBUG){} 255 else { 256 AMD_debug_init ("amd") ; 257 } 258 259 /* clear the Info array, if it exists */ 260 info = Info != cast(c_float *) NULL ; 261 if (info) 262 { 263 for (i = 0 ; i < AMD_INFO ; i++) 264 { 265 Info [i] = EMPTY ; 266 } 267 Info [AMD_N] = n ; 268 Info [AMD_STATUS] = AMD_OK ; 269 } 270 271 /* make sure inputs exist and n is >= 0 */ 272 if (Ai == cast(Int *) NULL || Ap == cast(Int *) NULL || P == cast(Int *) NULL || n < 0) 273 { 274 if (info) Info [AMD_STATUS] = AMD_INVALID ; 275 return (AMD_INVALID) ; /* arguments are invalid */ 276 } 277 278 if (n == 0) 279 { 280 return (AMD_OK) ; /* n is 0 so there's nothing to do */ 281 } 282 283 nz = Ap [n] ; 284 if (info) 285 { 286 Info [AMD_NZ] = nz ; 287 } 288 if (nz < 0) 289 { 290 if (info) Info [AMD_STATUS] = AMD_INVALID ; 291 return (AMD_INVALID) ; 292 } 293 294 /* check if n or nz will cause size_t overflow */ 295 if ((cast(size_t) n) >= SIZE_T_MAX / (Int.sizeof) 296 || (cast(size_t) nz) >= SIZE_T_MAX / (Int.sizeof)) 297 { 298 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; 299 return (AMD_OUT_OF_MEMORY) ; /* problem too large */ 300 } 301 302 /* check the input matrix: AMD_OK, AMD_INVALID, or AMD_OK_BUT_JUMBLED */ 303 status = AMD_valid (n, n, Ap, Ai) ; 304 305 if (status == AMD_INVALID) 306 { 307 if (info) Info [AMD_STATUS] = AMD_INVALID ; 308 return (AMD_INVALID) ; /* matrix is invalid */ 309 } 310 311 /* allocate two size-n integer workspaces */ 312 Len = cast(Int*)SuiteSparse_malloc (n, (Int.sizeof)) ; 313 Pinv = cast(Int*)SuiteSparse_malloc (n, (Int.sizeof)) ; 314 mem += n ; 315 mem += n ; 316 if (!Len || !Pinv) 317 { 318 /* :: out of memory :: */ 319 SuiteSparse_free (Len) ; 320 SuiteSparse_free (Pinv) ; 321 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; 322 return (AMD_OUT_OF_MEMORY) ; 323 } 324 325 if (status == AMD_OK_BUT_JUMBLED) 326 { 327 /* sort the input matrix and remove duplicate entries */ 328 AMD_DEBUG1 (("Matrix is jumbled\n")) ; 329 Rp = cast(Int*)SuiteSparse_malloc (n+1, (Int.sizeof)) ; 330 Ri = cast(Int*)SuiteSparse_malloc (nz, (Int.sizeof)) ; 331 mem += (n+1) ; 332 mem += MAX (nz,1) ; 333 if (!Rp || !Ri) 334 { 335 /* :: out of memory :: */ 336 SuiteSparse_free (Rp) ; 337 SuiteSparse_free (Ri) ; 338 SuiteSparse_free (Len) ; 339 SuiteSparse_free (Pinv) ; 340 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; 341 return (AMD_OUT_OF_MEMORY) ; 342 } 343 /* use Len and Pinv as workspace to create R = A' */ 344 AMD_preprocess (n, Ap, Ai, Rp, Ri, Len, Pinv) ; 345 Cp = Rp ; 346 Ci = Ri ; 347 } 348 else 349 { 350 /* order the input matrix as-is. No need to compute R = A' first */ 351 Rp = NULL ; 352 Ri = NULL ; 353 Cp = cast(Int *) Ap ; 354 Ci = cast(Int *) Ai ; 355 } 356 357 /* --------------------------------------------------------------------- */ 358 /* determine the symmetry and count off-diagonal nonzeros in A+A' */ 359 /* --------------------------------------------------------------------- */ 360 361 nzaat = AMD_aat (n, Cp, Ci, Len, P, Info) ; 362 AMD_DEBUG1 ("nzaat: %g\n", cast(c_float) nzaat) ; 363 ASSERT ((MAX (nz-n, 0) <= nzaat) && (nzaat <= 2 * cast(size_t) nz)) ; 364 365 /* --------------------------------------------------------------------- */ 366 /* allocate workspace for matrix, elbow room, and 6 size-n vectors */ 367 /* --------------------------------------------------------------------- */ 368 369 S = NULL ; 370 slen = nzaat ; /* space for matrix */ 371 ok = ((slen + nzaat/5) >= slen) ; /* check for size_t overflow */ 372 slen += nzaat/5 ; /* add elbow room */ 373 for (i = 0 ; ok && i < 7 ; i++) 374 { 375 ok = ((slen + n) > slen) ; /* check for size_t overflow */ 376 slen += n ; /* size-n elbow room, 6 size-n work */ 377 } 378 mem += slen ; 379 ok = ok && (slen < SIZE_T_MAX / (Int.sizeof)) ; /* check for overflow */ 380 ok = ok && (slen < Int_MAX) ; /* S[i] for Int i must be OK */ 381 if (ok) 382 { 383 S = cast(Int*)SuiteSparse_malloc (slen, (Int.sizeof)) ; 384 } 385 AMD_DEBUG1 ("slen %g\n", cast(c_float) slen) ; 386 if (!S) 387 { 388 /* :: out of memory :: (or problem too large) */ 389 SuiteSparse_free (Rp) ; 390 SuiteSparse_free (Ri) ; 391 SuiteSparse_free (Len) ; 392 SuiteSparse_free (Pinv) ; 393 if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ; 394 return (AMD_OUT_OF_MEMORY) ; 395 } 396 if (info) 397 { 398 /* memory usage, in bytes. */ 399 Info [AMD_MEMORY] = mem * (Int.sizeof) ; 400 } 401 402 /* --------------------------------------------------------------------- */ 403 /* order the matrix */ 404 /* --------------------------------------------------------------------- */ 405 406 AMD_1 (n, Cp, Ci, P, Pinv, Len, cast(Int)slen, S, Control, Info) ; 407 408 /* --------------------------------------------------------------------- */ 409 /* free the workspace */ 410 /* --------------------------------------------------------------------- */ 411 412 SuiteSparse_free (Rp) ; 413 SuiteSparse_free (Ri) ; 414 SuiteSparse_free (Len) ; 415 SuiteSparse_free (Pinv) ; 416 SuiteSparse_free (S) ; 417 if (info) Info [AMD_STATUS] = status ; 418 return (status) ; /* successful ordering */ 419 } 420 } // DLONG 421