[9] | 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 | } |
---|