Amd.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
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@inria.fr> 5 6/* 7 8NOTE: this routine has been adapted from the CSparse library: 9 10Copyright (c) 2006, Timothy A. Davis. 11http://www.cise.ufl.edu/research/sparse/CSparse 12 13CSparse is free software; you can redistribute it and/or 14modify it under the terms of the GNU Lesser General Public 15License as published by the Free Software Foundation; either 16version 2.1 of the License, or (at your option) any later version. 17 18CSparse is distributed in the hope that it will be useful, 19but WITHOUT ANY WARRANTY; without even the implied warranty of 20MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 21Lesser General Public License for more details. 22 23You should have received a copy of the GNU Lesser General Public 24License along with this Module; if not, write to the Free Software 25Foundation, 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 34namespace Eigen { 35 36namespace internal { 37 38template<typename T> inline T amd_flip(const T& i) { return -i-2; } 39template<typename T> inline T amd_unflip(const T& i) { return i<0 ? amd_flip(i) : i; } 40template<typename T0, typename T1> inline bool amd_marked(const T0* w, const T1& j) { return w[j]<0; } 41template<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 */ 44template<typename Index> 45static 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 */ 59template<typename Index> 60Index 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 */ 90template<typename Scalar, typename Index> 91void 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