Skip to content

Commit 4a89949

Browse files
fredrik-johanssonfredrik-johansson
authored
Matrix nullspace and rref refactoring
* Add gr_mat_nullspace_no_resize, gr_mat_nullspace_from_rref * Make fq_X_mat_nullspace, fq_X_mat_rref gr wrappers
1 parent e6562ca commit 4a89949

File tree

6 files changed

+111
-309
lines changed

6 files changed

+111
-309
lines changed

‎doc/source/gr_mat.rst

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -577,6 +577,23 @@ Nullspace
577577
basis over the fraction field. The result may be meaningless
578578
if the ring is not an integral domain.
579579

580+
.. function:: int gr_mat_nullspace_from_rref(gr_mat_t X, const gr_mat_t A, gr_srcptr Aden, slong rank, gr_ctx_t ctx)
581+
582+
Computes nullspace given the precomputed reduced row
583+
echelon form matrix *A* with rank *rank*.
584+
If *Aden* is not *NULL*, assume tha *A* has been multiplied
585+
by this common denominator as in the output of
586+
:func:`gr_mat_rref_den`.
587+
588+
.. function:: int gr_mat_nullspace_no_resize(slong * nullity, gr_mat_t X, const gr_mat_t A, gr_ctx_t ctx)
589+
590+
Similar to :func:`gr_mat_nullspace`, but does not resize the
591+
matrix *X*, instead zero-padding if needed and returning
592+
the nullity (number of basis columns) in a separate variable.
593+
The user must supply a output matrix *X* with at least as many
594+
columns as the nullity.
595+
596+
580597
Inverse and adjugate
581598
-------------------------------------------------------------------------------
582599

‎src/fmpz_mod_mat/nullspace.c

Lines changed: 6 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -11,71 +11,15 @@
1111
(at your option) any later version. See <https://www.gnu.org/licenses/>.
1212
*/
1313

14-
#include "fmpz.h"
15-
#include "fmpz_mod.h"
1614
#include "fmpz_mod_mat.h"
15+
#include "gr.h"
16+
#include "gr_mat.h"
1717

1818
slong fmpz_mod_mat_nullspace(fmpz_mod_mat_t X, const fmpz_mod_mat_t A, const fmpz_mod_ctx_t ctx)
1919
{
20-
slong i, j, k, m, n, rank, nullity;
21-
slong *p;
22-
slong *pivots;
23-
slong *nonpivots;
24-
fmpz_mod_mat_t tmp;
25-
26-
m = A->r;
27-
n = A->c;
28-
29-
p = flint_malloc(sizeof(slong) * FLINT_MAX(m, n));
30-
31-
fmpz_mod_mat_init_set(tmp, A, ctx);
32-
rank = fmpz_mod_mat_rref(tmp, tmp, ctx);
33-
nullity = n - rank;
34-
35-
fmpz_mod_mat_zero(X, ctx);
36-
37-
if (rank == 0)
38-
{
39-
for (i = 0; i < nullity; i++)
40-
fmpz_one(fmpz_mod_mat_entry(X, i, i));
41-
}
42-
else if (nullity)
43-
{
44-
pivots = p; /* length = rank */
45-
nonpivots = p + rank; /* length = nullity */
46-
47-
for (i = j = k = 0; i < rank; i++)
48-
{
49-
while (fmpz_is_zero(fmpz_mod_mat_entry(tmp, i, j)))
50-
{
51-
nonpivots[k] = j;
52-
k++;
53-
j++;
54-
}
55-
pivots[i] = j;
56-
j++;
57-
}
58-
while (k < nullity)
59-
{
60-
nonpivots[k] = j;
61-
k++;
62-
j++;
63-
}
64-
65-
for (i = 0; i < nullity; i++)
66-
{
67-
for (j = 0; j < rank; j++)
68-
{
69-
fmpz_mod_neg(fmpz_mod_mat_entry(X, pivots[j], i),
70-
fmpz_mod_mat_entry(tmp, j, nonpivots[i]), ctx);
71-
}
72-
73-
fmpz_one(fmpz_mod_mat_entry(X, nonpivots[i], i));
74-
}
75-
}
76-
77-
flint_free(p);
78-
fmpz_mod_mat_clear(tmp, ctx);
79-
20+
gr_ctx_t gr_ctx;
21+
slong nullity;
22+
_gr_ctx_init_fmpz_mod_from_ref(gr_ctx, ctx);
23+
GR_MUST_SUCCEED(gr_mat_nullspace_no_resize(&nullity, (gr_mat_struct *) X, (const gr_mat_struct *) A, gr_ctx));
8024
return nullity;
8125
}

‎src/fq_mat_templates/nullspace.c

Lines changed: 6 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -12,75 +12,18 @@
1212

1313
#ifdef T
1414

15+
#include "gr.h"
16+
#include "gr_mat.h"
1517
#include "templates.h"
1618

