1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
#include "matrix.hpp"
class MatrixSorter
{
const Ring *R;
int deg_ascending;
int ringorder_ascending;
int *sort_vals;
vec *sort_vecs;
int *sort_degs;
M2_arrayint result;
int sort_compare(int i, int j)
{
if (i == j) return 0;
vec v1 = sort_vecs[i];
vec v2 = sort_vecs[j];
if (v1 == NULL) return 1;
if (v2 == NULL) return -1;
if (deg_ascending != 0)
{
int d1 = sort_degs[i];
int d2 = sort_degs[j];
if (d1 > d2) return -deg_ascending;
if (d1 < d2) return deg_ascending;
}
int cmp = R->compare_vecs(v1, v2);
if (cmp > 0) return -ringorder_ascending;
if (cmp < 0) return ringorder_ascending;
return 0;
}
int sort_partition(int lo, int hi)
{
int pivot = sort_vals[lo];
int i = lo - 1;
int j = hi + 1;
for (;;)
{
do
{
j--;
}
while (sort_compare(sort_vals[j], pivot) < 0);
do
{
i++;
}
while (sort_compare(sort_vals[i], pivot) > 0);
if (i < j)
{
int tmp = sort_vals[j];
sort_vals[j] = sort_vals[i];
sort_vals[i] = tmp;
}
else
return j;
}
}
void sort_range(int lo, int hi)
{
if (lo < hi)
{
int q = sort_partition(lo, hi);
sort_range(lo, q);
sort_range(q + 1, hi);
}
}
public:
MatrixSorter(const Matrix *m, int degorder, int ringorder);
M2_arrayintOrNull value();
};
MatrixSorter::MatrixSorter(const Matrix *m, int degorder, int ringorder)
: deg_ascending(degorder), ringorder_ascending(ringorder)
{
R = m->get_ring();
int nelems = m->n_cols();
if (degorder != 0)
{
sort_degs = newarray_atomic(int, nelems);
for (int i = 0; i < nelems; i++)
sort_degs[i] = m->cols()->primary_degree(i);
}
result = M2_makearrayint(nelems);
sort_vals = result->array;
for (int i = 0; i < nelems; i++) sort_vals[i] = i;
sort_vecs = newarray(vec, nelems);
for (int i = 0; i < nelems; i++) sort_vecs[i] = m->elem(i);
}
M2_arrayintOrNull MatrixSorter::value()
{
sort_range(0, result->len - 1);
return result;
}
M2_arrayint Matrix::sort(int degorder, int ringorder) const
// Sort the columns of 'this': Place the column indices into 'result'.
// If degorder < 0, sort in descending degree order, if >0 ascending degree
// If ==0, or in the event that two columns have the same (simple) degree,
// use the ring order: ringorder > 0 means ascending, <0 means descending.
{
MatrixSorter sorter(this, degorder, ringorder);
return sorter.value();
}
// Local Variables:
// compile-command: "make -C $M2BUILDDIR/Macaulay2/e "
// indent-tabs-mode: nil
// End:
|
