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

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

• Generated files and doc/notes are removed
File size: 6.5 KB
Line
1/* ========================================================================= */
2/* === AMD_1 =============================================================== */
3/* ========================================================================= */
4
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/* ------------------------------------------------------------------------- */
11
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  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 */
27
28#include "amd_internal.h"
29
30GLOBAL 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 ;
47
48    /* --------------------------------------------------------------------- */
49    /* construct the matrix for AMD_2 */
50    /* --------------------------------------------------------------------- */
51
52    ASSERT (n > 0) ;
53
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 ;
63
64    ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
65
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    }
76
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) ;
82
83#ifndef NDEBUG
84    for (p = 0 ; p < iwlen ; p++) Iw [p] = EMPTY ;
85#endif
86
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] ;
92
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    }
152
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    }
167
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
172
173    /* Tp and Sp no longer needed ] */
174
175    /* --------------------------------------------------------------------- */
176    /* order the matrix */
177    /* --------------------------------------------------------------------- */
178
179    AMD_2 (n, Pe, Iw, Len, iwlen, pfree,
180        Nv, Pinv, P, Head, Elen, Degree, W, Control, Info) ;
181}
Note: See TracBrowser for help on using the repository browser.