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