1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/amd/amd_aat.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,185 @@
1.4 +/* ========================================================================= */
1.5 +/* === AMD_aat ============================================================= */
1.6 +/* ========================================================================= */
1.7 +
1.8 +/* ------------------------------------------------------------------------- */
1.9 +/* AMD, Copyright (c) Timothy A. Davis, */
1.10 +/* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */
1.11 +/* email: davis at cise.ufl.edu CISE Department, Univ. of Florida. */
1.12 +/* web: http://www.cise.ufl.edu/research/sparse/amd */
1.13 +/* ------------------------------------------------------------------------- */
1.14 +
1.15 +/* AMD_aat: compute the symmetry of the pattern of A, and count the number of
1.16 + * nonzeros each column of A+A' (excluding the diagonal). Assumes the input
1.17 + * matrix has no errors, with sorted columns and no duplicates
1.18 + * (AMD_valid (n, n, Ap, Ai) must be AMD_OK, but this condition is not
1.19 + * checked).
1.20 + */
1.21 +
1.22 +#include "amd_internal.h"
1.23 +
1.24 +GLOBAL size_t AMD_aat /* returns nz in A+A' */
1.25 +(
1.26 + Int n,
1.27 + const Int Ap [ ],
1.28 + const Int Ai [ ],
1.29 + Int Len [ ], /* Len [j]: length of column j of A+A', excl diagonal*/
1.30 + Int Tp [ ], /* workspace of size n */
1.31 + double Info [ ]
1.32 +)
1.33 +{
1.34 + Int p1, p2, p, i, j, pj, pj2, k, nzdiag, nzboth, nz ;
1.35 + double sym ;
1.36 + size_t nzaat ;
1.37 +
1.38 +#ifndef NDEBUG
1.39 + AMD_debug_init ("AMD AAT") ;
1.40 + for (k = 0 ; k < n ; k++) Tp [k] = EMPTY ;
1.41 + ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
1.42 +#endif
1.43 +
1.44 + if (Info != (double *) NULL)
1.45 + {
1.46 + /* clear the Info array, if it exists */
1.47 + for (i = 0 ; i < AMD_INFO ; i++)
1.48 + {
1.49 + Info [i] = EMPTY ;
1.50 + }
1.51 + Info [AMD_STATUS] = AMD_OK ;
1.52 + }
1.53 +
1.54 + for (k = 0 ; k < n ; k++)
1.55 + {
1.56 + Len [k] = 0 ;
1.57 + }
1.58 +
1.59 + nzdiag = 0 ;
1.60 + nzboth = 0 ;
1.61 + nz = Ap [n] ;
1.62 +
1.63 + for (k = 0 ; k < n ; k++)
1.64 + {
1.65 + p1 = Ap [k] ;
1.66 + p2 = Ap [k+1] ;
1.67 + AMD_DEBUG2 (("\nAAT Column: "ID" p1: "ID" p2: "ID"\n", k, p1, p2)) ;
1.68 +
1.69 + /* construct A+A' */
1.70 + for (p = p1 ; p < p2 ; )
1.71 + {
1.72 + /* scan the upper triangular part of A */
1.73 + j = Ai [p] ;
1.74 + if (j < k)
1.75 + {
1.76 + /* entry A (j,k) is in the strictly upper triangular part,
1.77 + * add both A (j,k) and A (k,j) to the matrix A+A' */
1.78 + Len [j]++ ;
1.79 + Len [k]++ ;
1.80 + AMD_DEBUG3 ((" upper ("ID","ID") ("ID","ID")\n", j,k, k,j));
1.81 + p++ ;
1.82 + }
1.83 + else if (j == k)
1.84 + {
1.85 + /* skip the diagonal */
1.86 + p++ ;
1.87 + nzdiag++ ;
1.88 + break ;
1.89 + }
1.90 + else /* j > k */
1.91 + {
1.92 + /* first entry below the diagonal */
1.93 + break ;
1.94 + }
1.95 + /* scan lower triangular part of A, in column j until reaching
1.96 + * row k. Start where last scan left off. */
1.97 + ASSERT (Tp [j] != EMPTY) ;
1.98 + ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
1.99 + pj2 = Ap [j+1] ;
1.100 + for (pj = Tp [j] ; pj < pj2 ; )
1.101 + {
1.102 + i = Ai [pj] ;
1.103 + if (i < k)
1.104 + {
1.105 + /* A (i,j) is only in the lower part, not in upper.
1.106 + * add both A (i,j) and A (j,i) to the matrix A+A' */
1.107 + Len [i]++ ;
1.108 + Len [j]++ ;
1.109 + AMD_DEBUG3 ((" lower ("ID","ID") ("ID","ID")\n",
1.110 + i,j, j,i)) ;
1.111 + pj++ ;
1.112 + }
1.113 + else if (i == k)
1.114 + {
1.115 + /* entry A (k,j) in lower part and A (j,k) in upper */
1.116 + pj++ ;
1.117 + nzboth++ ;
1.118 + break ;
1.119 + }
1.120 + else /* i > k */
1.121 + {
1.122 + /* consider this entry later, when k advances to i */
1.123 + break ;
1.124 + }
1.125 + }
1.126 + Tp [j] = pj ;
1.127 + }
1.128 + /* Tp [k] points to the entry just below the diagonal in column k */
1.129 + Tp [k] = p ;
1.130 + }
1.131 +
1.132 + /* clean up, for remaining mismatched entries */
1.133 + for (j = 0 ; j < n ; j++)
1.134 + {
1.135 + for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
1.136 + {
1.137 + i = Ai [pj] ;
1.138 + /* A (i,j) is only in the lower part, not in upper.
1.139 + * add both A (i,j) and A (j,i) to the matrix A+A' */
1.140 + Len [i]++ ;
1.141 + Len [j]++ ;
1.142 + AMD_DEBUG3 ((" lower cleanup ("ID","ID") ("ID","ID")\n",
1.143 + i,j, j,i)) ;
1.144 + }
1.145 + }
1.146 +
1.147 + /* --------------------------------------------------------------------- */
1.148 + /* compute the symmetry of the nonzero pattern of A */
1.149 + /* --------------------------------------------------------------------- */
1.150 +
1.151 + /* Given a matrix A, the symmetry of A is:
1.152 + * B = tril (spones (A), -1) + triu (spones (A), 1) ;
1.153 + * sym = nnz (B & B') / nnz (B) ;
1.154 + * or 1 if nnz (B) is zero.
1.155 + */
1.156 +
1.157 + if (nz == nzdiag)
1.158 + {
1.159 + sym = 1 ;
1.160 + }
1.161 + else
1.162 + {
1.163 + sym = (2 * (double) nzboth) / ((double) (nz - nzdiag)) ;
1.164 + }
1.165 +
1.166 + nzaat = 0 ;
1.167 + for (k = 0 ; k < n ; k++)
1.168 + {
1.169 + nzaat += Len [k] ;
1.170 + }
1.171 +
1.172 + AMD_DEBUG1 (("AMD nz in A+A', excluding diagonal (nzaat) = %g\n",
1.173 + (double) nzaat)) ;
1.174 + AMD_DEBUG1 ((" nzboth: "ID" nz: "ID" nzdiag: "ID" symmetry: %g\n",
1.175 + nzboth, nz, nzdiag, sym)) ;
1.176 +
1.177 + if (Info != (double *) NULL)
1.178 + {
1.179 + Info [AMD_STATUS] = AMD_OK ;
1.180 + Info [AMD_N] = n ;
1.181 + Info [AMD_NZ] = nz ;
1.182 + Info [AMD_SYMMETRY] = sym ; /* symmetry of pattern of A */
1.183 + Info [AMD_NZDIAG] = nzdiag ; /* nonzeros on diagonal of A */
1.184 + Info [AMD_NZ_A_PLUS_AT] = nzaat ; /* nonzeros in A+A' */
1.185 + }
1.186 +
1.187 + return (nzaat) ;
1.188 +}