1 /* ========================================================================= */ 2 /* === AMD_1 =============================================================== */ 3 /* ========================================================================= */ 4 5 /* ------------------------------------------------------------------------- */ 6 /* AMD, Copyright (c) Timothy A. Davis, */ 7 /* Patrick R. Amestoy, and Iain S. Duff. See ../README.txt for License. */ 8 /* email: davis at cise.ufl.edu CISE Department, Univ. of Florida. */ 9 /* web: http://www.cise.ufl.edu/research/sparse/amd */ 10 /* ------------------------------------------------------------------------- */ 11 12 /* AMD_1: Construct A+A' for a sparse matrix A and perform the AMD ordering. 13 * 14 * The n-by-n sparse matrix A can be unsymmetric. It is stored in MATLAB-style 15 * compressed-column form, with sorted row indices in each column, and no 16 * duplicate entries. Diagonal entries may be present, but they are ignored. 17 * Row indices of column j of A are stored in Ai [Ap [j] ... Ap [j+1]-1]. 18 * Ap [0] must be zero, and nz = Ap [n] is the number of entries in A. The 19 * size of the matrix, n, must be greater than or equal to zero. 20 * 21 * This routine must be preceded by a call to AMD_aat, which computes the 22 * number of entries in each row/column in A+A', excluding the diagonal. 23 * Len [j], on input, is the number of entries in row/column j of A+A'. This 24 * routine constructs the matrix A+A' and then calls AMD_2. No error checking 25 * is performed (this was done in AMD_valid). 26 */ 27 28 #include "amd_internal.h" 29 30 GLOBAL void AMD_1 31 ( 32 Int n, /* n > 0 */ 33 const Int Ap [ ], /* input of size n+1, not modified */ 34 const Int Ai [ ], /* input of size nz = Ap [n], not modified */ 35 Int P [ ], /* size n output permutation */ 36 Int Pinv [ ], /* size n output inverse permutation */ 37 Int Len [ ], /* size n input, undefined on output */ 38 Int slen, /* slen >= sum (Len [0..n-1]) + 7n, 39 * ideally slen = 1.2 * sum (Len) + 8n */ 40 Int S [ ], /* size slen workspace */ 41 double Control [ ], /* input array of size AMD_CONTROL */ 42 double Info [ ] /* output array of size AMD_INFO */ 43 ) 44 { 45 Int i, j, k, p, pfree, iwlen, pj, p1, p2, pj2, *Iw, *Pe, *Nv, *Head, 46 *Elen, *Degree, *s, *W, *Sp, *Tp ; 47 48 /* --------------------------------------------------------------------- */ 49 /* construct the matrix for AMD_2 */ 50 /* --------------------------------------------------------------------- */ 51 52 ASSERT (n > 0) ; 53 54 iwlen = slen - 6*n ; 55 s = S ; 56 Pe = s ; s += n ; 57 Nv = s ; s += n ; 58 Head = s ; s += n ; 59 Elen = s ; s += n ; 60 Degree = s ; s += n ; 61 W = s ; s += n ; 62 Iw = s ; s += iwlen ; 63 64 ASSERT (AMD_valid (n, n, Ap, Ai) == AMD_OK) ; 65 66 /* construct the pointers for A+A' */ 67 Sp = Nv ; /* use Nv and W as workspace for Sp and Tp [ */ 68 Tp = W ; 69 pfree = 0 ; 70 for (j = 0 ; j < n ; j++) 71 { 72 Pe [j] = pfree ; 73 Sp [j] = pfree ; 74 pfree += Len [j] ; 75 } 76 77 /* Note that this restriction on iwlen is slightly more restrictive than 78 * what is strictly required in AMD_2. AMD_2 can operate with no elbow 79 * room at all, but it will be very slow. For better performance, at 80 * least size-n elbow room is enforced. */ 81 ASSERT (iwlen >= pfree + n) ; 82 83 #ifndef NDEBUG 84 for (p = 0 ; p < iwlen ; p++) Iw [p] = EMPTY ; 85 #endif 86 87 for (k = 0 ; k < n ; k++) 88 { 89 AMD_DEBUG1 (("Construct row/column k= "ID" of A+A'\n", k)) ; 90 p1 = Ap [k] ; 91 p2 = Ap [k+1] ; 92 93 /* construct A+A' */ 94 for (p = p1 ; p < p2 ; ) 95 { 96 /* scan the upper triangular part of A */ 97 j = Ai [p] ; 98 ASSERT (j >= 0 && j < n) ; 99 if (j < k) 100 { 101 /* entry A (j,k) in the strictly upper triangular part */ 102 ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ; 103 ASSERT (Sp [k] < (k == n-1 ? pfree : Pe [k+1])) ; 104 Iw [Sp [j]++] = k ; 105 Iw [Sp [k]++] = j ; 106 p++ ; 107 } 108 else if (j == k) 109 { 110 /* skip the diagonal */ 111 p++ ; 112 break ; 113 } 114 else /* j > k */ 115 { 116 /* first entry below the diagonal */ 117 break ; 118 } 119 /* scan lower triangular part of A, in column j until reaching 120 * row k. Start where last scan left off. */ 121 ASSERT (Ap [j] <= Tp [j] && Tp [j] <= Ap [j+1]) ; 122 pj2 = Ap [j+1] ; 123 for (pj = Tp [j] ; pj < pj2 ; ) 124 { 125 i = Ai [pj] ; 126 ASSERT (i >= 0 && i < n) ; 127 if (i < k) 128 { 129 /* A (i,j) is only in the lower part, not in upper */ 130 ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ; 131 ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ; 132 Iw [Sp [i]++] = j ; 133 Iw [Sp [j]++] = i ; 134 pj++ ; 135 } 136 else if (i == k) 137 { 138 /* entry A (k,j) in lower part and A (j,k) in upper */ 139 pj++ ; 140 break ; 141 } 142 else /* i > k */ 143 { 144 /* consider this entry later, when k advances to i */ 145 break ; 146 } 147 } 148 Tp [j] = pj ; 149 } 150 Tp [k] = p ; 151 } 152 153 /* clean up, for remaining mismatched entries */ 154 for (j = 0 ; j < n ; j++) 155 { 156 for (pj = Tp [j] ; pj < Ap [j+1] ; pj++) 157 { 158 i = Ai [pj] ; 159 ASSERT (i >= 0 && i < n) ; 160 /* A (i,j) is only in the lower part, not in upper */ 161 ASSERT (Sp [i] < (i == n-1 ? pfree : Pe [i+1])) ; 162 ASSERT (Sp [j] < (j == n-1 ? pfree : Pe [j+1])) ; 163 Iw [Sp [i]++] = j ; 164 Iw [Sp [j]++] = i ; 165 } 166 } 167 168 #ifndef NDEBUG 169 for (j = 0 ; j < n-1 ; j++) ASSERT (Sp [j] == Pe [j+1]) ; 170 ASSERT (Sp [n-1] == pfree) ; 171 #endif 172 173 /* Tp and Sp no longer needed ] */ 174 175 /* --------------------------------------------------------------------- */ 176 /* order the matrix */ 177 /* --------------------------------------------------------------------- */ 178 179 AMD_2 (n, Pe, Iw, Len, iwlen, pfree, 180 Nv, Pinv, P, Head, Elen, Degree, W, Control, Info) ; 181 } 182