lemon-project-template-glpk

annotate deps/glpk/src/amd/amd_order.c @ 9:33de93886c88

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