alpar@1: /* ========================================================================= */
alpar@1: /* === AMD_1 =============================================================== */
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_1: Construct A+A' for a sparse matrix A and perform the AMD ordering.
alpar@1:  *
alpar@1:  * The n-by-n sparse matrix A can be unsymmetric.  It is stored in MATLAB-style
alpar@1:  * compressed-column form, with sorted row indices in each column, and no
alpar@1:  * duplicate entries.  Diagonal entries may be present, but they are ignored.
alpar@1:  * Row indices of column j of A are stored in Ai [Ap [j] ... Ap [j+1]-1].
alpar@1:  * Ap [0] must be zero, and nz = Ap [n] is the number of entries in A.  The
alpar@1:  * size of the matrix, n, must be greater than or equal to zero.
alpar@1:  *
alpar@1:  * This routine must be preceded by a call to AMD_aat, which computes the
alpar@1:  * number of entries in each row/column in A+A', excluding the diagonal.
alpar@1:  * Len [j], on input, is the number of entries in row/column j of A+A'.  This
alpar@1:  * routine constructs the matrix A+A' and then calls AMD_2.  No error checking
alpar@1:  * is performed (this was done in AMD_valid).
alpar@1:  */
alpar@1: 
alpar@1: #include "amd_internal.h"
alpar@1: 
alpar@1: GLOBAL void AMD_1
alpar@1: (
alpar@1:     Int n,              /* n > 0 */
alpar@1:     const Int Ap [ ],   /* input of size n+1, not modified */
alpar@1:     const Int Ai [ ],   /* input of size nz = Ap [n], not modified */
alpar@1:     Int P [ ],          /* size n output permutation */
alpar@1:     Int Pinv [ ],       /* size n output inverse permutation */
alpar@1:     Int Len [ ],        /* size n input, undefined on output */
alpar@1:     Int slen,           /* slen >= sum (Len [0..n-1]) + 7n,
alpar@1:                          * ideally slen = 1.2 * sum (Len) + 8n */
alpar@1:     Int S [ ],          /* size slen workspace */
alpar@1:     double Control [ ], /* input array of size AMD_CONTROL */
alpar@1:     double Info [ ]     /* output array of size AMD_INFO */
alpar@1: )
alpar@1: {
alpar@1:     Int i, j, k, p, pfree, iwlen, pj, p1, p2, pj2, *Iw, *Pe, *Nv, *Head,
alpar@1:         *Elen, *Degree, *s, *W, *Sp, *Tp ;
alpar@1: 
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1:     /* construct the matrix for AMD_2 */
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1: 
alpar@1:     ASSERT (n > 0) ;
alpar@1: 
alpar@1:     iwlen = slen - 6*n ;
alpar@1:     s = S ;
alpar@1:     Pe = s ;        s += n ;
alpar@1:     Nv = s ;        s += n ;
alpar@1:     Head = s ;      s += n ;
alpar@1:     Elen = s ;      s += n ;
alpar@1:     Degree = s ;    s += n ;
alpar@1:     W = s ;         s += n ;
alpar@1:     Iw = s ;        s += iwlen ;
alpar@1: 
alpar@1:     ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ;
alpar@1: 
alpar@1:     /* construct the pointers for A+A' */
alpar@1:     Sp = Nv ;                   /* use Nv and W as workspace for Sp and Tp [ */
alpar@1:     Tp = W ;
alpar@1:     pfree = 0 ;
alpar@1:     for (j = 0 ; j < n ; j++)
alpar@1:     {
alpar@1:         Pe [j] = pfree ;
alpar@1:         Sp [j] = pfree ;
alpar@1:         pfree += Len [j] ;
alpar@1:     }
alpar@1: 
alpar@1:     /* Note that this restriction on iwlen is slightly more restrictive than
alpar@1:      * what is strictly required in AMD_2.  AMD_2 can operate with no elbow
alpar@1:      * room at all, but it will be very slow.  For better performance, at
alpar@1:      * least size-n elbow room is enforced. */
alpar@1:     ASSERT (iwlen >= pfree + n) ;
alpar@1: 
alpar@1: #ifndef NDEBUG
alpar@1:     for (p = 0 ; p < iwlen ; p++) Iw [p] = EMPTY ;
alpar@1: #endif
alpar@1: 
alpar@1:     for (k = 0 ; k < n ; k++)
alpar@1:     {
alpar@1:         AMD_DEBUG1 (("Construct row/column k= "ID" of A+A'\n", k))  ;
alpar@1:         p1 = Ap [k] ;
alpar@1:         p2 = Ap [k+1] ;
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:             ASSERT (j >= 0 && j < n) ;
alpar@1:             if (j < k)
alpar@1:             {
alpar@1:                 /* entry A (j,k) in the strictly upper triangular part */
alpar@1:                 ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
alpar@1:                 ASSERT (Sp [k] < (k == n-1 ? pfree : Pe [k+1])) ;
alpar@1:                 Iw [Sp [j]++] = k ;
alpar@1:                 Iw [Sp [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:                 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 (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:                 ASSERT (i >= 0 && i < n) ;
alpar@1:                 if (i < k)
alpar@1:                 {
alpar@1:                     /* A (i,j) is only in the lower part, not in upper */
alpar@1:                     ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;
alpar@1:                     ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
alpar@1:                     Iw [Sp [i]++] = j ;
alpar@1:                     Iw [Sp [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:                     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] = 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:             ASSERT (i >= 0 && i < n) ;
alpar@1:             /* A (i,j) is only in the lower part, not in upper */
alpar@1:             ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;
alpar@1:             ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
alpar@1:             Iw [Sp [i]++] = j ;
alpar@1:             Iw [Sp [j]++] = i ;
alpar@1:         }
alpar@1:     }
alpar@1: 
alpar@1: #ifndef NDEBUG
alpar@1:     for (j = 0 ; j < n-1 ; j++) ASSERT (Sp [j] == Pe [j+1]) ;
alpar@1:     ASSERT (Sp [n-1] == pfree) ;
alpar@1: #endif
alpar@1: 
alpar@1:     /* Tp and Sp no longer needed ] */
alpar@1: 
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1:     /* order the matrix */
alpar@1:     /* --------------------------------------------------------------------- */
alpar@1: 
alpar@1:     AMD_2 (n, Pe, Iw, Len, iwlen, pfree,
alpar@1:         Nv, Pinv, P, Head, Elen, Degree, W, Control, Info) ;
alpar@1: }