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 bool has_diag = false; 148 for(p = Cp[i]; p<Cp[i+1]; ++p) 149 if(Ci[p]==i) 150 { 151 has_diag = true; 152 break; 153 } 154 155 d = degree[i]; 156 if(d == 1) /* node i is empty */ 157 { 158 elen[i] = -2; /* element i is dead */ 159 nel++; 160 Cp[i] = -1; /* i is a root of assembly tree */ 161 w[i] = 0; 162 } 163 else if(d > dense || !has_diag) /* node i is dense or has no structural diagonal element */ 164 { 165 nv[i] = 0; /* absorb i into element n */ 166 elen[i] = -1; /* node i is dead */ 167 nel++; 168 Cp[i] = amd_flip (n); 169 nv[n]++; 170 } 171 else 172 { 173 if(head[d] != -1) last[head[d]] = i; 174 next[i] = head[d]; /* put node i in degree list d */ 175 head[d] = i; 176 } 177 } 178 179 elen[n] = -2; /* n is a dead element */ 180 Cp[n] = -1; /* n is a root of assembly tree */ 181 w[n] = 0; /* n is a dead element */ 182 183 while (nel < n) /* while (selecting pivots) do */ 184 { 185 /* --- Select node of minimum approximate degree -------------------- */ 186 for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {} 187 if(next[k] != -1) last[next[k]] = -1; 188 head[mindeg] = next[k]; /* remove k from degree list */ 189 elenk = elen[k]; /* elenk = |Ek| */ 190 nvk = nv[k]; /* # of nodes k represents */ 191 nel += nvk; /* nv[k] nodes of A eliminated */ 192 193 /* --- Garbage collection ------------------------------------------- */ 194 if(elenk > 0 && cnz + mindeg >= nzmax) 195 { 196 for(j = 0; j < n; j++) 197 { 198 if((p = Cp[j]) >= 0) /* j is a live node or element */ 199 { 200 Cp[j] = Ci[p]; /* save first entry of object */ 201 Ci[p] = amd_flip (j); /* first entry is now amd_flip(j) */ 202 } 203 } 204 for(q = 0, p = 0; p < cnz; ) /* scan all of memory */ 205 { 206 if((j = amd_flip (Ci[p++])) >= 0) /* found object j */ 207 { 208 Ci[q] = Cp[j]; /* restore first entry of object */ 209 Cp[j] = q++; /* new pointer to object j */ 210 for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++]; 211 } 212 } 213 cnz = q; /* Ci[cnz...nzmax-1] now free */ 214 } 215 216 /* --- Construct new element ---------------------------------------- */ 217 dk = 0; 218 nv[k] = -nvk; /* flag k as in Lk */ 219 p = Cp[k]; 220 pk1 = (elenk == 0) ? p : cnz; /* do in place if elen[k] == 0 */ 221 pk2 = pk1; 222 for(k1 = 1; k1 <= elenk + 1; k1++) 223 { 224 if(k1 > elenk) 225 { 226 e = k; /* search the nodes in k */ 227 pj = p; /* list of nodes starts at Ci[pj]*/ 228 ln = len[k] - elenk; /* length of list of nodes in k */ 229 } 230 else 231 { 232 e = Ci[p++]; /* search the nodes in e */ 233 pj = Cp[e]; 234 ln = len[e]; /* length of list of nodes in e */ 235 } 236 for(k2 = 1; k2 <= ln; k2++) 237 { 238 i = Ci[pj++]; 239 if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */ 240 dk += nvi; /* degree[Lk] += size of node i */ 241 nv[i] = -nvi; /* negate nv[i] to denote i in Lk*/ 242 Ci[pk2++] = i; /* place i in Lk */ 243 if(next[i] != -1) last[next[i]] = last[i]; 244 if(last[i] != -1) /* remove i from degree list */ 245 { 246 next[last[i]] = next[i]; 247 } 248 else 249 { 250 head[degree[i]] = next[i]; 251 } 252 } 253 if(e != k) 254 { 255 Cp[e] = amd_flip (k); /* absorb e into k */ 256 w[e] = 0; /* e is now a dead element */ 257 } 258 } 259 if(elenk != 0) cnz = pk2; /* Ci[cnz...nzmax] is free */ 260 degree[k] = dk; /* external degree of k - |Lk\i| */ 261 Cp[k] = pk1; /* element k is in Ci[pk1..pk2-1] */ 262 len[k] = pk2 - pk1; 263 elen[k] = -2; /* k is now an element */ 264 265 /* --- Find set differences ----------------------------------------- */ 266 mark = internal::cs_wclear<Index>(mark, lemax, w, n); /* clear w if necessary */ 267 for(pk = pk1; pk < pk2; pk++) /* scan 1: find |Le\Lk| */ 268 { 269 i = Ci[pk]; 270 if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */ 271 nvi = -nv[i]; /* nv[i] was negated */ 272 wnvi = mark - nvi; 273 for(p = Cp[i]; p <= Cp[i] + eln - 1; p++) /* scan Ei */ 274 { 275 e = Ci[p]; 276 if(w[e] >= mark) 277 { 278 w[e] -= nvi; /* decrement |Le\Lk| */ 279 } 280 else if(w[e] != 0) /* ensure e is a live element */ 281 { 282 w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */ 283 } 284 } 285 } 286 287 /* --- Degree update ------------------------------------------------ */ 288 for(pk = pk1; pk < pk2; pk++) /* scan2: degree update */ 289 { 290 i = Ci[pk]; /* consider node i in Lk */ 291 p1 = Cp[i]; 292 p2 = p1 + elen[i] - 1; 293 pn = p1; 294 for(h = 0, d = 0, p = p1; p <= p2; p++) /* scan Ei */ 295 { 296 e = Ci[p]; 297 if(w[e] != 0) /* e is an unabsorbed element */ 298 { 299 dext = w[e] - mark; /* dext = |Le\Lk| */ 300 if(dext > 0) 301 { 302 d += dext; /* sum up the set differences */ 303 Ci[pn++] = e; /* keep e in Ei */ 304 h += e; /* compute the hash of node i */ 305 } 306 else 307 { 308 Cp[e] = amd_flip (k); /* aggressive absorb. e->k */ 309 w[e] = 0; /* e is a dead element */ 310 } 311 } 312 } 313 elen[i] = pn - p1 + 1; /* elen[i] = |Ei| */ 314 p3 = pn; 315 p4 = p1 + len[i]; 316 for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */ 317 { 318 j = Ci[p]; 319 if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */ 320 d += nvj; /* degree(i) += |j| */ 321 Ci[pn++] = j; /* place j in node list of i */ 322 h += j; /* compute hash for node i */ 323 } 324 if(d == 0) /* check for mass elimination */ 325 { 326 Cp[i] = amd_flip (k); /* absorb i into k */ 327 nvi = -nv[i]; 328 dk -= nvi; /* |Lk| -= |i| */ 329 nvk += nvi; /* |k| += nv[i] */ 330 nel += nvi; 331 nv[i] = 0; 332 elen[i] = -1; /* node i is dead */ 333 } 334 else 335 { 336 degree[i] = std::min<Index> (degree[i], d); /* update degree(i) */ 337 Ci[pn] = Ci[p3]; /* move first node to end */ 338 Ci[p3] = Ci[p1]; /* move 1st el. to end of Ei */ 339 Ci[p1] = k; /* add k as 1st element in of Ei */ 340 len[i] = pn - p1 + 1; /* new len of adj. list of node i */ 341 h %= n; /* finalize hash of i */ 342 next[i] = hhead[h]; /* place i in hash bucket */ 343 hhead[h] = i; 344 last[i] = h; /* save hash of i in last[i] */ 345 } 346 } /* scan2 is done */ 347 degree[k] = dk; /* finalize |Lk| */ 348 lemax = std::max<Index>(lemax, dk); 349 mark = internal::cs_wclear<Index>(mark+lemax, lemax, w, n); /* clear w */ 350 351 /* --- Supernode detection ------------------------------------------ */ 352 for(pk = pk1; pk < pk2; pk++) 353 { 354 i = Ci[pk]; 355 if(nv[i] >= 0) continue; /* skip if i is dead */ 356 h = last[i]; /* scan hash bucket of node i */ 357 i = hhead[h]; 358 hhead[h] = -1; /* hash bucket will be empty */ 359 for(; i != -1 && next[i] != -1; i = next[i], mark++) 360 { 361 ln = len[i]; 362 eln = elen[i]; 363 for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark; 364 jlast = i; 365 for(j = next[i]; j != -1; ) /* compare i with all j */ 366 { 367 ok = (len[j] == ln) && (elen[j] == eln); 368 for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++) 369 { 370 if(w[Ci[p]] != mark) ok = 0; /* compare i and j*/ 371 } 372 if(ok) /* i and j are identical */ 373 { 374 Cp[j] = amd_flip (i); /* absorb j into i */ 375 nv[i] += nv[j]; 376 nv[j] = 0; 377 elen[j] = -1; /* node j is dead */ 378 j = next[j]; /* delete j from hash bucket */ 379 next[jlast] = j; 380 } 381 else 382 { 383 jlast = j; /* j and i are different */ 384 j = next[j]; 385 } 386 } 387 } 388 } 389 390 /* --- Finalize new element------------------------------------------ */ 391 for(p = pk1, pk = pk1; pk < pk2; pk++) /* finalize Lk */ 392 { 393 i = Ci[pk]; 394 if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */ 395 nv[i] = nvi; /* restore nv[i] */ 396 d = degree[i] + dk - nvi; /* compute external degree(i) */ 397 d = std::min<Index> (d, n - nel - nvi); 398 if(head[d] != -1) last[head[d]] = i; 399 next[i] = head[d]; /* put i back in degree list */ 400 last[i] = -1; 401 head[d] = i; 402 mindeg = std::min<Index> (mindeg, d); /* find new minimum degree */ 403 degree[i] = d; 404 Ci[p++] = i; /* place i in Lk */ 405 } 406 nv[k] = nvk; /* # nodes absorbed into k */ 407 if((len[k] = p-pk1) == 0) /* length of adj list of element k*/ 408 { 409 Cp[k] = -1; /* k is a root of the tree */ 410 w[k] = 0; /* k is now a dead element */ 411 } 412 if(elenk != 0) cnz = p; /* free unused space in Lk */ 413 } 414 415 /* --- Postordering ----------------------------------------------------- */ 416 for(i = 0; i < n; i++) Cp[i] = amd_flip (Cp[i]);/* fix assembly tree */ 417 for(j = 0; j <= n; j++) head[j] = -1; 418 for(j = n; j >= 0; j--) /* place unordered nodes in lists */ 419 { 420 if(nv[j] > 0) continue; /* skip if j is an element */ 421 next[j] = head[Cp[j]]; /* place j in list of its parent */ 422 head[Cp[j]] = j; 423 } 424 for(e = n; e >= 0; e--) /* place elements in lists */ 425 { 426 if(nv[e] <= 0) continue; /* skip unless e is an element */ 427 if(Cp[e] != -1) 428 { 429 next[e] = head[Cp[e]]; /* place e in list of its parent */ 430 head[Cp[e]] = e; 431 } 432 } 433 for(k = 0, i = 0; i <= n; i++) /* postorder the assembly tree */ 434 { 435 if(Cp[i] == -1) k = internal::cs_tdfs<Index>(i, k, head, next, perm.indices().data(), w); 436 } 437 438 perm.indices().conservativeResize(n); 439 440 delete[] W; 441 } 442 443 } // namespace internal 444 445 } // end namespace Eigen 446 447 #endif // EIGEN_SPARSE_AMD_H 448
