00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 #ifndef EIGEN_SPARSE_AMD_H
00049 #define EIGEN_SPARSE_AMD_H
00050
00051 namespace internal {
00052
00053
00054 #define CS_FLIP(i) (-(i)-2)
00055 #define CS_UNFLIP(i) (((i) < 0) ? CS_FLIP(i) : (i))
00056 #define CS_MARKED(w,j) (w[j] < 0)
00057 #define CS_MARK(w,j) { w[j] = CS_FLIP (w[j]); }
00058
00059
00060 template<typename Index>
00061 static int cs_wclear (Index mark, Index lemax, Index *w, Index n)
00062 {
00063 Index k;
00064 if(mark < 2 || (mark + lemax < 0))
00065 {
00066 for(k = 0; k < n; k++)
00067 if(w[k] != 0)
00068 w[k] = 1;
00069 mark = 2;
00070 }
00071 return (mark);
00072 }
00073
00074
00075 template<typename Index>
00076 Index cs_tdfs(Index j, Index k, Index *head, const Index *next, Index *post, Index *stack)
00077 {
00078 int i, p, top = 0;
00079 if(!head || !next || !post || !stack) return (-1);
00080 stack[0] = j;
00081 while (top >= 0)
00082 {
00083 p = stack[top];
00084 i = head[p];
00085 if(i == -1)
00086 {
00087 top--;
00088 post[k++] = p;
00089 }
00090 else
00091 {
00092 head[p] = next[i];
00093 stack[++top] = i;
00094 }
00095 }
00096 return k;
00097 }
00098
00099
00105 template<typename Scalar, typename Index>
00106 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,Index>& C, PermutationMatrix<Dynamic>& perm)
00107 {
00108 typedef SparseMatrix<Scalar,ColMajor,Index> CCS;
00109
00110 int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
00111 k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
00112 ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t;
00113 unsigned int h;
00114
00115 Index n = C.cols();
00116 dense = std::max<Index> (16, 10 * sqrt ((double) n));
00117 dense = std::min<Index> (n-2, dense);
00118
00119 Index cnz = C.nonZeros();
00120 perm.resize(n+1);
00121 t = cnz + cnz/5 + 2*n;
00122 C.resizeNonZeros(t);
00123
00124 Index* W = new Index[8*(n+1)];
00125 Index* len = W;
00126 Index* nv = W + (n+1);
00127 Index* next = W + 2*(n+1);
00128 Index* head = W + 3*(n+1);
00129 Index* elen = W + 4*(n+1);
00130 Index* degree = W + 5*(n+1);
00131 Index* w = W + 6*(n+1);
00132 Index* hhead = W + 7*(n+1);
00133 Index* last = perm.indices().data();
00134
00135
00136 Index* Cp = C._outerIndexPtr();
00137 Index* Ci = C._innerIndexPtr();
00138 for(k = 0; k < n; k++)
00139 len[k] = Cp[k+1] - Cp[k];
00140 len[n] = 0;
00141 nzmax = t;
00142
00143 for(i = 0; i <= n; i++)
00144 {
00145 head[i] = -1;
00146 last[i] = -1;
00147 next[i] = -1;
00148 hhead[i] = -1;
00149 nv[i] = 1;
00150 w[i] = 1;
00151 elen[i] = 0;
00152 degree[i] = len[i];
00153 }
00154 mark = cs_wclear (0, 0, w, n);
00155 elen[n] = -2;
00156 Cp[n] = -1;
00157 w[n] = 0;
00158
00159
00160 for(i = 0; i < n; i++)
00161 {
00162 d = degree[i];
00163 if(d == 0)
00164 {
00165 elen[i] = -2;
00166 nel++;
00167 Cp[i] = -1;
00168 w[i] = 0;
00169 }
00170 else if(d > dense)
00171 {
00172 nv[i] = 0;
00173 elen[i] = -1;
00174 nel++;
00175 Cp[i] = CS_FLIP (n);
00176 nv[n]++;
00177 }
00178 else
00179 {
00180 if(head[d] != -1) last[head[d]] = i;
00181 next[i] = head[d];
00182 head[d] = i;
00183 }
00184 }
00185
00186 while (nel < n)
00187 {
00188
00189 for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
00190 if(next[k] != -1) last[next[k]] = -1;
00191 head[mindeg] = next[k];
00192 elenk = elen[k];
00193 nvk = nv[k];
00194 nel += nvk;
00195
00196
00197 if(elenk > 0 && cnz + mindeg >= nzmax)
00198 {
00199 for(j = 0; j < n; j++)
00200 {
00201 if((p = Cp[j]) >= 0)
00202 {
00203 Cp[j] = Ci[p];
00204 Ci[p] = CS_FLIP (j);
00205 }
00206 }
00207 for(q = 0, p = 0; p < cnz; )
00208 {
00209 if((j = CS_FLIP (Ci[p++])) >= 0)
00210 {
00211 Ci[q] = Cp[j];
00212 Cp[j] = q++;
00213 for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
00214 }
00215 }
00216 cnz = q;
00217 }
00218
00219
00220 dk = 0;
00221 nv[k] = -nvk;
00222 p = Cp[k];
00223 pk1 = (elenk == 0) ? p : cnz;
00224 pk2 = pk1;
00225 for(k1 = 1; k1 <= elenk + 1; k1++)
00226 {
00227 if(k1 > elenk)
00228 {
00229 e = k;
00230 pj = p;
00231 ln = len[k] - elenk;
00232 }
00233 else
00234 {
00235 e = Ci[p++];
00236 pj = Cp[e];
00237 ln = len[e];
00238 }
00239 for(k2 = 1; k2 <= ln; k2++)
00240 {
00241 i = Ci[pj++];
00242 if((nvi = nv[i]) <= 0) continue;
00243 dk += nvi;
00244 nv[i] = -nvi;
00245 Ci[pk2++] = i;
00246 if(next[i] != -1) last[next[i]] = last[i];
00247 if(last[i] != -1)
00248 {
00249 next[last[i]] = next[i];
00250 }
00251 else
00252 {
00253 head[degree[i]] = next[i];
00254 }
00255 }
00256 if(e != k)
00257 {
00258 Cp[e] = CS_FLIP (k);
00259 w[e] = 0;
00260 }
00261 }
00262 if(elenk != 0) cnz = pk2;
00263 degree[k] = dk;
00264 Cp[k] = pk1;
00265 len[k] = pk2 - pk1;
00266 elen[k] = -2;
00267
00268
00269 mark = cs_wclear (mark, lemax, w, n);
00270 for(pk = pk1; pk < pk2; pk++)
00271 {
00272 i = Ci[pk];
00273 if((eln = elen[i]) <= 0) continue;
00274 nvi = -nv[i];
00275 wnvi = mark - nvi;
00276 for(p = Cp[i]; p <= Cp[i] + eln - 1; p++)
00277 {
00278 e = Ci[p];
00279 if(w[e] >= mark)
00280 {
00281 w[e] -= nvi;
00282 }
00283 else if(w[e] != 0)
00284 {
00285 w[e] = degree[e] + wnvi;
00286 }
00287 }
00288 }
00289
00290
00291 for(pk = pk1; pk < pk2; pk++)
00292 {
00293 i = Ci[pk];
00294 p1 = Cp[i];
00295 p2 = p1 + elen[i] - 1;
00296 pn = p1;
00297 for(h = 0, d = 0, p = p1; p <= p2; p++)
00298 {
00299 e = Ci[p];
00300 if(w[e] != 0)
00301 {
00302 dext = w[e] - mark;
00303 if(dext > 0)
00304 {
00305 d += dext;
00306 Ci[pn++] = e;
00307 h += e;
00308 }
00309 else
00310 {
00311 Cp[e] = CS_FLIP (k);
00312 w[e] = 0;
00313 }
00314 }
00315 }
00316 elen[i] = pn - p1 + 1;
00317 p3 = pn;
00318 p4 = p1 + len[i];
00319 for(p = p2 + 1; p < p4; p++)
00320 {
00321 j = Ci[p];
00322 if((nvj = nv[j]) <= 0) continue;
00323 d += nvj;
00324 Ci[pn++] = j;
00325 h += j;
00326 }
00327 if(d == 0)
00328 {
00329 Cp[i] = CS_FLIP (k);
00330 nvi = -nv[i];
00331 dk -= nvi;
00332 nvk += nvi;
00333 nel += nvi;
00334 nv[i] = 0;
00335 elen[i] = -1;
00336 }
00337 else
00338 {
00339 degree[i] = std::min<Index> (degree[i], d);
00340 Ci[pn] = Ci[p3];
00341 Ci[p3] = Ci[p1];
00342 Ci[p1] = k;
00343 len[i] = pn - p1 + 1;
00344 h %= n;
00345 next[i] = hhead[h];
00346 hhead[h] = i;
00347 last[i] = h;
00348 }
00349 }
00350 degree[k] = dk;
00351 lemax = std::max<Index>(lemax, dk);
00352 mark = cs_wclear (mark+lemax, lemax, w, n);
00353
00354
00355 for(pk = pk1; pk < pk2; pk++)
00356 {
00357 i = Ci[pk];
00358 if(nv[i] >= 0) continue;
00359 h = last[i];
00360 i = hhead[h];
00361 hhead[h] = -1;
00362 for(; i != -1 && next[i] != -1; i = next[i], mark++)
00363 {
00364 ln = len[i];
00365 eln = elen[i];
00366 for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
00367 jlast = i;
00368 for(j = next[i]; j != -1; )
00369 {
00370 ok = (len[j] == ln) && (elen[j] == eln);
00371 for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
00372 {
00373 if(w[Ci[p]] != mark) ok = 0;
00374 }
00375 if(ok)
00376 {
00377 Cp[j] = CS_FLIP (i);
00378 nv[i] += nv[j];
00379 nv[j] = 0;
00380 elen[j] = -1;
00381 j = next[j];
00382 next[jlast] = j;
00383 }
00384 else
00385 {
00386 jlast = j;
00387 j = next[j];
00388 }
00389 }
00390 }
00391 }
00392
00393
00394 for(p = pk1, pk = pk1; pk < pk2; pk++)
00395 {
00396 i = Ci[pk];
00397 if((nvi = -nv[i]) <= 0) continue;
00398 nv[i] = nvi;
00399 d = degree[i] + dk - nvi;
00400 d = std::min<Index> (d, n - nel - nvi);
00401 if(head[d] != -1) last[head[d]] = i;
00402 next[i] = head[d];
00403 last[i] = -1;
00404 head[d] = i;
00405 mindeg = std::min<Index> (mindeg, d);
00406 degree[i] = d;
00407 Ci[p++] = i;
00408 }
00409 nv[k] = nvk;
00410 if((len[k] = p-pk1) == 0)
00411 {
00412 Cp[k] = -1;
00413 w[k] = 0;
00414 }
00415 if(elenk != 0) cnz = p;
00416 }
00417
00418
00419 for(i = 0; i < n; i++) Cp[i] = CS_FLIP (Cp[i]);
00420 for(j = 0; j <= n; j++) head[j] = -1;
00421 for(j = n; j >= 0; j--)
00422 {
00423 if(nv[j] > 0) continue;
00424 next[j] = head[Cp[j]];
00425 head[Cp[j]] = j;
00426 }
00427 for(e = n; e >= 0; e--)
00428 {
00429 if(nv[e] <= 0) continue;
00430 if(Cp[e] != -1)
00431 {
00432 next[e] = head[Cp[e]];
00433 head[Cp[e]] = e;
00434 }
00435 }
00436 for(k = 0, i = 0; i <= n; i++)
00437 {
00438 if(Cp[i] == -1) k = cs_tdfs (i, k, head, next, perm.indices().data(), w);
00439 }
00440
00441 perm.indices().conservativeResize(n);
00442
00443 delete[] W;
00444 }
00445
00446 }
00447
00448 #endif // EIGEN_SPARSE_AMD_H