48 #ifndef EIGEN_SPARSE_AMD_H
49 #define EIGEN_SPARSE_AMD_H
55 template<
typename T>
inline T
amd_flip(
const T& i) {
return -i-2; }
57 template<
typename T0,
typename T1>
inline bool amd_marked(
const T0* w,
const T1& j) {
return w[j]<0; }
58 template<
typename T0,
typename T1>
inline void amd_mark(
const T0* w,
const T1& j) {
return w[j] =
amd_flip(w[j]); }
61 template<
typename Index>
62 static int cs_wclear (Index mark, Index lemax, Index *w, Index n)
65 if(mark < 2 || (mark + lemax < 0))
67 for(k = 0; k < n; k++)
76 template<
typename Index>
77 Index
cs_tdfs(Index j, Index k, Index *head,
const Index *next, Index *post, Index *stack)
80 if(!head || !next || !post || !stack)
return (-1);
106 template<
typename Scalar,
typename Index>
112 int d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
113 k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
114 ok, nel = 0,
p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn,
q, t;
118 dense = std::max<Index> (16, Index(10 *
sqrt(
double(n))));
119 dense = std::min<Index> (n-2, dense);
123 t = cnz + cnz/5 + 2*n;
126 Index* W =
new Index[8*(n+1)];
128 Index* nv = W + (n+1);
129 Index* next = W + 2*(n+1);
130 Index* head = W + 3*(n+1);
131 Index* elen = W + 4*(n+1);
132 Index* degree = W + 5*(n+1);
133 Index* w = W + 6*(n+1);
134 Index* hhead = W + 7*(n+1);
135 Index* last = perm.
indices().data();
140 for(k = 0; k < n; k++)
141 len[k] = Cp[k+1] - Cp[k];
145 for(i = 0; i <= n; i++)
156 mark = internal::cs_wclear<Index>(0, 0, w, n);
162 for(i = 0; i < n; i++)
182 if(head[d] != -1) last[head[d]] = i;
191 for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
192 if(next[k] != -1) last[next[k]] = -1;
193 head[mindeg] = next[k];
199 if(elenk > 0 && cnz + mindeg >= nzmax)
201 for(j = 0; j < n; j++)
209 for(q = 0,
p = 0;
p < cnz; )
215 for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[
p++];
225 pk1 = (elenk == 0) ?
p : cnz;
227 for(k1 = 1; k1 <= elenk + 1; k1++)
241 for(k2 = 1; k2 <= ln; k2++)
244 if((nvi = nv[i]) <= 0)
continue;
248 if(next[i] != -1) last[next[i]] = last[i];
251 next[last[i]] = next[i];
255 head[degree[i]] = next[i];
264 if(elenk != 0) cnz = pk2;
271 mark = internal::cs_wclear<Index>(mark, lemax, w, n);
272 for(pk = pk1; pk < pk2; pk++)
275 if((eln = elen[i]) <= 0)
continue;
278 for(
p = Cp[i];
p <= Cp[i] + eln - 1;
p++)
287 w[e] = degree[e] + wnvi;
293 for(pk = pk1; pk < pk2; pk++)
297 p2 = p1 + elen[i] - 1;
299 for(h = 0, d = 0,
p = p1;
p <= p2;
p++)
318 elen[i] = pn - p1 + 1;
321 for(
p = p2 + 1;
p < p4;
p++)
324 if((nvj = nv[j]) <= 0)
continue;
341 degree[i] = std::min<Index> (degree[i], d);
345 len[i] = pn - p1 + 1;
353 lemax = std::max<Index>(lemax, dk);
354 mark = internal::cs_wclear<Index>(mark+lemax, lemax, w, n);
357 for(pk = pk1; pk < pk2; pk++)
360 if(nv[i] >= 0)
continue;
364 for(; i != -1 && next[i] != -1; i = next[i], mark++)
368 for(
p = Cp[i]+1;
p <= Cp[i] + ln-1;
p++) w[Ci[
p]] = mark;
370 for(j = next[i]; j != -1; )
372 ok = (len[j] == ln) && (elen[j] == eln);
373 for(
p = Cp[j] + 1; ok &&
p <= Cp[j] + ln - 1; p++)
375 if(w[Ci[p]] != mark) ok = 0;
396 for(
p = pk1, pk = pk1; pk < pk2; pk++)
399 if((nvi = -nv[i]) <= 0)
continue;
401 d = degree[i] + dk - nvi;
402 d = std::min<Index> (d, n - nel - nvi);
403 if(head[d] != -1) last[head[d]] = i;
407 mindeg = std::min<Index> (mindeg, d);
412 if((len[k] =
p-pk1) == 0)
417 if(elenk != 0) cnz =
p;
421 for(i = 0; i < n; i++) Cp[i] =
amd_flip (Cp[i]);
422 for(j = 0; j <= n; j++) head[j] = -1;
423 for(j = n; j >= 0; j--)
425 if(nv[j] > 0)
continue;
426 next[j] = head[Cp[j]];
429 for(e = n; e >= 0; e--)
431 if(nv[e] <= 0)
continue;
434 next[e] = head[Cp[e]];
438 for(k = 0, i = 0; i <= n; i++)
440 if(Cp[i] == -1) k = internal::cs_tdfs<Index>(i, k, head, next, perm.
indices().data(), w);
443 perm.
indices().conservativeResize(n);
452 #endif // EIGEN_SPARSE_AMD_H