amesos_amd_l_aat.c

Go to the documentation of this file.

/* ========================================================================= */
/* === AMD_aat ============================================================= */
/* ========================================================================= */

/* ------------------------------------------------------------------------- */
/* AMD, Copyright (c) Timothy A. Davis,              */
/* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
/* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
/* web: http://www.cise.ufl.edu/research/sparse/amd                          */
/* ------------------------------------------------------------------------- */

/* AMD_aat:  compute the symmetry of the pattern of A, and count the number of
 * nonzeros each column of A+A' (excluding the diagonal).  Assumes the input
 * matrix has no errors, with sorted columns and no duplicates
 * (AMD_valid (n, n, Ap, Ai) must be AMD_OK, but this condition is not
 * checked).
 */

/* This file should make the long int version of AMD */
#define DLONG 1

#include "amesos_amd_internal.h"

GLOBAL size_t AMD_aat /* returns nz in A+A' */
(
    Int n,
    const Int Ap [ ],
    const Int Ai [ ],
    Int Len [ ],  /* Len [j]: length of column j of A+A', excl diagonal*/
    Int Tp [ ],   /* workspace of size n */
    double Info [ ]
)
{
    Int p1, p2, p, i, j, pj, pj2, k, nzdiag, nzboth, nz ;
    double sym ;
    size_t nzaat ;

#ifndef NDEBUG
    AMD_debug_init ("AMD AAT") ;
    for (k = 0 ; k < n ; k++) Tp [k] = EMPTY ;
    ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
#endif

    if (Info != (double *) NULL)
    {
  /* clear the Info array, if it exists */
  for (i = 0 ; i < AMD_INFO ; i++)
  {
      Info [i] = EMPTY ;
  }
  Info [AMD_STATUS] = AMD_OK ;
    }

    for (k = 0 ; k < n ; k++)
    {
  Len [k] = 0 ;
    }

    nzdiag = 0 ;
    nzboth = 0 ;
    nz = Ap [n] ;

    for (k = 0 ; k < n ; k++)
    {
  p1 = Ap [k] ;
  p2 = Ap [k+1] ;
  AMD_DEBUG2 (("\nAAT Column: "ID" p1: "ID" p2: "ID"\n", k, p1, p2)) ;

  /* construct A+A' */
  for (p = p1 ; p < p2 ; )
  {
      /* scan the upper triangular part of A */
      j = Ai [p] ;
      if (j < k)
      {
    /* entry A (j,k) is in the strictly upper triangular part,
     * add both A (j,k) and A (k,j) to the matrix A+A' */
    Len [j]++ ;
    Len [k]++ ;
    AMD_DEBUG3 (("    upper ("ID","ID") ("ID","ID")\n", j,k, k,j));
    p++ ;
      }
      else if (j == k)
      {
    /* skip the diagonal */
    p++ ;
    nzdiag++ ;
    break ;
      }
      else /* j > k */
      {
    /* first entry below the diagonal */
    break ;
      }
      /* scan lower triangular part of A, in column j until reaching
       * row k.  Start where last scan left off. */
      ASSERT (Tp [j] != EMPTY) ;
      ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
      pj2 = Ap [j+1] ;
      for (pj = Tp [j] ; pj < pj2 ; )
      {
    i = Ai [pj] ;
    if (i < k)
    {
        /* A (i,j) is only in the lower part, not in upper.
         * add both A (i,j) and A (j,i) to the matrix A+A' */
        Len [i]++ ;
        Len [j]++ ;
        AMD_DEBUG3 (("    lower ("ID","ID") ("ID","ID")\n",
      i,j, j,i)) ;
        pj++ ;
    }
    else if (i == k)
    {
        /* entry A (k,j) in lower part and A (j,k) in upper */
        pj++ ;
        nzboth++ ;
        break ;
    }
    else /* i > k */
    {
        /* consider this entry later, when k advances to i */
        break ;
    }
      }
      Tp [j] = pj ;
  }
  /* Tp [k] points to the entry just below the diagonal in column k */
  Tp [k] = p ;
    }

    /* clean up, for remaining mismatched entries */
    for (j = 0 ; j < n ; j++)
    {
  for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
  {
      i = Ai [pj] ;
      /* A (i,j) is only in the lower part, not in upper.
       * add both A (i,j) and A (j,i) to the matrix A+A' */
      Len [i]++ ;
      Len [j]++ ;
      AMD_DEBUG3 (("    lower cleanup ("ID","ID") ("ID","ID")\n",
    i,j, j,i)) ;
  }
    }

    /* --------------------------------------------------------------------- */
    /* compute the symmetry of the nonzero pattern of A */
    /* --------------------------------------------------------------------- */

    /* Given a matrix A, the symmetry of A is:
     *  B = tril (spones (A), -1) + triu (spones (A), 1) ;
     *  sym = nnz (B & B') / nnz (B) ;
     *  or 1 if nnz (B) is zero.
     */

    if (nz == nzdiag)
    {
  sym = 1 ;
    }
    else
    {
  sym = (2 * (double) nzboth) / ((double) (nz - nzdiag)) ;
    }

    nzaat = 0 ;
    for (k = 0 ; k < n ; k++)
    {
  nzaat += Len [k] ;
    }

    AMD_DEBUG1 (("AMD nz in A+A', excluding diagonal (nzaat) = %g\n",
  (double) nzaat)) ;
    AMD_DEBUG1 (("   nzboth: "ID" nz: "ID" nzdiag: "ID" symmetry: %g\n",
    nzboth, nz, nzdiag, sym)) ;

    if (Info != (double *) NULL)
    {
  Info [AMD_STATUS] = AMD_OK ;
  Info [AMD_N] = n ;
  Info [AMD_NZ] = nz ;
  Info [AMD_SYMMETRY] = sym ;     /* symmetry of pattern of A */
  Info [AMD_NZDIAG] = nzdiag ;      /* nonzeros on diagonal of A */
  Info [AMD_NZ_A_PLUS_AT] = nzaat ;   /* nonzeros in A+A' */
    }

    return (nzaat) ;
}