Amesos Package Browser (Single Doxygen Collection) Development
amesos_camd_1.c
Go to the documentation of this file.
```00001 /* ========================================================================= */
00002 /* === CAMD_1 ============================================================== */
00003 /* ========================================================================= */
00004
00005 /* ------------------------------------------------------------------------- */
00006 /* CAMD, Copyright (c) Timothy A. Davis, Yanqing Chen,           */
00007 /* Patrick R. Amestoy, and Iain S. Duff.  See ../README.txt for License.     */
00008 /* email: davis at cise.ufl.edu    CISE Department, Univ. of Florida.        */
00009 /* web: http://www.cise.ufl.edu/research/sparse/camd                         */
00010 /* ------------------------------------------------------------------------- */
00011
00012 /* CAMD_1: Construct A+A' for a sparse matrix A and perform the CAMD ordering.
00013  *
00014  * The n-by-n sparse matrix A can be unsymmetric.  It is stored in MATLAB-style
00015  * compressed-column form, with sorted row indices in each column, and no
00016  * duplicate entries.  Diagonal entries may be present, but they are ignored.
00017  * Row indices of column j of A are stored in Ai [Ap [j] ... Ap [j+1]-1].
00018  * Ap [0] must be zero, and nz = Ap [n] is the number of entries in A.  The
00019  * size of the matrix, n, must be greater than or equal to zero.
00020  *
00021  * This routine must be preceded by a call to CAMD_aat, which computes the
00022  * number of entries in each row/column in A+A', excluding the diagonal.
00023  * Len [j], on input, is the number of entries in row/column j of A+A'.  This
00024  * routine constructs the matrix A+A' and then calls CAMD_2.  No error checking
00025  * is performed (this was done in CAMD_valid).
00026  */
00027
00028 #include "amesos_camd_internal.h"
00029
00030 GLOBAL void CAMD_1
00031 (
00032     Int n,    /* n > 0 */
00033     const Int Ap [ ], /* input of size n+1, not modified */
00034     const Int Ai [ ], /* input of size nz = Ap [n], not modified */
00035     Int P [ ],    /* size n output permutation */
00036     Int Pinv [ ], /* size n output inverse permutation */
00037     Int Len [ ],  /* size n input, undefined on output */
00038     Int slen,   /* slen >= sum (Len [0..n-1]) + 7n+2,
00039        * ideally slen = 1.2 * sum (Len) + 8n+2 */
00040     Int S [ ],    /* size slen workspace */
00041     double Control [ ], /* input array of size CAMD_CONTROL */
00042     double Info [ ],  /* output array of size CAMD_INFO */
00043     const Int C [ ] /* Constraint set of size n */
00044 )
00045 {
00046     Int i, j, k, p, pfree, iwlen, pj, p1, p2, pj2, *Iw, *Pe, *Nv, *Head,
00047   *Elen, *Degree, *s, *W, *Sp, *Tp, *BucketSet ;
00048
00049     /* --------------------------------------------------------------------- */
00050     /* construct the matrix for CAMD_2 */
00051     /* --------------------------------------------------------------------- */
00052
00053     ASSERT (n > 0) ;
00054
00055     iwlen = slen - (7*n+2) ;  /* allocate 7*n+2 workspace from S */
00056     s = S ;
00057     Pe = s ;      s += n ;
00058     Nv = s ;      s += n ;
00059     Head = s ;      s += n+1 ;  /* NOTE: was size n in AMD; size n+1 in CAMD */
00060     Elen = s ;      s += n ;
00061     Degree = s ;    s += n ;
00062     W = s ;     s += n+1 ;  /* NOTE: was size n in AMD; size n+1 in CAMD */
00063     BucketSet = s ; s += n ;
00064     Iw = s ;      s += iwlen ;
00065
00066     ASSERT (CAMD_valid (n, n, Ap, Ai) == CAMD_OK) ;
00067     ASSERT (CAMD_cvalid (n, C)) ;
00068
00069     /* construct the pointers for A+A' */
00070     Sp = Nv ;     /* use Nv and W as workspace for Sp and Tp [ */
00071     Tp = W ;
00072     pfree = 0 ;
00073     for (j = 0 ; j < n ; j++)
00074     {
00075   Pe [j] = pfree ;
00076   Sp [j] = pfree ;
00077   pfree += Len [j] ;
00078     }
00079
00080     /* Note that this restriction on iwlen is slightly more restrictive than
00081      * what is strictly required in CAMD_2.  CAMD_2 can operate with no elbow
00082      * room at all, but it will be very slow.  For better performance, at
00083      * least size-n elbow room is enforced. */
00084     ASSERT (iwlen >= pfree + n) ;
00085
00086 #ifndef NDEBUG
00087     for (p = 0 ; p < iwlen ; p++) Iw [p] = EMPTY ;
00088 #endif
00089
00090     for (k = 0 ; k < n ; k++)
00091     {
00092   CAMD_DEBUG1 (("Construct row/column k= "ID" of A+A'\n", k))  ;
00093   p1 = Ap [k] ;
00094   p2 = Ap [k+1] ;
00095
00096   /* construct A+A' */
00097   for (p = p1 ; p < p2 ; )
00098   {
00099       /* scan the upper triangular part of A */
00100       j = Ai [p] ;
00101       ASSERT (j >= 0 && j < n) ;
00102       if (j < k)
00103       {
00104     /* entry A (j,k) in the strictly upper triangular part */
00105     ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
00106     ASSERT (Sp [k] < (k == n-1 ? pfree : Pe [k+1])) ;
00107     Iw [Sp [j]++] = k ;
00108     Iw [Sp [k]++] = j ;
00109     p++ ;
00110       }
00111       else if (j == k)
00112       {
00113     /* skip the diagonal */
00114     p++ ;
00115     break ;
00116       }
00117       else /* j > k */
00118       {
00119     /* first entry below the diagonal */
00120     break ;
00121       }
00122       /* scan lower triangular part of A, in column j until reaching
00123        * row k.  Start where last scan left off. */
00124       ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ;
00125       pj2 = Ap [j+1] ;
00126       for (pj = Tp [j] ; pj < pj2 ; )
00127       {
00128     i = Ai [pj] ;
00129     ASSERT (i >= 0 && i < n) ;
00130     if (i < k)
00131     {
00132         /* A (i,j) is only in the lower part, not in upper */
00133         ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;
00134         ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
00135         Iw [Sp [i]++] = j ;
00136         Iw [Sp [j]++] = i ;
00137         pj++ ;
00138     }
00139     else if (i == k)
00140     {
00141         /* entry A (k,j) in lower part and A (j,k) in upper */
00142         pj++ ;
00143         break ;
00144     }
00145     else /* i > k */
00146     {
00147         /* consider this entry later, when k advances to i */
00148         break ;
00149     }
00150       }
00151       Tp [j] = pj ;
00152   }
00153   Tp [k] = p ;
00154     }
00155
00156     /* clean up, for remaining mismatched entries */
00157     for (j = 0 ; j < n ; j++)
00158     {
00159   for (pj = Tp [j] ; pj < Ap [j+1] ; pj++)
00160   {
00161       i = Ai [pj] ;
00162       ASSERT (i >= 0 && i < n) ;
00163       /* A (i,j) is only in the lower part, not in upper */
00164       ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ;
00165       ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ;
00166       Iw [Sp [i]++] = j ;
00167       Iw [Sp [j]++] = i ;
00168   }
00169     }
00170
00171 #ifndef NDEBUG
00172     for (j = 0 ; j < n-1 ; j++) ASSERT (Sp [j] == Pe [j+1]) ;
00173     ASSERT (Sp [n-1] == pfree) ;
00174 #endif
00175
00176     /* Tp and Sp no longer needed ] */
00177
00178     /* --------------------------------------------------------------------- */
00179     /* order the matrix */
00180     /* --------------------------------------------------------------------- */
00181
00182     CAMD_2 (n, Pe, Iw, Len, iwlen, pfree,
00183   Nv, Pinv, P, Head, Elen, Degree, W, Control, Info, C, BucketSet) ;
00184 }
```