lemon-project-template-glpk
view deps/glpk/src/amd/amd_aat.c @ 11:4fc6ad2fb8a6
Test GLPK in src/main.cc
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 21:43:29 +0100 |
| parents | |
| children |
line source
1 /* ========================================================================= */
2 /* === AMD_aat ============================================================= */
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_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 */
19 #include "amd_internal.h"
21 GLOBAL 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 ;
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
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 }
51 for (k = 0 ; k < n ; k++)
52 {
53 Len [k] = 0 ;
54 }
56 nzdiag = 0 ;
57 nzboth = 0 ;
58 nz = Ap [n] ;
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)) ;
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 }
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 }
144 /* --------------------------------------------------------------------- */
145 /* compute the symmetry of the nonzero pattern of A */
146 /* --------------------------------------------------------------------- */
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 */
154 if (nz == nzdiag)
155 {
156 sym = 1 ;
157 }
158 else
159 {
160 sym = (2 * (double) nzboth) / ((double) (nz - nzdiag)) ;
161 }
163 nzaat = 0 ;
164 for (k = 0 ; k < n ; k++)
165 {
166 nzaat += Len [k] ;
167 }
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)) ;
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 }
184 return (nzaat) ;
185 }