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

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

     2 /* === AMD_1 =============================================================== */

     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 /* AMD_1: Construct A+A' for a sparse matrix A and perform the AMD ordering.

    13  *

    14  * The n-by-n sparse matrix A can be unsymmetric.  It is stored in MATLAB-style

    15  * compressed-column form, with sorted row indices in each column, and no

    16  * duplicate entries.  Diagonal entries may be present, but they are ignored.

    17  * Row indices of column j of A are stored in Ai [Ap [j] ... Ap [j+1]-1].

    18  * Ap [0] must be zero, and nz = Ap [n] is the number of entries in A.  The

    19  * size of the matrix, n, must be greater than or equal to zero.

    20  *

    21  * This routine must be preceded by a call to AMD_aat, which computes the

    22  * number of entries in each row/column in A+A', excluding the diagonal.

    23  * Len [j], on input, is the number of entries in row/column j of A+A'.  This

    24  * routine constructs the matrix A+A' and then calls AMD_2.  No error checking

    25  * is performed (this was done in AMD_valid).

    26  */

    28 #include "amd_internal.h"

    30 GLOBAL void AMD_1

    31 (

    32     Int n,              /* n > 0 */

    33     const Int Ap [ ],   /* input of size n+1, not modified */

    34     const Int Ai [ ],   /* input of size nz = Ap [n], not modified */

    35     Int P [ ],          /* size n output permutation */

    36     Int Pinv [ ],       /* size n output inverse permutation */

    37     Int Len [ ],        /* size n input, undefined on output */

    38     Int slen,           /* slen >= sum (Len [0..n-1]) + 7n,

    39                          * ideally slen = 1.2 * sum (Len) + 8n */

    40     Int S [ ],          /* size slen workspace */

    41     double Control [ ], /* input array of size AMD_CONTROL */

    42     double Info [ ]     /* output array of size AMD_INFO */

    43 )

    44 {

    45     Int i, j, k, p, pfree, iwlen, pj, p1, p2, pj2, *Iw, *Pe, *Nv, *Head,

    46         *Elen, *Degree, *s, *W, *Sp, *Tp ;

    48     /* --------------------------------------------------------------------- */

    49     /* construct the matrix for AMD_2 */

    50     /* --------------------------------------------------------------------- */

    52     ASSERT (n > 0) ;

    54     iwlen = slen - 6*n ;

    55     s = S ;

    56     Pe = s ;        s += n ;

    57     Nv = s ;        s += n ;

    58     Head = s ;      s += n ;

    59     Elen = s ;      s += n ;

    60     Degree = s ;    s += n ;

    61     W = s ;         s += n ;

    62     Iw = s ;        s += iwlen ;

    64     ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;

    66     /* construct the pointers for A+A' */

    67     Sp = Nv ;                   /* use Nv and W as workspace for Sp and Tp [ */

    68     Tp = W ;

    69     pfree = 0 ;

    70     for (j = 0 ; j < n ; j++)

    71     {

    72         Pe [j] = pfree ;

    73         Sp [j] = pfree ;

    74         pfree += Len [j] ;

    75     }

    77     /* Note that this restriction on iwlen is slightly more restrictive than

    78      * what is strictly required in AMD_2.  AMD_2 can operate with no elbow

    79      * room at all, but it will be very slow.  For better performance, at

    80      * least size-n elbow room is enforced. */

    81     ASSERT (iwlen >= pfree + n) ;

    83 #ifndef NDEBUG

    84     for (p = 0 ; p < iwlen ; p++) Iw [p] = EMPTY ;

    85 #endif

    87     for (k = 0 ; k < n ; k++)

    88     {

    89         AMD_DEBUG1 (("Construct row/column k= "ID" of A+A'\n", k))  ;

    90         p1 = Ap [k] ;

    91         p2 = Ap [k+1] ;

    93         /* construct A+A' */

    94         for (p = p1 ; p < p2 ; )

    95         {

    96             /* scan the upper triangular part of A */

    97             j = Ai [p] ;

    98             ASSERT (j >= 0 && j < n) ;

    99             if (j < k)

   100             {

   101                 /* entry A (j,k) in the strictly upper triangular part */

   102                 ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;

   103                 ASSERT (Sp [k] < (k == n-1 ? pfree : Pe [k+1])) ;

   104                 Iw [Sp [j]++] = k ;

   105                 Iw [Sp [k]++] = j ;

   106                 p++ ;

   107             }

   108             else if (j == k)

   109             {

   110                 /* skip the diagonal */

   111                 p++ ;

   112                 break ;

   113             }

   114             else /* j > k */

   115             {

   116                 /* first entry below the diagonal */

   117                 break ;

   118             }

   119             /* scan lower triangular part of A, in column j until reaching

   120              * row k.  Start where last scan left off. */

   121             ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;

   122             pj2 = Ap [j+1] ;

   123             for (pj = Tp [j] ; pj < pj2 ; )

   124             {

   125                 i = Ai [pj] ;

   126                 ASSERT (i >= 0 && i < n) ;

   127                 if (i < k)

   128                 {

   129                     /* A (i,j) is only in the lower part, not in upper */

   130                     ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;

   131                     ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;

   132                     Iw [Sp [i]++] = j ;

   133                     Iw [Sp [j]++] = i ;

   134                     pj++ ;

   135                 }

   136                 else if (i == k)

   137                 {

   138                     /* entry A (k,j) in lower part and A (j,k) in upper */

   139                     pj++ ;

   140                     break ;

   141                 }

   142                 else /* i > k */

   143                 {

   144                     /* consider this entry later, when k advances to i */

   145                     break ;

   146                 }

   147             }

   148             Tp [j] = pj ;

   149         }

   150         Tp [k] = p ;

   151     }

   153     /* clean up, for remaining mismatched entries */

   154     for (j = 0 ; j < n ; j++)

   155     {

   156         for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)

   157         {

   158             i = Ai [pj] ;

   159             ASSERT (i >= 0 && i < n) ;

   160             /* A (i,j) is only in the lower part, not in upper */

   161             ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;

   162             ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;

   163             Iw [Sp [i]++] = j ;

   164             Iw [Sp [j]++] = i ;

   165         }

   166     }

   168 #ifndef NDEBUG

   169     for (j = 0 ; j < n-1 ; j++) ASSERT (Sp [j] == Pe [j+1]) ;

   170     ASSERT (Sp [n-1] == pfree) ;

   171 #endif

   173     /* Tp and Sp no longer needed ] */

   175     /* --------------------------------------------------------------------- */

   176     /* order the matrix */

   177     /* --------------------------------------------------------------------- */

   179     AMD_2 (n, Pe, Iw, Len, iwlen, pfree,

   180         Nv, Pinv, P, Head, Elen, Degree, W, Control, Info) ;

   181 }

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