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