glpk-cmake: src/amd/amd_order.c@c445c931472f

     1 /* ========================================================================= */

     2 /* === AMD_order =========================================================== */

     3 /* ========================================================================= */

     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 /* ------------------------------------------------------------------------- */

    12 /* User-callable AMD minimum degree ordering routine.  See amd.h for

    13  * documentation.

    14  */

    16 #include "amd_internal.h"

    18 /* ========================================================================= */

    19 /* === AMD_order =========================================================== */

    20 /* ========================================================================= */

    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 ;

    36 #ifndef NDEBUG

    37     AMD_debug_init ("amd") ;

    38 #endif

    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     }

    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     }

    59     if (n == 0)

    60     {

    61         return (AMD_OK) ;           /* n is 0 so there's nothing to do */

    62     }

    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     }

    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     }

    83     /* check the input matrix:  AMD_OK, AMD_INVALID, or AMD_OK_BUT_JUMBLED */

    84     status = AMD_valid (n, n, Ap, Ai) ;

    86     if (status == AMD_INVALID)

    87     {

    88         if (info) Info [AMD_STATUS] = AMD_INVALID ;

    89         return (AMD_INVALID) ;      /* matrix is invalid */

    90     }

    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     }

   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     }

   138     /* --------------------------------------------------------------------- */

   139     /* determine the symmetry and count off-diagonal nonzeros in A+A' */

   140     /* --------------------------------------------------------------------- */

   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)) ;

   146     /* --------------------------------------------------------------------- */

   147     /* allocate workspace for matrix, elbow room, and 6 size-n vectors */

   148     /* --------------------------------------------------------------------- */

   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     }

   183     /* --------------------------------------------------------------------- */

   184     /* order the matrix */

   185     /* --------------------------------------------------------------------- */

   187     AMD_1 (n, Cp, Ci, P, Pinv, Len, slen, S, Control, Info) ;

   189     /* --------------------------------------------------------------------- */

   190     /* free the workspace */

   191     /* --------------------------------------------------------------------- */

   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 }

author	Alpar Juttner <alpar@cs.elte.hu>
	Mon, 06 Dec 2010 13:09:21 +0100
changeset 1	c445c931472f
permissions	-rw-r--r--