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