1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/amd/amd_order.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,200 @@
1.4 +/* ========================================================================= */
1.5 +/* === AMD_order =========================================================== */
1.6 +/* ========================================================================= */
1.7 +
1.8 +/* ------------------------------------------------------------------------- */
1.9 +/* AMD, Copyright (c) Timothy A. Davis, */
1.10 +/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
1.11 +/* email: davis at cise.ufl.edu CISE Department, Univ. of Florida. */
1.12 +/* web: http://www.cise.ufl.edu/research/sparse/amd */
1.13 +/* ------------------------------------------------------------------------- */
1.14 +
1.15 +/* User-callable AMD minimum degree ordering routine. See amd.h for
1.16 + * documentation.
1.17 + */
1.18 +
1.19 +#include "amd_internal.h"
1.20 +
1.21 +/* ========================================================================= */
1.22 +/* === AMD_order =========================================================== */
1.23 +/* ========================================================================= */
1.24 +
1.25 +GLOBAL Int AMD_order
1.26 +(
1.27 + Int n,
1.28 + const Int Ap [ ],
1.29 + const Int Ai [ ],
1.30 + Int P [ ],
1.31 + double Control [ ],
1.32 + double Info [ ]
1.33 +)
1.34 +{
1.35 + Int *Len, *S, nz, i, *Pinv, info, status, *Rp, *Ri, *Cp, *Ci, ok ;
1.36 + size_t nzaat, slen ;
1.37 + double mem = 0 ;
1.38 +
1.39 +#ifndef NDEBUG
1.40 + AMD_debug_init ("amd") ;
1.41 +#endif
1.42 +
1.43 + /* clear the Info array, if it exists */
1.44 + info = Info != (double *) NULL ;
1.45 + if (info)
1.46 + {
1.47 + for (i = 0 ; i < AMD_INFO ; i++)
1.48 + {
1.49 + Info [i] = EMPTY ;
1.50 + }
1.51 + Info [AMD_N] = n ;
1.52 + Info [AMD_STATUS] = AMD_OK ;
1.53 + }
1.54 +
1.55 + /* make sure inputs exist and n is >= 0 */
1.56 + if (Ai == (Int *) NULL || Ap == (Int *) NULL || P == (Int *) NULL || n < 0)
1.57 + {
1.58 + if (info) Info [AMD_STATUS] = AMD_INVALID ;
1.59 + return (AMD_INVALID) ; /* arguments are invalid */
1.60 + }
1.61 +
1.62 + if (n == 0)
1.63 + {
1.64 + return (AMD_OK) ; /* n is 0 so there's nothing to do */
1.65 + }
1.66 +
1.67 + nz = Ap [n] ;
1.68 + if (info)
1.69 + {
1.70 + Info [AMD_NZ] = nz ;
1.71 + }
1.72 + if (nz < 0)
1.73 + {
1.74 + if (info) Info [AMD_STATUS] = AMD_INVALID ;
1.75 + return (AMD_INVALID) ;
1.76 + }
1.77 +
1.78 + /* check if n or nz will cause size_t overflow */
1.79 + if (((size_t) n) >= SIZE_T_MAX / sizeof (Int)
1.80 + || ((size_t) nz) >= SIZE_T_MAX / sizeof (Int))
1.81 + {
1.82 + if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
1.83 + return (AMD_OUT_OF_MEMORY) ; /* problem too large */
1.84 + }
1.85 +
1.86 + /* check the input matrix: AMD_OK, AMD_INVALID, or AMD_OK_BUT_JUMBLED */
1.87 + status = AMD_valid (n, n, Ap, Ai) ;
1.88 +
1.89 + if (status == AMD_INVALID)
1.90 + {
1.91 + if (info) Info [AMD_STATUS] = AMD_INVALID ;
1.92 + return (AMD_INVALID) ; /* matrix is invalid */
1.93 + }
1.94 +
1.95 + /* allocate two size-n integer workspaces */
1.96 + Len = amd_malloc (n * sizeof (Int)) ;
1.97 + Pinv = amd_malloc (n * sizeof (Int)) ;
1.98 + mem += n ;
1.99 + mem += n ;
1.100 + if (!Len || !Pinv)
1.101 + {
1.102 + /* :: out of memory :: */
1.103 + amd_free (Len) ;
1.104 + amd_free (Pinv) ;
1.105 + if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
1.106 + return (AMD_OUT_OF_MEMORY) ;
1.107 + }
1.108 +
1.109 + if (status == AMD_OK_BUT_JUMBLED)
1.110 + {
1.111 + /* sort the input matrix and remove duplicate entries */
1.112 + AMD_DEBUG1 (("Matrix is jumbled\n")) ;
1.113 + Rp = amd_malloc ((n+1) * sizeof (Int)) ;
1.114 + Ri = amd_malloc (MAX (nz,1) * sizeof (Int)) ;
1.115 + mem += (n+1) ;
1.116 + mem += MAX (nz,1) ;
1.117 + if (!Rp || !Ri)
1.118 + {
1.119 + /* :: out of memory :: */
1.120 + amd_free (Rp) ;
1.121 + amd_free (Ri) ;
1.122 + amd_free (Len) ;
1.123 + amd_free (Pinv) ;
1.124 + if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
1.125 + return (AMD_OUT_OF_MEMORY) ;
1.126 + }
1.127 + /* use Len and Pinv as workspace to create R = A' */
1.128 + AMD_preprocess (n, Ap, Ai, Rp, Ri, Len, Pinv) ;
1.129 + Cp = Rp ;
1.130 + Ci = Ri ;
1.131 + }
1.132 + else
1.133 + {
1.134 + /* order the input matrix as-is. No need to compute R = A' first */
1.135 + Rp = NULL ;
1.136 + Ri = NULL ;
1.137 + Cp = (Int *) Ap ;
1.138 + Ci = (Int *) Ai ;
1.139 + }
1.140 +
1.141 + /* --------------------------------------------------------------------- */
1.142 + /* determine the symmetry and count off-diagonal nonzeros in A+A' */
1.143 + /* --------------------------------------------------------------------- */
1.144 +
1.145 + nzaat = AMD_aat (n, Cp, Ci, Len, P, Info) ;
1.146 + AMD_DEBUG1 (("nzaat: %g\n", (double) nzaat)) ;
1.147 + ASSERT ((MAX (nz-n, 0) <= nzaat) && (nzaat <= 2 * (size_t) nz)) ;
1.148 +
1.149 + /* --------------------------------------------------------------------- */
1.150 + /* allocate workspace for matrix, elbow room, and 6 size-n vectors */
1.151 + /* --------------------------------------------------------------------- */
1.152 +
1.153 + S = NULL ;
1.154 + slen = nzaat ; /* space for matrix */
1.155 + ok = ((slen + nzaat/5) >= slen) ; /* check for size_t overflow */
1.156 + slen += nzaat/5 ; /* add elbow room */
1.157 + for (i = 0 ; ok && i < 7 ; i++)
1.158 + {
1.159 + ok = ((slen + n) > slen) ; /* check for size_t overflow */
1.160 + slen += n ; /* size-n elbow room, 6 size-n work */
1.161 + }
1.162 + mem += slen ;
1.163 + ok = ok && (slen < SIZE_T_MAX / sizeof (Int)) ; /* check for overflow */
1.164 + ok = ok && (slen < Int_MAX) ; /* S[i] for Int i must be OK */
1.165 + if (ok)
1.166 + {
1.167 + S = amd_malloc (slen * sizeof (Int)) ;
1.168 + }
1.169 + AMD_DEBUG1 (("slen %g\n", (double) slen)) ;
1.170 + if (!S)
1.171 + {
1.172 + /* :: out of memory :: (or problem too large) */
1.173 + amd_free (Rp) ;
1.174 + amd_free (Ri) ;
1.175 + amd_free (Len) ;
1.176 + amd_free (Pinv) ;
1.177 + if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
1.178 + return (AMD_OUT_OF_MEMORY) ;
1.179 + }
1.180 + if (info)
1.181 + {
1.182 + /* memory usage, in bytes. */
1.183 + Info [AMD_MEMORY] = mem * sizeof (Int) ;
1.184 + }
1.185 +
1.186 + /* --------------------------------------------------------------------- */
1.187 + /* order the matrix */
1.188 + /* --------------------------------------------------------------------- */
1.189 +
1.190 + AMD_1 (n, Cp, Ci, P, Pinv, Len, slen, S, Control, Info) ;
1.191 +
1.192 + /* --------------------------------------------------------------------- */
1.193 + /* free the workspace */
1.194 + /* --------------------------------------------------------------------- */
1.195 +
1.196 + amd_free (Rp) ;
1.197 + amd_free (Ri) ;
1.198 + amd_free (Len) ;
1.199 + amd_free (Pinv) ;
1.200 + amd_free (S) ;
1.201 + if (info) Info [AMD_STATUS] = status ;
1.202 + return (status) ; /* successful ordering */
1.203 +}