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