1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud (at) inria.fr> 5 6 /* 7 8 NOTE: this routine has been adapted from the CSparse library: 9 10 Copyright (c) 2006, Timothy A. Davis. 11 http://www.cise.ufl.edu/research/sparse/CSparse 12 13 CSparse is free software; you can redistribute it and/or 14 modify it under the terms of the GNU Lesser General Public 15 License as published by the Free Software Foundation; either 16 version 2.1 of the License, or (at your option) any later version. 17 18 CSparse is distributed in the hope that it will be useful, 19 but WITHOUT ANY WARRANTY; without even the implied warranty of 20 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 21 Lesser General Public License for more details. 22 23 You should have received a copy of the GNU Lesser General Public 24 License along with this Module; if not, write to the Free Software 25 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 26 27 */ 28 29 #include "../Core/util/NonMPL2.h" 30 31 #ifndef EIGEN_SPARSE_AMD_H 32 #define EIGEN_SPARSE_AMD_H 33 34 namespace Eigen { 35 36 namespace internal { 37 38 template<typename T> inline T amd_flip(const T& i) { return -i-2; } 39 template<typename T> inline T amd_unflip(const T& i) { return i<0 ? amd_flip(i) : i; } 40 template<typename T0, typename T1> inline bool amd_marked(const T0* w, const T1& j) { return w[j]<0; } 41 template<typename T0, typename T1> inline void amd_mark(const T0* w, const T1& j) { return w[j] = amd_flip(w[j]); } 42 43 /* clear w */ 44 template<typename Index> 45 static int cs_wclear (Index mark, Index lemax, Index *w, Index n) 46 { 47 Index k; 48 if(mark < 2 || (mark + lemax < 0)) 49 { 50 for(k = 0; k < n; k++) 51 if(w[k] != 0) 52 w[k] = 1; 53 mark = 2; 54 } 55 return (mark); /* at this point, w[0..n-1] < mark holds */ 56 } 57 58 /* depth-first search and postorder of a tree rooted at node j */ 59 template<typename Index> 60 Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack) 61 { 62 int i, p, top = 0; 63 if(!head || !next || !post || !stack) return (-1); /* check inputs */ 64 stack[0] = j; /* place j on the stack */ 65 while (top >= 0) /* while (stack is not empty) */ 66 { 67 p = stack[top]; /* p = top of stack */ 68 i = head[p]; /* i = youngest child of p */ 69 if(i == -1) 70 { 71 top--; /* p has no unordered children left */ 72 post[k++] = p; /* node p is the kth postordered node */ 73 } 74 else 75 { 76 head[p] = next[i]; /* remove i from children of p */ 77 stack[++top] = i; /* start dfs on child node i */ 78 } 79 } 80 return k; 81 } 82 83 84 /** \internal 85 * \ingroup OrderingMethods_Module 86 * Approximate minimum degree ordering algorithm. 87 * \returns the permutation P reducing the fill-in of the input matrix \a C 88 * The input matrix \a C must be a selfadjoint compressed column major SparseMatrix object. Both the upper and lower parts have to be stored, but the diagonal entries are optional. 89 * On exit the values of C are destroyed */ 90 template<typename Scalar, typename Index> 91 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic,Dynamic,Index>& perm) 92 { 93 using std::sqrt; 94 95 int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1, 96 k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi, 97 ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t; 98 unsigned int h; 99 100 Index n = C.cols(); 101 dense = std::max<Index> (16, Index(10 * sqrt(double(n)))); /* find dense threshold */ 102 dense = std::min<Index> (n-2, dense); 103 104 Index cnz = C.nonZeros(); 105 perm.resize(n+1); 106 t = cnz + cnz/5 + 2*n; /* add elbow room to C */ 107 C.resizeNonZeros(t); 108 109 Index* W = new Index[8*(n+1)]; /* get workspace */ 110 Index* len = W; 111 Index* nv = W + (n+1); 112 Index* next = W + 2*(n+1); 113 Index* head = W + 3*(n+1); 114 Index* elen = W + 4*(n+1); 115 Index* degree = W + 5*(n+1); 116 Index* w = W + 6*(n+1); 117 Index* hhead = W + 7*(n+1); 118 Index* last = perm.indices().data(); /* use P as workspace for last */ 119 120 /* --- Initialize quotient graph ---------------------------------------- */ 121 Index* Cp = C.outerIndexPtr(); 122 Index* Ci = C.innerIndexPtr(); 123 for(k = 0; k < n; k++) 124 len[k] = Cp[k+1] - Cp[k]; 125 len[n] = 0; 126 nzmax = t; 127 128 for(i = 0; i <= n; i++) 129 { 130 head[i] = -1; // degree list i is empty 131 last[i] = -1; 132 next[i] = -1; 133 hhead[i] = -1; // hash list i is empty 134 nv[i] = 1; // node i is just one node 135 w[i] = 1; // node i is alive 136 elen[i] = 0; // Ek of node i is empty 137 degree[i] = len[i]; // degree of node i 138 } 139 mark = internal::cs_wclear<Index>(0, 0, w, n); /* clear w */ 140 elen[n] = -2; /* n is a dead element */ 141 Cp[n] = -1; /* n is a root of assembly tree */ 142 w[n] = 0; /* n is a dead element */ 143 144 /* --- Initialize degree lists ------------------------------------------ */ 145 for(i = 0; i < n; i++) 146 { 147 d = degree[i]; 148 if(d == 0) /* node i is empty */ 149 { 150 elen[i] = -2; /* element i is dead */ 151 nel++; 152 Cp[i] = -1; /* i is a root of assembly tree */ 153 w[i] = 0; 154 } 155 else if(d > dense) /* node i is dense */ 156 { 157 nv[i] = 0; /* absorb i into element n */ 158 elen[i] = -1; /* node i is dead */ 159 nel++; 160 Cp[i] = amd_flip (n); 161 nv[n]++; 162 } 163 else 164 { 165 if(head[d] != -1) last[head[d]] = i; 166 next[i] = head[d]; /* put node i in degree list d */ 167 head[d] = i; 168 } 169 } 170 171 while (nel < n) /* while (selecting pivots) do */ 172 { 173 /* --- Select node of minimum approximate degree -------------------- */ 174 for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {} 175 if(next[k] != -1) last[next[k]] = -1; 176 head[mindeg] = next[k]; /* remove k from degree list */ 177 elenk = elen[k]; /* elenk = |Ek| */ 178 nvk = nv[k]; /* # of nodes k represents */ 179 nel += nvk; /* nv[k] nodes of A eliminated */ 180 181 /* --- Garbage collection ------------------------------------------- */ 182 if(elenk > 0 && cnz + mindeg >= nzmax) 183 { 184 for(j = 0; j < n; j++) 185 { 186 if((p = Cp[j]) >= 0) /* j is a live node or element */ 187 { 188 Cp[j] = Ci[p]; /* save first entry of object */ 189 Ci[p] = amd_flip (j); /* first entry is now amd_flip(j) */ 190 } 191 } 192 for(q = 0, p = 0; p < cnz; ) /* scan all of memory */ 193 { 194 if((j = amd_flip (Ci[p++])) >= 0) /* found object j */ 195 { 196 Ci[q] = Cp[j]; /* restore first entry of object */ 197 Cp[j] = q++; /* new pointer to object j */ 198 for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++]; 199 } 200 } 201 cnz = q; /* Ci[cnz...nzmax-1] now free */ 202 } 203 204 /* --- Construct new element ---------------------------------------- */ 205 dk = 0; 206 nv[k] = -nvk; /* flag k as in Lk */ 207 p = Cp[k]; 208 pk1 = (elenk == 0) ? p : cnz; /* do in place if elen[k] == 0 */ 209 pk2 = pk1; 210 for(k1 = 1; k1 <= elenk + 1; k1++) 211 { 212 if(k1 > elenk) 213 { 214 e = k; /* search the nodes in k */ 215 pj = p; /* list of nodes starts at Ci[pj]*/ 216 ln = len[k] - elenk; /* length of list of nodes in k */ 217 } 218 else 219 { 220 e = Ci[p++]; /* search the nodes in e */ 221 pj = Cp[e]; 222 ln = len[e]; /* length of list of nodes in e */ 223 } 224 for(k2 = 1; k2 <= ln; k2++) 225 { 226 i = Ci[pj++]; 227 if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */ 228 dk += nvi; /* degree[Lk] += size of node i */ 229 nv[i] = -nvi; /* negate nv[i] to denote i in Lk*/ 230 Ci[pk2++] = i; /* place i in Lk */ 231 if(next[i] != -1) last[next[i]] = last[i]; 232 if(last[i] != -1) /* remove i from degree list */ 233 { 234 next[last[i]] = next[i]; 235 } 236 else 237 { 238 head[degree[i]] = next[i]; 239 } 240 } 241 if(e != k) 242 { 243 Cp[e] = amd_flip (k); /* absorb e into k */ 244 w[e] = 0; /* e is now a dead element */ 245 } 246 } 247 if(elenk != 0) cnz = pk2; /* Ci[cnz...nzmax] is free */ 248 degree[k] = dk; /* external degree of k - |Lk\i| */ 249 Cp[k] = pk1; /* element k is in Ci[pk1..pk2-1] */ 250 len[k] = pk2 - pk1; 251 elen[k] = -2; /* k is now an element */ 252 253 /* --- Find set differences ----------------------------------------- */ 254 mark = internal::cs_wclear<Index>(mark, lemax, w, n); /* clear w if necessary */ 255 for(pk = pk1; pk < pk2; pk++) /* scan 1: find |Le\Lk| */ 256 { 257 i = Ci[pk]; 258 if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */ 259 nvi = -nv[i]; /* nv[i] was negated */ 260 wnvi = mark - nvi; 261 for(p = Cp[i]; p <= Cp[i] + eln - 1; p++) /* scan Ei */ 262 { 263 e = Ci[p]; 264 if(w[e] >= mark) 265 { 266 w[e] -= nvi; /* decrement |Le\Lk| */ 267 } 268 else if(w[e] != 0) /* ensure e is a live element */ 269 { 270 w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */ 271 } 272 } 273 } 274 275 /* --- Degree update ------------------------------------------------ */ 276 for(pk = pk1; pk < pk2; pk++) /* scan2: degree update */ 277 { 278 i = Ci[pk]; /* consider node i in Lk */ 279 p1 = Cp[i]; 280 p2 = p1 + elen[i] - 1; 281 pn = p1; 282 for(h = 0, d = 0, p = p1; p <= p2; p++) /* scan Ei */ 283 { 284 e = Ci[p]; 285 if(w[e] != 0) /* e is an unabsorbed element */ 286 { 287 dext = w[e] - mark; /* dext = |Le\Lk| */ 288 if(dext > 0) 289 { 290 d += dext; /* sum up the set differences */ 291 Ci[pn++] = e; /* keep e in Ei */ 292 h += e; /* compute the hash of node i */ 293 } 294 else 295 { 296 Cp[e] = amd_flip (k); /* aggressive absorb. e->k */ 297 w[e] = 0; /* e is a dead element */ 298 } 299 } 300 } 301 elen[i] = pn - p1 + 1; /* elen[i] = |Ei| */ 302 p3 = pn; 303 p4 = p1 + len[i]; 304 for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */ 305 { 306 j = Ci[p]; 307 if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */ 308 d += nvj; /* degree(i) += |j| */ 309 Ci[pn++] = j; /* place j in node list of i */ 310 h += j; /* compute hash for node i */ 311 } 312 if(d == 0) /* check for mass elimination */ 313 { 314 Cp[i] = amd_flip (k); /* absorb i into k */ 315 nvi = -nv[i]; 316 dk -= nvi; /* |Lk| -= |i| */ 317 nvk += nvi; /* |k| += nv[i] */ 318 nel += nvi; 319 nv[i] = 0; 320 elen[i] = -1; /* node i is dead */ 321 } 322 else 323 { 324 degree[i] = std::min<Index> (degree[i], d); /* update degree(i) */ 325 Ci[pn] = Ci[p3]; /* move first node to end */ 326 Ci[p3] = Ci[p1]; /* move 1st el. to end of Ei */ 327 Ci[p1] = k; /* add k as 1st element in of Ei */ 328 len[i] = pn - p1 + 1; /* new len of adj. list of node i */ 329 h %= n; /* finalize hash of i */ 330 next[i] = hhead[h]; /* place i in hash bucket */ 331 hhead[h] = i; 332 last[i] = h; /* save hash of i in last[i] */ 333 } 334 } /* scan2 is done */ 335 degree[k] = dk; /* finalize |Lk| */ 336 lemax = std::max<Index>(lemax, dk); 337 mark = internal::cs_wclear<Index>(mark+lemax, lemax, w, n); /* clear w */ 338 339 /* --- Supernode detection ------------------------------------------ */ 340 for(pk = pk1; pk < pk2; pk++) 341 { 342 i = Ci[pk]; 343 if(nv[i] >= 0) continue; /* skip if i is dead */ 344 h = last[i]; /* scan hash bucket of node i */ 345 i = hhead[h]; 346 hhead[h] = -1; /* hash bucket will be empty */ 347 for(; i != -1 && next[i] != -1; i = next[i], mark++) 348 { 349 ln = len[i]; 350 eln = elen[i]; 351 for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark; 352 jlast = i; 353 for(j = next[i]; j != -1; ) /* compare i with all j */ 354 { 355 ok = (len[j] == ln) && (elen[j] == eln); 356 for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++) 357 { 358 if(w[Ci[p]] != mark) ok = 0; /* compare i and j*/ 359 } 360 if(ok) /* i and j are identical */ 361 { 362 Cp[j] = amd_flip (i); /* absorb j into i */ 363 nv[i] += nv[j]; 364 nv[j] = 0; 365 elen[j] = -1; /* node j is dead */ 366 j = next[j]; /* delete j from hash bucket */ 367 next[jlast] = j; 368 } 369 else 370 { 371 jlast = j; /* j and i are different */ 372 j = next[j]; 373 } 374 } 375 } 376 } 377 378 /* --- Finalize new element------------------------------------------ */ 379 for(p = pk1, pk = pk1; pk < pk2; pk++) /* finalize Lk */ 380 { 381 i = Ci[pk]; 382 if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */ 383 nv[i] = nvi; /* restore nv[i] */ 384 d = degree[i] + dk - nvi; /* compute external degree(i) */ 385 d = std::min<Index> (d, n - nel - nvi); 386 if(head[d] != -1) last[head[d]] = i; 387 next[i] = head[d]; /* put i back in degree list */ 388 last[i] = -1; 389 head[d] = i; 390 mindeg = std::min<Index> (mindeg, d); /* find new minimum degree */ 391 degree[i] = d; 392 Ci[p++] = i; /* place i in Lk */ 393 } 394 nv[k] = nvk; /* # nodes absorbed into k */ 395 if((len[k] = p-pk1) == 0) /* length of adj list of element k*/ 396 { 397 Cp[k] = -1; /* k is a root of the tree */ 398 w[k] = 0; /* k is now a dead element */ 399 } 400 if(elenk != 0) cnz = p; /* free unused space in Lk */ 401 } 402 403 /* --- Postordering ----------------------------------------------------- */ 404 for(i = 0; i < n; i++) Cp[i] = amd_flip (Cp[i]);/* fix assembly tree */ 405 for(j = 0; j <= n; j++) head[j] = -1; 406 for(j = n; j >= 0; j--) /* place unordered nodes in lists */ 407 { 408 if(nv[j] > 0) continue; /* skip if j is an element */ 409 next[j] = head[Cp[j]]; /* place j in list of its parent */ 410 head[Cp[j]] = j; 411 } 412 for(e = n; e >= 0; e--) /* place elements in lists */ 413 { 414 if(nv[e] <= 0) continue; /* skip unless e is an element */ 415 if(Cp[e] != -1) 416 { 417 next[e] = head[Cp[e]]; /* place e in list of its parent */ 418 head[Cp[e]] = e; 419 } 420 } 421 for(k = 0, i = 0; i <= n; i++) /* postorder the assembly tree */ 422 { 423 if(Cp[i] == -1) k = internal::cs_tdfs<Index>(i, k, head, next, perm.indices().data(), w); 424 } 425 426 perm.indices().conservativeResize(n); 427 428 delete[] W; 429 } 430 431 } // namespace internal 432 433 } // end namespace Eigen 434 435 #endif // EIGEN_SPARSE_AMD_H 436