1719
slong
1820
TEMPLATE(T, mat_nullspace) (TEMPLATE(T, mat_t) X, const TEMPLATE(T, mat_t) A,
1921
const TEMPLATE(T, ctx_t) ctx)
2022
{
21-
slong i, j, k, m, n, rank, nullity;
22-
slong *p;
23-
slong *pivots;
24-
slong *nonpivots;
25-
TEMPLATE(T, mat_t) tmp;
26-
27-
m = A->r;
28-
n = A->c;
29-
30-
p = flint_malloc(sizeof(slong) * FLINT_MAX(m, n));
31-
32-
TEMPLATE(T, mat_init_set) (tmp, A, ctx);
33-
rank = TEMPLATE(T, mat_rref) (tmp, tmp, ctx);
34-
nullity = n - rank;
35-
36-
TEMPLATE(T, mat_zero) (X, ctx);
37-
38-
if (rank == 0)
39-
{
40-
for (i = 0; i < nullity; i++)
41-
TEMPLATE(T, one) (TEMPLATE(T, mat_entry) (X, i, i), ctx);
42-
}
43-
else if (nullity)
44-
{
45-
pivots = p; /* length = rank */
46-
nonpivots = p + rank; /* length = nullity */
47-
48-
for (i = j = k = 0; i < rank; i++)
49-
{
50-
while (TEMPLATE(T, is_zero)
51-
(TEMPLATE(T, mat_entry) (tmp, i, j), ctx))
52-
{
53-
nonpivots[k] = j;
54-
k++;
55-
j++;
56-
}
57-
pivots[i] = j;
58-
j++;
59-
}
60-
while (k < nullity)
61-
{
62-
nonpivots[k] = j;
63-
k++;
64-
j++;
65-
}
66-
67-
for (i = 0; i < nullity; i++)
68-
{
69-
for (j = 0; j < rank; j++)
70-
{
71-
TEMPLATE(T, neg) (TEMPLATE(T, mat_entry) (X, pivots[j], i),
72-
TEMPLATE(T, mat_entry) (tmp, j,
73-
nonpivots[i]), ctx);
74-
}
75-
76-
TEMPLATE(T, one) (TEMPLATE(T, mat_entry) (X, nonpivots[i], i),
77-
ctx);
78-
}
79-
}
80-
81-
flint_free(p);
82-
TEMPLATE(T, mat_clear) (tmp, ctx);
83-
23+
gr_ctx_t gr_ctx;
24+
slong nullity;
25+
TEMPLATE3(_gr_ctx_init, T, from_ref)(gr_ctx, ctx);
26+
GR_MUST_SUCCEED(gr_mat_nullspace_no_resize(&nullity, (gr_mat_struct *) X, (const gr_mat_struct *) A, gr_ctx));
8427
return nullity;
8528
}
8629

‎src/fq_mat_templates/rref.c

Lines changed: 6 additions & 139 deletions
Original file line numberDiff line numberDiff line change
@@ -12,151 +12,18 @@
1212

1313
#ifdef T
1414

15+
#include "gr.h"
16+
#include "gr_mat.h"
1517
#include "templates.h"
1618

