# source:glpk-cmake/src/amd/amd_aat.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: 5.8 KB
Line
1/* ========================================================================= */
2/* === AMD_aat ============================================================= */
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_aat:  compute the symmetry of the pattern of A, and count the number of
13 * nonzeros each column of A+A' (excluding the diagonal).  Assumes the input
14 * matrix has no errors, with sorted columns and no duplicates
15 * (AMD_valid (n, n, Ap, Ai) must be AMD_OK, but this condition is not
16 * checked).
17 */
18
19#include "amd_internal.h"
20
21GLOBAL size_t AMD_aat   /* returns nz in A+A' */
22(
23    Int n,
24    const Int Ap [ ],
25    const Int Ai [ ],
26    Int Len [ ],        /* Len [j]: length of column j of A+A', excl diagonal*/
27    Int Tp [ ],         /* workspace of size n */
28    double Info [ ]
29)
30{
31    Int p1, p2, p, i, j, pj, pj2, k, nzdiag, nzboth, nz ;
32    double sym ;
33    size_t nzaat ;
34
35#ifndef NDEBUG
36    AMD_debug_init ("AMD AAT") ;
37    for (k = 0 ; k < n ; k++) Tp [k] = EMPTY ;
38    ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
39#endif
40
41    if (Info != (double *) NULL)
42    {
43        /* clear the Info array, if it exists */
44        for (i = 0 ; i < AMD_INFO ; i++)
45        {
46            Info [i] = EMPTY ;
47        }
48        Info [AMD_STATUS] = AMD_OK ;
49    }
50
51    for (k = 0 ; k < n ; k++)
52    {
53        Len [k] = 0 ;
54    }
55
56    nzdiag = 0 ;
57    nzboth = 0 ;
58    nz = Ap [n] ;
59
60    for (k = 0 ; k < n ; k++)
61    {
62        p1 = Ap [k] ;
63        p2 = Ap [k+1] ;
64        AMD_DEBUG2 (("\nAAT Column: "ID" p1: "ID" p2: "ID"\n", k, p1, p2)) ;
65
66        /* construct A+A' */
67        for (p = p1 ; p < p2 ; )
68        {
69            /* scan the upper triangular part of A */
70            j = Ai [p] ;
71            if (j < k)
72            {
73                /* entry A (j,k) is in the strictly upper triangular part,
74                 * add both A (j,k) and A (k,j) to the matrix A+A' */
75                Len [j]++ ;
76                Len [k]++ ;
77                AMD_DEBUG3 (("    upper ("ID","ID") ("ID","ID")\n", j,k, k,j));
78                p++ ;
79            }
80            else if (j == k)
81            {
82                /* skip the diagonal */
83                p++ ;
84                nzdiag++ ;
85                break ;
86            }
87            else /* j > k */
88            {
89                /* first entry below the diagonal */
90                break ;
91            }
92            /* scan lower triangular part of A, in column j until reaching
93             * row k.  Start where last scan left off. */
94            ASSERT (Tp [j] != EMPTY) ;
95            ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
96            pj2 = Ap [j+1] ;
97            for (pj = Tp [j] ; pj < pj2 ; )
98            {
99                i = Ai [pj] ;
100                if (i < k)
101                {
102                    /* A (i,j) is only in the lower part, not in upper.
103                     * add both A (i,j) and A (j,i) to the matrix A+A' */
104                    Len [i]++ ;
105                    Len [j]++ ;
106                    AMD_DEBUG3 (("    lower ("ID","ID") ("ID","ID")\n",
107                        i,j, j,i)) ;
108                    pj++ ;
109                }
110                else if (i == k)
111                {
112                    /* entry A (k,j) in lower part and A (j,k) in upper */
113                    pj++ ;
114                    nzboth++ ;
115                    break ;
116                }
117                else /* i > k */
118                {
119                    /* consider this entry later, when k advances to i */
120                    break ;
121                }
122            }
123            Tp [j] = pj ;
124        }
125        /* Tp [k] points to the entry just below the diagonal in column k */
126        Tp [k] = p ;
127    }
128
129    /* clean up, for remaining mismatched entries */
130    for (j = 0 ; j < n ; j++)
131    {
132        for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
133        {
134            i = Ai [pj] ;
135            /* A (i,j) is only in the lower part, not in upper.
136             * add both A (i,j) and A (j,i) to the matrix A+A' */
137            Len [i]++ ;
138            Len [j]++ ;
139            AMD_DEBUG3 (("    lower cleanup ("ID","ID") ("ID","ID")\n",
140                i,j, j,i)) ;
141        }
142    }
143
144    /* --------------------------------------------------------------------- */
145    /* compute the symmetry of the nonzero pattern of A */
146    /* --------------------------------------------------------------------- */
147
148    /* Given a matrix A, the symmetry of A is:
149     *  B = tril (spones (A), -1) + triu (spones (A), 1) ;
150     *  sym = nnz (B & B') / nnz (B) ;
151     *  or 1 if nnz (B) is zero.
152     */
153
154    if (nz == nzdiag)
155    {
156        sym = 1 ;
157    }
158    else
159    {
160        sym = (2 * (double) nzboth) / ((double) (nz - nzdiag)) ;
161    }
162
163    nzaat = 0 ;
164    for (k = 0 ; k < n ; k++)
165    {
166        nzaat += Len [k] ;
167    }
168
169    AMD_DEBUG1 (("AMD nz in A+A', excluding diagonal (nzaat) = %g\n",
170        (double) nzaat)) ;
171    AMD_DEBUG1 (("   nzboth: "ID" nz: "ID" nzdiag: "ID" symmetry: %g\n",
172                nzboth, nz, nzdiag, sym)) ;
173
174    if (Info != (double *) NULL)
175    {
176        Info [AMD_STATUS] = AMD_OK ;
177        Info [AMD_N] = n ;
178        Info [AMD_NZ] = nz ;
179        Info [AMD_SYMMETRY] = sym ;         /* symmetry of pattern of A */
180        Info [AMD_NZDIAG] = nzdiag ;        /* nonzeros on diagonal of A */
181        Info [AMD_NZ_A_PLUS_AT] = nzaat ;   /* nonzeros in A+A' */
182    }
183
184    return (nzaat) ;
185}
Note: See TracBrowser for help on using the repository browser.