$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2010 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 00006 /* 00007 00008 NOTE: this routine has been adapted from the CSparse library: 00009 00010 Copyright (c) 2006, Timothy A. Davis. 00011 http://www.cise.ufl.edu/research/sparse/CSparse 00012 00013 CSparse is free software; you can redistribute it and/or 00014 modify it under the terms of the GNU Lesser General Public 00015 License as published by the Free Software Foundation; either 00016 version 2.1 of the License, or (at your option) any later version. 00017 00018 CSparse is distributed in the hope that it will be useful, 00019 but WITHOUT ANY WARRANTY; without even the implied warranty of 00020 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00021 Lesser General Public License for more details. 00022 00023 You should have received a copy of the GNU Lesser General Public 00024 License along with this Module; if not, write to the Free Software 00025 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 00026 00027 */ 00028 00029 #include "../Core/util/NonMPL2.h" 00030 00031 #ifndef EIGEN_SPARSE_AMD_H 00032 #define EIGEN_SPARSE_AMD_H 00033 00034 namespace Eigen { 00035 00036 namespace internal { 00037 00038 template<typename T> inline T amd_flip(const T& i) { return -i-2; } 00039 template<typename T> inline T amd_unflip(const T& i) { return i<0 ? amd_flip(i) : i; } 00040 template<typename T0, typename T1> inline bool amd_marked(const T0* w, const T1& j) { return w[j]<0; } 00041 template<typename T0, typename T1> inline void amd_mark(const T0* w, const T1& j) { return w[j] = amd_flip(w[j]); } 00042 00043 /* clear w */ 00044 template<typename Index> 00045 static int cs_wclear (Index mark, Index lemax, Index *w, Index n) 00046 { 00047 Index k; 00048 if(mark < 2 || (mark + lemax < 0)) 00049 { 00050 for(k = 0; k < n; k++) 00051 if(w[k] != 0) 00052 w[k] = 1; 00053 mark = 2; 00054 } 00055 return (mark); /* at this point, w[0..n-1] < mark holds */ 00056 } 00057 00058 /* depth-first search and postorder of a tree rooted at node j */ 00059 template<typename Index> 00060 Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack) 00061 { 00062 int i, p, top = 0; 00063 if(!head || !next || !post || !stack) return (-1); /* check inputs */ 00064 stack[0] = j; /* place j on the stack */ 00065 while (top >= 0) /* while (stack is not empty) */ 00066 { 00067 p = stack[top]; /* p = top of stack */ 00068 i = head[p]; /* i = youngest child of p */ 00069 if(i == -1) 00070 { 00071 top--; /* p has no unordered children left */ 00072 post[k++] = p; /* node p is the kth postordered node */ 00073 } 00074 else 00075 { 00076 head[p] = next[i]; /* remove i from children of p */ 00077 stack[++top] = i; /* start dfs on child node i */ 00078 } 00079 } 00080 return k; 00081 } 00082 00083 00090 template<typename Scalar, typename Index> 00091 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic,Dynamic,Index>& perm) 00092 { 00093 using std::sqrt; 00094 00095 int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1, 00096 k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi, 00097 ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t; 00098 unsigned int h; 00099 00100 Index n = C.cols(); 00101 dense = std::max<Index> (16, Index(10 * sqrt(double(n)))); /* find dense threshold */ 00102 dense = std::min<Index> (n-2, dense); 00103 00104 Index cnz = C.nonZeros(); 00105 perm.resize(n+1); 00106 t = cnz + cnz/5 + 2*n; /* add elbow room to C */ 00107 C.resizeNonZeros(t); 00108 00109 Index* W = new Index[8*(n+1)]; /* get workspace */ 00110 Index* len = W; 00111 Index* nv = W + (n+1); 00112 Index* next = W + 2*(n+1); 00113 Index* head = W + 3*(n+1); 00114 Index* elen = W + 4*(n+1); 00115 Index* degree = W + 5*(n+1); 00116 Index* w = W + 6*(n+1); 00117 Index* hhead = W + 7*(n+1); 00118 Index* last = perm.indices().data(); /* use P as workspace for last */ 00119 00120 /* --- Initialize quotient graph ---------------------------------------- */ 00121 Index* Cp = C.outerIndexPtr(); 00122 Index* Ci = C.innerIndexPtr(); 00123 for(k = 0; k < n; k++) 00124 len[k] = Cp[k+1] - Cp[k]; 00125 len[n] = 0; 00126 nzmax = t; 00127 00128 for(i = 0; i <= n; i++) 00129 { 00130 head[i] = -1; // degree list i is empty 00131 last[i] = -1; 00132 next[i] = -1; 00133 hhead[i] = -1; // hash list i is empty 00134 nv[i] = 1; // node i is just one node 00135 w[i] = 1; // node i is alive 00136 elen[i] = 0; // Ek of node i is empty 00137 degree[i] = len[i]; // degree of node i 00138 } 00139 mark = internal::cs_wclear<Index>(0, 0, w, n); /* clear w */ 00140 elen[n] = -2; /* n is a dead element */ 00141 Cp[n] = -1; /* n is a root of assembly tree */ 00142 w[n] = 0; /* n is a dead element */ 00143 00144 /* --- Initialize degree lists ------------------------------------------ */ 00145 for(i = 0; i < n; i++) 00146 { 00147 bool has_diag = false; 00148 for(p = Cp[i]; p<Cp[i+1]; ++p) 00149 if(Ci[p]==i) 00150 { 00151 has_diag = true; 00152 break; 00153 } 00154 00155 d = degree[i]; 00156 if(d == 1) /* node i is empty */ 00157 { 00158 elen[i] = -2; /* element i is dead */ 00159 nel++; 00160 Cp[i] = -1; /* i is a root of assembly tree */ 00161 w[i] = 0; 00162 } 00163 else if(d > dense || !has_diag) /* node i is dense or has no structural diagonal element */ 00164 { 00165 nv[i] = 0; /* absorb i into element n */ 00166 elen[i] = -1; /* node i is dead */ 00167 nel++; 00168 Cp[i] = amd_flip (n); 00169 nv[n]++; 00170 } 00171 else 00172 { 00173 if(head[d] != -1) last[head[d]] = i; 00174 next[i] = head[d]; /* put node i in degree list d */ 00175 head[d] = i; 00176 } 00177 } 00178 00179 elen[n] = -2; /* n is a dead element */ 00180 Cp[n] = -1; /* n is a root of assembly tree */ 00181 w[n] = 0; /* n is a dead element */ 00182 00183 while (nel < n) /* while (selecting pivots) do */ 00184 { 00185 /* --- Select node of minimum approximate degree -------------------- */ 00186 for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {} 00187 if(next[k] != -1) last[next[k]] = -1; 00188 head[mindeg] = next[k]; /* remove k from degree list */ 00189 elenk = elen[k]; /* elenk = |Ek| */ 00190 nvk = nv[k]; /* # of nodes k represents */ 00191 nel += nvk; /* nv[k] nodes of A eliminated */ 00192 00193 /* --- Garbage collection ------------------------------------------- */ 00194 if(elenk > 0 && cnz + mindeg >= nzmax) 00195 { 00196 for(j = 0; j < n; j++) 00197 { 00198 if((p = Cp[j]) >= 0) /* j is a live node or element */ 00199 { 00200 Cp[j] = Ci[p]; /* save first entry of object */ 00201 Ci[p] = amd_flip (j); /* first entry is now amd_flip(j) */ 00202 } 00203 } 00204 for(q = 0, p = 0; p < cnz; ) /* scan all of memory */ 00205 { 00206 if((j = amd_flip (Ci[p++])) >= 0) /* found object j */ 00207 { 00208 Ci[q] = Cp[j]; /* restore first entry of object */ 00209 Cp[j] = q++; /* new pointer to object j */ 00210 for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++]; 00211 } 00212 } 00213 cnz = q; /* Ci[cnz...nzmax-1] now free */ 00214 } 00215 00216 /* --- Construct new element ---------------------------------------- */ 00217 dk = 0; 00218 nv[k] = -nvk; /* flag k as in Lk */ 00219 p = Cp[k]; 00220 pk1 = (elenk == 0) ? p : cnz; /* do in place if elen[k] == 0 */ 00221 pk2 = pk1; 00222 for(k1 = 1; k1 <= elenk + 1; k1++) 00223 { 00224 if(k1 > elenk) 00225 { 00226 e = k; /* search the nodes in k */ 00227 pj = p; /* list of nodes starts at Ci[pj]*/ 00228 ln = len[k] - elenk; /* length of list of nodes in k */ 00229 } 00230 else 00231 { 00232 e = Ci[p++]; /* search the nodes in e */ 00233 pj = Cp[e]; 00234 ln = len[e]; /* length of list of nodes in e */ 00235 } 00236 for(k2 = 1; k2 <= ln; k2++) 00237 { 00238 i = Ci[pj++]; 00239 if((nvi = nv[i]) <= 0) continue; /* node i dead, or seen */ 00240 dk += nvi; /* degree[Lk] += size of node i */ 00241 nv[i] = -nvi; /* negate nv[i] to denote i in Lk*/ 00242 Ci[pk2++] = i; /* place i in Lk */ 00243 if(next[i] != -1) last[next[i]] = last[i]; 00244 if(last[i] != -1) /* remove i from degree list */ 00245 { 00246 next[last[i]] = next[i]; 00247 } 00248 else 00249 { 00250 head[degree[i]] = next[i]; 00251 } 00252 } 00253 if(e != k) 00254 { 00255 Cp[e] = amd_flip (k); /* absorb e into k */ 00256 w[e] = 0; /* e is now a dead element */ 00257 } 00258 } 00259 if(elenk != 0) cnz = pk2; /* Ci[cnz...nzmax] is free */ 00260 degree[k] = dk; /* external degree of k - |Lk\i| */ 00261 Cp[k] = pk1; /* element k is in Ci[pk1..pk2-1] */ 00262 len[k] = pk2 - pk1; 00263 elen[k] = -2; /* k is now an element */ 00264 00265 /* --- Find set differences ----------------------------------------- */ 00266 mark = internal::cs_wclear<Index>(mark, lemax, w, n); /* clear w if necessary */ 00267 for(pk = pk1; pk < pk2; pk++) /* scan 1: find |Le\Lk| */ 00268 { 00269 i = Ci[pk]; 00270 if((eln = elen[i]) <= 0) continue;/* skip if elen[i] empty */ 00271 nvi = -nv[i]; /* nv[i] was negated */ 00272 wnvi = mark - nvi; 00273 for(p = Cp[i]; p <= Cp[i] + eln - 1; p++) /* scan Ei */ 00274 { 00275 e = Ci[p]; 00276 if(w[e] >= mark) 00277 { 00278 w[e] -= nvi; /* decrement |Le\Lk| */ 00279 } 00280 else if(w[e] != 0) /* ensure e is a live element */ 00281 { 00282 w[e] = degree[e] + wnvi; /* 1st time e seen in scan 1 */ 00283 } 00284 } 00285 } 00286 00287 /* --- Degree update ------------------------------------------------ */ 00288 for(pk = pk1; pk < pk2; pk++) /* scan2: degree update */ 00289 { 00290 i = Ci[pk]; /* consider node i in Lk */ 00291 p1 = Cp[i]; 00292 p2 = p1 + elen[i] - 1; 00293 pn = p1; 00294 for(h = 0, d = 0, p = p1; p <= p2; p++) /* scan Ei */ 00295 { 00296 e = Ci[p]; 00297 if(w[e] != 0) /* e is an unabsorbed element */ 00298 { 00299 dext = w[e] - mark; /* dext = |Le\Lk| */ 00300 if(dext > 0) 00301 { 00302 d += dext; /* sum up the set differences */ 00303 Ci[pn++] = e; /* keep e in Ei */ 00304 h += e; /* compute the hash of node i */ 00305 } 00306 else 00307 { 00308 Cp[e] = amd_flip (k); /* aggressive absorb. e->k */ 00309 w[e] = 0; /* e is a dead element */ 00310 } 00311 } 00312 } 00313 elen[i] = pn - p1 + 1; /* elen[i] = |Ei| */ 00314 p3 = pn; 00315 p4 = p1 + len[i]; 00316 for(p = p2 + 1; p < p4; p++) /* prune edges in Ai */ 00317 { 00318 j = Ci[p]; 00319 if((nvj = nv[j]) <= 0) continue; /* node j dead or in Lk */ 00320 d += nvj; /* degree(i) += |j| */ 00321 Ci[pn++] = j; /* place j in node list of i */ 00322 h += j; /* compute hash for node i */ 00323 } 00324 if(d == 0) /* check for mass elimination */ 00325 { 00326 Cp[i] = amd_flip (k); /* absorb i into k */ 00327 nvi = -nv[i]; 00328 dk -= nvi; /* |Lk| -= |i| */ 00329 nvk += nvi; /* |k| += nv[i] */ 00330 nel += nvi; 00331 nv[i] = 0; 00332 elen[i] = -1; /* node i is dead */ 00333 } 00334 else 00335 { 00336 degree[i] = std::min<Index> (degree[i], d); /* update degree(i) */ 00337 Ci[pn] = Ci[p3]; /* move first node to end */ 00338 Ci[p3] = Ci[p1]; /* move 1st el. to end of Ei */ 00339 Ci[p1] = k; /* add k as 1st element in of Ei */ 00340 len[i] = pn - p1 + 1; /* new len of adj. list of node i */ 00341 h %= n; /* finalize hash of i */ 00342 next[i] = hhead[h]; /* place i in hash bucket */ 00343 hhead[h] = i; 00344 last[i] = h; /* save hash of i in last[i] */ 00345 } 00346 } /* scan2 is done */ 00347 degree[k] = dk; /* finalize |Lk| */ 00348 lemax = std::max<Index>(lemax, dk); 00349 mark = internal::cs_wclear<Index>(mark+lemax, lemax, w, n); /* clear w */ 00350 00351 /* --- Supernode detection ------------------------------------------ */ 00352 for(pk = pk1; pk < pk2; pk++) 00353 { 00354 i = Ci[pk]; 00355 if(nv[i] >= 0) continue; /* skip if i is dead */ 00356 h = last[i]; /* scan hash bucket of node i */ 00357 i = hhead[h]; 00358 hhead[h] = -1; /* hash bucket will be empty */ 00359 for(; i != -1 && next[i] != -1; i = next[i], mark++) 00360 { 00361 ln = len[i]; 00362 eln = elen[i]; 00363 for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark; 00364 jlast = i; 00365 for(j = next[i]; j != -1; ) /* compare i with all j */ 00366 { 00367 ok = (len[j] == ln) && (elen[j] == eln); 00368 for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++) 00369 { 00370 if(w[Ci[p]] != mark) ok = 0; /* compare i and j*/ 00371 } 00372 if(ok) /* i and j are identical */ 00373 { 00374 Cp[j] = amd_flip (i); /* absorb j into i */ 00375 nv[i] += nv[j]; 00376 nv[j] = 0; 00377 elen[j] = -1; /* node j is dead */ 00378 j = next[j]; /* delete j from hash bucket */ 00379 next[jlast] = j; 00380 } 00381 else 00382 { 00383 jlast = j; /* j and i are different */ 00384 j = next[j]; 00385 } 00386 } 00387 } 00388 } 00389 00390 /* --- Finalize new element------------------------------------------ */ 00391 for(p = pk1, pk = pk1; pk < pk2; pk++) /* finalize Lk */ 00392 { 00393 i = Ci[pk]; 00394 if((nvi = -nv[i]) <= 0) continue;/* skip if i is dead */ 00395 nv[i] = nvi; /* restore nv[i] */ 00396 d = degree[i] + dk - nvi; /* compute external degree(i) */ 00397 d = std::min<Index> (d, n - nel - nvi); 00398 if(head[d] != -1) last[head[d]] = i; 00399 next[i] = head[d]; /* put i back in degree list */ 00400 last[i] = -1; 00401 head[d] = i; 00402 mindeg = std::min<Index> (mindeg, d); /* find new minimum degree */ 00403 degree[i] = d; 00404 Ci[p++] = i; /* place i in Lk */ 00405 } 00406 nv[k] = nvk; /* # nodes absorbed into k */ 00407 if((len[k] = p-pk1) == 0) /* length of adj list of element k*/ 00408 { 00409 Cp[k] = -1; /* k is a root of the tree */ 00410 w[k] = 0; /* k is now a dead element */ 00411 } 00412 if(elenk != 0) cnz = p; /* free unused space in Lk */ 00413 } 00414 00415 /* --- Postordering ----------------------------------------------------- */ 00416 for(i = 0; i < n; i++) Cp[i] = amd_flip (Cp[i]);/* fix assembly tree */ 00417 for(j = 0; j <= n; j++) head[j] = -1; 00418 for(j = n; j >= 0; j--) /* place unordered nodes in lists */ 00419 { 00420 if(nv[j] > 0) continue; /* skip if j is an element */ 00421 next[j] = head[Cp[j]]; /* place j in list of its parent */ 00422 head[Cp[j]] = j; 00423 } 00424 for(e = n; e >= 0; e--) /* place elements in lists */ 00425 { 00426 if(nv[e] <= 0) continue; /* skip unless e is an element */ 00427 if(Cp[e] != -1) 00428 { 00429 next[e] = head[Cp[e]]; /* place e in list of its parent */ 00430 head[Cp[e]] = e; 00431 } 00432 } 00433 for(k = 0, i = 0; i <= n; i++) /* postorder the assembly tree */ 00434 { 00435 if(Cp[i] == -1) k = internal::cs_tdfs<Index>(i, k, head, next, perm.indices().data(), w); 00436 } 00437 00438 perm.indices().conservativeResize(n); 00439 00440 delete[] W; 00441 } 00442 00443 } // namespace internal 00444 00445 } // end namespace Eigen 00446 00447 #endif // EIGEN_SPARSE_AMD_H