17-
#include "perm.h"
18-
1919
slong
2020
TEMPLATE(T, mat_rref) (TEMPLATE(T, mat_t) B, const TEMPLATE(T, mat_t) A, const TEMPLATE(T, ctx_t) ctx)
2121
{
22-
slong i, j, k, n, rank;
23-
slong *pivots;
24-
slong *nonpivots;
25-
slong *P;
26-
TEMPLATE(T, struct) * e;
27-
TEMPLATE(T, mat_t) U, V;
28-
29-
TEMPLATE(T, mat_set)(B, A, ctx);
30-
31-
if (TEMPLATE(T, mat_is_zero)(B, ctx))
32-
return 0;
33-
34-
if (B->r == 1)
35-
{
36-
TEMPLATE(T, struct) * c;
37-
slong i, j;
38-
slong r = 0;
39-
40-
for (i = 0; i < B->c; i++)
41-
{
42-
c = TEMPLATE(T, mat_entry)(B, 0, i);
43-
if (!TEMPLATE(T, is_zero)(c, ctx))
44-
{
45-
r = 1;
46-
if (TEMPLATE(T, is_one)(c, ctx))
47-
break;
48-
49-
TEMPLATE(T, inv)(c, c, ctx);
50-
for (j = i + 1;j < B->c; j++)
51-
{
52-
TEMPLATE(T, mul)(TEMPLATE(T, mat_entry)(B, 0, j), TEMPLATE(T, mat_entry)(B, 0, j), c, ctx);
53-
}
54-
TEMPLATE(T, one)(c, ctx);
55-
break;
56-
}
57-
}
58-
return r;
59-
}
60-
61-
n = A->c;
62-
63-
P = _perm_init(TEMPLATE(T, mat_nrows) (B, ctx));
64-
rank = TEMPLATE(T, mat_lu) (P, B, 0, ctx);
65-
_perm_clear(P);
66-
67-
if (rank == 0)
68-
return rank;
69-
70-
/* Clear L */
71-
for (i = 0; i < B->r; i++)
72-
for (j = 0; j < FLINT_MIN(i, rank); j++)
73-
TEMPLATE(T, zero) (TEMPLATE(T, mat_entry) (B, i, j), ctx);
74-
75-
/* We now reorder U to proper upper triangular form U | V
76-
with U full-rank triangular, set V = U^(-1) V, and then
77-
put the column back in the original order.
78-
79-
An improvement for some matrices would be to compress V by
80-
discarding columns containing nothing but zeros. */
81-
82-
TEMPLATE(T, mat_init) (U, rank, rank, ctx);
83-
TEMPLATE(T, mat_init) (V, rank, n - rank, ctx);
84-
85-
pivots = flint_malloc(sizeof(slong) * rank);
86-
nonpivots = flint_malloc(sizeof(slong) * (n - rank));
87-
88-
for (i = j = k = 0; i < rank; i++)
89-
{
90-
while (TEMPLATE(T, is_zero) (TEMPLATE(T, mat_entry) (B, i, j), ctx))
91-
{
92-
nonpivots[k] = j;
93-
k++;
94-
j++;
95-
}
96-
pivots[i] = j;
97-
j++;
98-
}
99-
while (k < n - rank)
100-
{
101-
nonpivots[k] = j;
102-
k++;
103-
j++;
104-
}
105-
106-
for (i = 0; i < rank; i++)
107-
{
108-
for (j = 0; j <= i; j++)
109-
{
110-
e = TEMPLATE(T, mat_entry) (B, j, pivots[i]);
111-
TEMPLATE(T, mat_entry_set) (U, j, i, e, ctx);
112-
}
113-
}
114-
115-
for (i = 0; i < n - rank; i++)
116-
{
117-
for (j = 0; j < rank; j++)
118-
{
119-
e = TEMPLATE(T, mat_entry) (B, j, nonpivots[i]);
120-
TEMPLATE(T, mat_entry_set) (V, j, i, e, ctx);
121-
}
122-
}
123-
124-
TEMPLATE(T, mat_solve_triu) (V, U, V, 0, ctx);
125-
126-
/* Clear pivot columns */
127-
for (i = 0; i < rank; i++)
128-
{
129-
for (j = 0; j <= i; j++)
130-
{
131-
if (i == j)
132-
{
133-
TEMPLATE(T, one) (TEMPLATE(T, mat_entry) (B, j, pivots[i]),
134-
ctx);
135-
}
136-
else
137-
{
138-
TEMPLATE(T, zero) (TEMPLATE(T, mat_entry) (B, j, pivots[i]),
139-
ctx);
140-
}
141-
}
142-
}
143-
144-
/* Write back the actual content */
145-
for (i = 0; i < n - rank; i++)
146-
{
147-
for (j = 0; j < rank; j++)
148-
TEMPLATE(T, mat_entry_set) (B, j, nonpivots[i],
149-
TEMPLATE(T, mat_entry) (V, j, i), ctx);
150-
}
151-
152-
TEMPLATE(T, mat_clear) (U, ctx);
153-
TEMPLATE(T, mat_clear) (V, ctx);
154-
155-
flint_free(pivots);
156-
flint_free(nonpivots);
157-
22+
gr_ctx_t gr_ctx;
23+
slong rank;
24+
TEMPLATE3(_gr_ctx_init, T, from_ref)(gr_ctx, ctx);
25+
GR_MUST_SUCCEED(gr_mat_rref_lu(&rank, (gr_mat_struct *) B, (const gr_mat_struct *) A, gr_ctx));
15826
return rank;
15927
}
16028

161-
16229
#endif

‎src/gr_mat.h

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -255,6 +255,8 @@ WARN_UNUSED_RESULT int gr_mat_rref(slong * res_rank, gr_mat_t R, const gr_mat_t
255255
WARN_UNUSED_RESULT int gr_mat_rref_den_fflu(slong * res_rank, gr_mat_t R, gr_ptr den, const gr_mat_t A, gr_ctx_t ctx);
256256
WARN_UNUSED_RESULT int gr_mat_rref_den(slong * res_rank, gr_mat_t R, gr_ptr den, const gr_mat_t A, gr_ctx_t ctx);
257257

258+
WARN_UNUSED_RESULT int gr_mat_nullspace_from_rref(gr_mat_t X, const gr_mat_t A, gr_srcptr Aden, slong rank, gr_ctx_t ctx);
259+
WARN_UNUSED_RESULT int gr_mat_nullspace_no_resize(slong * nullity, gr_mat_t X, const gr_mat_t A, gr_ctx_t ctx);
258260
WARN_UNUSED_RESULT int gr_mat_nullspace(gr_mat_t X, const gr_mat_t A, gr_ctx_t ctx);
259261

260262
WARN_UNUSED_RESULT int gr_mat_ones(gr_mat_t mat, gr_ctx_t ctx);

0 commit comments

Comments
 (0)