Mercurial > octave-nkf

comparison liboctave/UMFPACK/AMD/Source/amd_aat.c @ 5164:57077d0ddc8e

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