29 #include "../Core/util/NonMPL2.h" 31 #ifndef EIGEN_SPARSE_AMD_H 32 #define EIGEN_SPARSE_AMD_H 38 template<
typename T>
inline T amd_flip(
const T& i) {
return -i-2; }
39 template<
typename T>
inline T amd_unflip(
const T& i) {
return i<0 ? amd_flip(i) : i; }
40 template<
typename T0,
typename T1>
inline bool amd_marked(
const T0* w,
const T1& j) {
return w[j]<0; }
41 template<
typename T0,
typename T1>
inline void amd_mark(
const T0* w,
const T1& j) {
return w[j] = amd_flip(w[j]); }
44 template<
typename StorageIndex>
45 static StorageIndex cs_wclear (StorageIndex mark, StorageIndex lemax, StorageIndex *w, StorageIndex n)
48 if(mark < 2 || (mark + lemax < 0))
50 for(k = 0; k < n; k++)
59 template<
typename StorageIndex>
60 StorageIndex cs_tdfs(StorageIndex j, StorageIndex k, StorageIndex *head,
const StorageIndex *next, StorageIndex *post, StorageIndex *stack)
62 StorageIndex i, p, top = 0;
63 if(!head || !next || !post || !stack)
return (-1);
93 template<
typename Scalar,
typename StorageIndex>
94 void minimum_degree_ordering(SparseMatrix<Scalar,ColMajor,StorageIndex>& C, PermutationMatrix<Dynamic,Dynamic,StorageIndex>& perm)
98 StorageIndex d, dk, dext, lemax = 0, e, elenk, eln, i, j, k, k1,
99 k2, k3, jlast, ln, dense, nzmax, mindeg = 0, nvi, nvj, nvk, mark, wnvi,
100 ok, nel = 0, p, p1, p2, p3, p4, pj, pk, pk1, pk2, pn, q, t, h;
102 StorageIndex n = StorageIndex(C.cols());
103 dense = std::max<StorageIndex> (16, StorageIndex(10 *
sqrt(
double(n))));
104 dense = (std::min)(n-2, dense);
106 StorageIndex cnz = StorageIndex(C.nonZeros());
108 t = cnz + cnz/5 + 2*n;
112 ei_declare_aligned_stack_constructed_variable(StorageIndex,W,8*(n+1),0);
113 StorageIndex* len = W;
114 StorageIndex* nv = W + (n+1);
115 StorageIndex* next = W + 2*(n+1);
116 StorageIndex* head = W + 3*(n+1);
117 StorageIndex* elen = W + 4*(n+1);
118 StorageIndex* degree = W + 5*(n+1);
119 StorageIndex* w = W + 6*(n+1);
120 StorageIndex* hhead = W + 7*(n+1);
121 StorageIndex* last = perm.indices().data();
124 StorageIndex* Cp = C.outerIndexPtr();
125 StorageIndex* Ci = C.innerIndexPtr();
126 for(k = 0; k < n; k++)
127 len[k] = Cp[k+1] - Cp[k];
131 for(i = 0; i <= n; i++)
142 mark = internal::cs_wclear<StorageIndex>(0, 0, w, n);
145 for(i = 0; i < n; i++)
147 bool has_diag =
false;
148 for(p = Cp[i]; p<Cp[i+1]; ++p)
156 if(d == 1 && has_diag)
163 else if(d > dense || !has_diag)
168 Cp[i] = amd_flip (n);
173 if(head[d] != -1) last[head[d]] = i;
186 for(k = -1; mindeg < n && (k = head[mindeg]) == -1; mindeg++) {}
187 if(next[k] != -1) last[next[k]] = -1;
188 head[mindeg] = next[k];
194 if(elenk > 0 && cnz + mindeg >= nzmax)
196 for(j = 0; j < n; j++)
201 Ci[p] = amd_flip (j);
204 for(q = 0, p = 0; p < cnz; )
206 if((j = amd_flip (Ci[p++])) >= 0)
210 for(k3 = 0; k3 < len[j]-1; k3++) Ci[q++] = Ci[p++];
220 pk1 = (elenk == 0) ? p : cnz;
222 for(k1 = 1; k1 <= elenk + 1; k1++)
236 for(k2 = 1; k2 <= ln; k2++)
239 if((nvi = nv[i]) <= 0)
continue;
243 if(next[i] != -1) last[next[i]] = last[i];
246 next[last[i]] = next[i];
250 head[degree[i]] = next[i];
255 Cp[e] = amd_flip (k);
259 if(elenk != 0) cnz = pk2;
266 mark = internal::cs_wclear<StorageIndex>(mark, lemax, w, n);
267 for(pk = pk1; pk < pk2; pk++)
270 if((eln = elen[i]) <= 0)
continue;
273 for(p = Cp[i]; p <= Cp[i] + eln - 1; p++)
282 w[e] = degree[e] + wnvi;
288 for(pk = pk1; pk < pk2; pk++)
292 p2 = p1 + elen[i] - 1;
294 for(h = 0, d = 0, p = p1; p <= p2; p++)
308 Cp[e] = amd_flip (k);
313 elen[i] = pn - p1 + 1;
316 for(p = p2 + 1; p < p4; p++)
319 if((nvj = nv[j]) <= 0)
continue;
326 Cp[i] = amd_flip (k);
336 degree[i] = std::min<StorageIndex> (degree[i], d);
340 len[i] = pn - p1 + 1;
348 lemax = std::max<StorageIndex>(lemax, dk);
349 mark = internal::cs_wclear<StorageIndex>(mark+lemax, lemax, w, n);
352 for(pk = pk1; pk < pk2; pk++)
355 if(nv[i] >= 0)
continue;
359 for(; i != -1 && next[i] != -1; i = next[i], mark++)
363 for(p = Cp[i]+1; p <= Cp[i] + ln-1; p++) w[Ci[p]] = mark;
365 for(j = next[i]; j != -1; )
367 ok = (len[j] == ln) && (elen[j] == eln);
368 for(p = Cp[j] + 1; ok && p <= Cp[j] + ln - 1; p++)
370 if(w[Ci[p]] != mark) ok = 0;
374 Cp[j] = amd_flip (i);
391 for(p = pk1, pk = pk1; pk < pk2; pk++)
394 if((nvi = -nv[i]) <= 0)
continue;
396 d = degree[i] + dk - nvi;
397 d = std::min<StorageIndex> (d, n - nel - nvi);
398 if(head[d] != -1) last[head[d]] = i;
402 mindeg = std::min<StorageIndex> (mindeg, d);
407 if((len[k] = p-pk1) == 0)
412 if(elenk != 0) cnz = p;
416 for(i = 0; i < n; i++) Cp[i] = amd_flip (Cp[i]);
417 for(j = 0; j <= n; j++) head[j] = -1;
418 for(j = n; j >= 0; j--)
420 if(nv[j] > 0)
continue;
421 next[j] = head[Cp[j]];
424 for(e = n; e >= 0; e--)
426 if(nv[e] <= 0)
continue;
429 next[e] = head[Cp[e]];
433 for(k = 0, i = 0; i <= n; i++)
435 if(Cp[i] == -1) k = internal::cs_tdfs<StorageIndex>(i, k, head, next, perm.indices().data(), w);
438 perm.indices().conservativeResize(n);
445 #endif // EIGEN_SPARSE_AMD_H const Eigen::CwiseUnaryOp< Eigen::internal::scalar_sqrt_op< typename Derived::Scalar >, const Derived > sqrt(const Eigen::ArrayBase< Derived > &x)
Namespace containing all symbols from the Eigen library.
Definition: Core:287
Definition: Eigen_Colamd.h:50