alpar@1: /* ========================================================================= */
alpar@1: /* === AMD_order =========================================================== */
alpar@1: /* ========================================================================= */
alpar@1: 
alpar@1: /* ------------------------------------------------------------------------- */
alpar@1: /* AMD, Copyright (c) Timothy A. Davis,                                      */
alpar@1: /* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
alpar@1: /* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
alpar@1: /* web: http://www.cise.ufl.edu/research/sparse/amd                          */
alpar@1: /* ------------------------------------------------------------------------- */
alpar@1: 
alpar@1: /* User-callable AMD minimum degree ordering routine.  See amd.h for
alpar@1:  * documentation.
alpar@1:  */
alpar@1: 
alpar@1: #include "amd_internal.h"
alpar@1: 
alpar@1: /* ========================================================================= */
alpar@1: /* === AMD_order =========================================================== */
alpar@1: /* ========================================================================= */
alpar@1: 
alpar@1: GLOBAL Int AMD_order
alpar@1: (
alpar@1:     Int n,
alpar@1:     const Int Ap [ ],
alpar@1:     const Int Ai [ ],
alpar@1:     Int P [ ],
alpar@1:     double Control [ ],
alpar@1:     double Info [ ]
alpar@1: )
alpar@1: {
alpar@1:     Int *Len, *S, nz, i, *Pinv, info, status, *Rp, *Ri, *Cp, *Ci, ok ;
alpar@1:     size_t nzaat, slen ;
alpar@1:     double mem = 0 ;
alpar@1: 
alpar@1: #ifndef NDEBUG
alpar@1:     AMD_debug_init ("amd") ;
alpar@1: #endif
alpar@1: 
alpar@1:     /* clear the Info array, if it exists */
alpar@1:     info = Info != (double *) NULL ;
alpar@1:     if (info)
alpar@1:     {
alpar@1:         for (i = 0 ; i < AMD_INFO ; i++)
alpar@1:         {
alpar@1:             Info [i] = EMPTY ;
alpar@1:         }
alpar@1:         Info [AMD_N] = n ;
alpar@1:         Info [AMD_STATUS] = AMD_OK ;
alpar@1:     }
alpar@1: 
alpar@1:     /* make sure inputs exist and n is >= 0 */
alpar@1:     if (Ai == (Int *) NULL || Ap == (Int *) NULL || P == (Int *) NULL || n < 0)
alpar@1:     {
alpar@1:         if (info) Info [AMD_STATUS] = AMD_INVALID ;
alpar@1:         return (AMD_INVALID) ;      /* arguments are invalid */
alpar@1:     }
alpar@1: 
alpar@1:     if (n == 0)
alpar@1:     {
alpar@1:         return (AMD_OK) ;           /* n is 0 so there's nothing to do */
alpar@1:     }
alpar@1: 
alpar@1:     nz = Ap [n] ;
alpar@1:     if (info)
alpar@1:     {
alpar@1:         Info [AMD_NZ] = nz ;
alpar@1:     }
alpar@1:     if (nz < 0)
alpar@1:     {
alpar@1:         if (info) Info [AMD_STATUS] = AMD_INVALID ;
alpar@1:         return (AMD_INVALID) ;
alpar@1:     }
alpar@1: 
alpar@1:     /* check if n or nz will cause size_t overflow */
alpar@1:     if (((size_t) n) >= SIZE_T_MAX / sizeof (Int)
alpar@1:      || ((size_t) nz) >= SIZE_T_MAX / sizeof (Int))
alpar@1:     {
alpar@1:         if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
alpar@1:         return (AMD_OUT_OF_MEMORY) ;        /* problem too large */
alpar@1:     }
alpar@1: 
alpar@1:     /* check the input matrix:  AMD_OK, AMD_INVALID, or AMD_OK_BUT_JUMBLED */
alpar@1:     status = AMD_valid (n, n, Ap, Ai) ;
alpar@1: 
alpar@1:     if (status == AMD_INVALID)
alpar@1:     {
alpar@1:         if (info) Info [AMD_STATUS] = AMD_INVALID ;
alpar@1:         return (AMD_INVALID) ;      /* matrix is invalid */
alpar@1:     }
alpar@1: 
alpar@1:     /* allocate two size-n integer workspaces */
alpar@1:     Len = amd_malloc (n * sizeof (Int)) ;
alpar@1:     Pinv = amd_malloc (n * sizeof (Int)) ;
alpar@1:     mem += n ;
alpar@1:     mem += n ;
alpar@1:     if (!Len || !Pinv)
alpar@1:     {
alpar@1:         /* :: out of memory :: */
alpar@1:         amd_free (Len) ;
alpar@1:         amd_free (Pinv) ;
alpar@1:         if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
alpar@1:         return (AMD_OUT_OF_MEMORY) ;
alpar@1:     }
alpar@1: 
alpar@1:     if (status == AMD_OK_BUT_JUMBLED)
alpar@1:     {
alpar@1:         /* sort the input matrix and remove duplicate entries */
alpar@1:         AMD_DEBUG1 (("Matrix is jumbled\n")) ;
alpar@1:         Rp = amd_malloc ((n+1) * sizeof (Int)) ;
alpar@1:         Ri = amd_malloc (MAX (nz,1) * sizeof (Int)) ;
alpar@1:         mem += (n+1) ;
alpar@1:         mem += MAX (nz,1) ;
alpar@1:         if (!Rp || !Ri)
alpar@1:         {
alpar@1:             /* :: out of memory :: */
alpar@1:             amd_free (Rp) ;
alpar@1:             amd_free (Ri) ;
alpar@1:             amd_free (Len) ;
alpar@1:             amd_free (Pinv) ;
alpar@1:             if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
alpar@1:             return (AMD_OUT_OF_MEMORY) ;
alpar@1:         }
alpar@1:         /* use Len and Pinv as workspace to create R = A' */
alpar@1:         AMD_preprocess (n, Ap, Ai, Rp, Ri, Len, Pinv) ;
alpar@1:         Cp = Rp ;
alpar@1:         Ci = Ri ;
alpar@1:     }
alpar@1:     else
alpar@1:     {
alpar@1:         /* order the input matrix as-is.  No need to compute R = A' first */
alpar@1:         Rp = NULL ;
alpar@1:         Ri = NULL ;
alpar@1:         Cp = (Int *) Ap ;
alpar@1:         Ci = (Int *) Ai ;
alpar@1:     }
alpar@1: 
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1:     /* determine the symmetry and count off-diagonal nonzeros in A+A' */
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1: 
alpar@1:     nzaat = AMD_aat (n, Cp, Ci, Len, P, Info) ;
alpar@1:     AMD_DEBUG1 (("nzaat: %g\n", (double) nzaat)) ;
alpar@1:     ASSERT ((MAX (nz-n, 0) <= nzaat) && (nzaat <= 2 * (size_t) nz)) ;
alpar@1: 
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1:     /* allocate workspace for matrix, elbow room, and 6 size-n vectors */
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1: 
alpar@1:     S = NULL ;
alpar@1:     slen = nzaat ;                      /* space for matrix */
alpar@1:     ok = ((slen + nzaat/5) >= slen) ;   /* check for size_t overflow */
alpar@1:     slen += nzaat/5 ;                   /* add elbow room */
alpar@1:     for (i = 0 ; ok && i < 7 ; i++)
alpar@1:     {
alpar@1:         ok = ((slen + n) > slen) ;      /* check for size_t overflow */
alpar@1:         slen += n ;                     /* size-n elbow room, 6 size-n work */
alpar@1:     }
alpar@1:     mem += slen ;
alpar@1:     ok = ok && (slen < SIZE_T_MAX / sizeof (Int)) ; /* check for overflow */
alpar@1:     ok = ok && (slen < Int_MAX) ;       /* S[i] for Int i must be OK */
alpar@1:     if (ok)
alpar@1:     {
alpar@1:         S = amd_malloc (slen * sizeof (Int)) ;
alpar@1:     }
alpar@1:     AMD_DEBUG1 (("slen %g\n", (double) slen)) ;
alpar@1:     if (!S)
alpar@1:     {
alpar@1:         /* :: out of memory :: (or problem too large) */
alpar@1:         amd_free (Rp) ;
alpar@1:         amd_free (Ri) ;
alpar@1:         amd_free (Len) ;
alpar@1:         amd_free (Pinv) ;
alpar@1:         if (info) Info [AMD_STATUS] = AMD_OUT_OF_MEMORY ;
alpar@1:         return (AMD_OUT_OF_MEMORY) ;
alpar@1:     }
alpar@1:     if (info)
alpar@1:     {
alpar@1:         /* memory usage, in bytes. */
alpar@1:         Info [AMD_MEMORY] = mem * sizeof (Int) ;
alpar@1:     }
alpar@1: 
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1:     /* order the matrix */
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1: 
alpar@1:     AMD_1 (n, Cp, Ci, P, Pinv, Len, slen, S, Control, Info) ;
alpar@1: 
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1:     /* free the workspace */
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1: 
alpar@1:     amd_free (Rp) ;
alpar@1:     amd_free (Ri) ;
alpar@1:     amd_free (Len) ;
alpar@1:     amd_free (Pinv) ;
alpar@1:     amd_free (S) ;
alpar@1:     if (info) Info [AMD_STATUS] = status ;
alpar@1:     return (status) ;       /* successful ordering */
alpar@1: }