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