1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/amd/amd_aat.c	Sun Nov 06 20:59:10 2011 +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 +}