alpar@1: /* ========================================================================= */
alpar@1: /* === AMD_aat ============================================================= */
alpar@1: /* ========================================================================= */
alpar@1: 
alpar@1: /* ------------------------------------------------------------------------- */
alpar@1: /* AMD, Copyright (c) Timothy A. Davis,                                      */
alpar@1: /* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
alpar@1: /* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
alpar@1: /* web: http://www.cise.ufl.edu/research/sparse/amd                          */
alpar@1: /* ------------------------------------------------------------------------- */
alpar@1: 
alpar@1: /* AMD_aat:  compute the symmetry of the pattern of A, and count the number of
alpar@1:  * nonzeros each column of A+A' (excluding the diagonal).  Assumes the input
alpar@1:  * matrix has no errors, with sorted columns and no duplicates
alpar@1:  * (AMD_valid (n, n, Ap, Ai) must be AMD_OK, but this condition is not
alpar@1:  * checked).
alpar@1:  */
alpar@1: 
alpar@1: #include "amd_internal.h"
alpar@1: 
alpar@1: GLOBAL size_t AMD_aat   /* returns nz in A+A' */
alpar@1: (
alpar@1:     Int n,
alpar@1:     const Int Ap [ ],
alpar@1:     const Int Ai [ ],
alpar@1:     Int Len [ ],        /* Len [j]: length of column j of A+A', excl diagonal*/
alpar@1:     Int Tp [ ],         /* workspace of size n */
alpar@1:     double Info [ ]
alpar@1: )
alpar@1: {
alpar@1:     Int p1, p2, p, i, j, pj, pj2, k, nzdiag, nzboth, nz ;
alpar@1:     double sym ;
alpar@1:     size_t nzaat ;
alpar@1: 
alpar@1: #ifndef NDEBUG
alpar@1:     AMD_debug_init ("AMD AAT") ;
alpar@1:     for (k = 0 ; k < n ; k++) Tp [k] = EMPTY ;
alpar@1:     ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
alpar@1: #endif
alpar@1: 
alpar@1:     if (Info != (double *) NULL)
alpar@1:     {
alpar@1:         /* clear the Info array, if it exists */
alpar@1:         for (i = 0 ; i < AMD_INFO ; i++)
alpar@1:         {
alpar@1:             Info [i] = EMPTY ;
alpar@1:         }
alpar@1:         Info [AMD_STATUS] = AMD_OK ;
alpar@1:     }
alpar@1: 
alpar@1:     for (k = 0 ; k < n ; k++)
alpar@1:     {
alpar@1:         Len [k] = 0 ;
alpar@1:     }
alpar@1: 
alpar@1:     nzdiag = 0 ;
alpar@1:     nzboth = 0 ;
alpar@1:     nz = Ap [n] ;
alpar@1: 
alpar@1:     for (k = 0 ; k < n ; k++)
alpar@1:     {
alpar@1:         p1 = Ap [k] ;
alpar@1:         p2 = Ap [k+1] ;
alpar@1:         AMD_DEBUG2 (("\nAAT Column: "ID" p1: "ID" p2: "ID"\n", k, p1, p2)) ;
alpar@1: 
alpar@1:         /* construct A+A' */
alpar@1:         for (p = p1 ; p < p2 ; )
alpar@1:         {
alpar@1:             /* scan the upper triangular part of A */
alpar@1:             j = Ai [p] ;
alpar@1:             if (j < k)
alpar@1:             {
alpar@1:                 /* entry A (j,k) is in the strictly upper triangular part,
alpar@1:                  * add both A (j,k) and A (k,j) to the matrix A+A' */
alpar@1:                 Len [j]++ ;
alpar@1:                 Len [k]++ ;
alpar@1:                 AMD_DEBUG3 (("    upper ("ID","ID") ("ID","ID")\n", j,k, k,j));
alpar@1:                 p++ ;
alpar@1:             }
alpar@1:             else if (j == k)
alpar@1:             {
alpar@1:                 /* skip the diagonal */
alpar@1:                 p++ ;
alpar@1:                 nzdiag++ ;
alpar@1:                 break ;
alpar@1:             }
alpar@1:             else /* j > k */
alpar@1:             {
alpar@1:                 /* first entry below the diagonal */
alpar@1:                 break ;
alpar@1:             }
alpar@1:             /* scan lower triangular part of A, in column j until reaching
alpar@1:              * row k.  Start where last scan left off. */
alpar@1:             ASSERT (Tp [j] != EMPTY) ;
alpar@1:             ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
alpar@1:             pj2 = Ap [j+1] ;
alpar@1:             for (pj = Tp [j] ; pj < pj2 ; )
alpar@1:             {
alpar@1:                 i = Ai [pj] ;
alpar@1:                 if (i < k)
alpar@1:                 {
alpar@1:                     /* A (i,j) is only in the lower part, not in upper.
alpar@1:                      * add both A (i,j) and A (j,i) to the matrix A+A' */
alpar@1:                     Len [i]++ ;
alpar@1:                     Len [j]++ ;
alpar@1:                     AMD_DEBUG3 (("    lower ("ID","ID") ("ID","ID")\n",
alpar@1:                         i,j, j,i)) ;
alpar@1:                     pj++ ;
alpar@1:                 }
alpar@1:                 else if (i == k)
alpar@1:                 {
alpar@1:                     /* entry A (k,j) in lower part and A (j,k) in upper */
alpar@1:                     pj++ ;
alpar@1:                     nzboth++ ;
alpar@1:                     break ;
alpar@1:                 }
alpar@1:                 else /* i > k */
alpar@1:                 {
alpar@1:                     /* consider this entry later, when k advances to i */
alpar@1:                     break ;
alpar@1:                 }
alpar@1:             }
alpar@1:             Tp [j] = pj ;
alpar@1:         }
alpar@1:         /* Tp [k] points to the entry just below the diagonal in column k */
alpar@1:         Tp [k] = p ;
alpar@1:     }
alpar@1: 
alpar@1:     /* clean up, for remaining mismatched entries */
alpar@1:     for (j = 0 ; j < n ; j++)
alpar@1:     {
alpar@1:         for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
alpar@1:         {
alpar@1:             i = Ai [pj] ;
alpar@1:             /* A (i,j) is only in the lower part, not in upper.
alpar@1:              * add both A (i,j) and A (j,i) to the matrix A+A' */
alpar@1:             Len [i]++ ;
alpar@1:             Len [j]++ ;
alpar@1:             AMD_DEBUG3 (("    lower cleanup ("ID","ID") ("ID","ID")\n",
alpar@1:                 i,j, j,i)) ;
alpar@1:         }
alpar@1:     }
alpar@1: 
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1:     /* compute the symmetry of the nonzero pattern of A */
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1: 
alpar@1:     /* Given a matrix A, the symmetry of A is:
alpar@1:      *  B = tril (spones (A), -1) + triu (spones (A), 1) ;
alpar@1:      *  sym = nnz (B & B') / nnz (B) ;
alpar@1:      *  or 1 if nnz (B) is zero.
alpar@1:      */
alpar@1: 
alpar@1:     if (nz == nzdiag)
alpar@1:     {
alpar@1:         sym = 1 ;
alpar@1:     }
alpar@1:     else
alpar@1:     {
alpar@1:         sym = (2 * (double) nzboth) / ((double) (nz - nzdiag)) ;
alpar@1:     }
alpar@1: 
alpar@1:     nzaat = 0 ;
alpar@1:     for (k = 0 ; k < n ; k++)
alpar@1:     {
alpar@1:         nzaat += Len [k] ;
alpar@1:     }
alpar@1: 
alpar@1:     AMD_DEBUG1 (("AMD nz in A+A', excluding diagonal (nzaat) = %g\n",
alpar@1:         (double) nzaat)) ;
alpar@1:     AMD_DEBUG1 (("   nzboth: "ID" nz: "ID" nzdiag: "ID" symmetry: %g\n",
alpar@1:                 nzboth, nz, nzdiag, sym)) ;
alpar@1: 
alpar@1:     if (Info != (double *) NULL)
alpar@1:     {
alpar@1:         Info [AMD_STATUS] = AMD_OK ;
alpar@1:         Info [AMD_N] = n ;
alpar@1:         Info [AMD_NZ] = nz ;
alpar@1:         Info [AMD_SYMMETRY] = sym ;         /* symmetry of pattern of A */
alpar@1:         Info [AMD_NZDIAG] = nzdiag ;        /* nonzeros on diagonal of A */
alpar@1:         Info [AMD_NZ_A_PLUS_AT] = nzaat ;   /* nonzeros in A+A' */
alpar@1:     }
alpar@1: 
alpar@1:     return (nzaat) ;
alpar@1: }