1 /*
2 Copyright (C) 2010-2012 Fredrik Johansson
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <gmp.h>
15 #include "flint.h"
16 #include "fmpz.h"
17 #include "fmpz_vec.h"
18 #include "fmpz_mat.h"
19 #include "perm.h"
20 #include "ulong_extras.h"
21
22 int
main(void)23 main(void)
24 {
25 slong iter;
26 FLINT_TEST_INIT(state);
27
28 flint_printf("rref_fflu....");
29 fflush(stdout);
30
31
32
33 for (iter = 0; iter < 10000 * flint_test_multiplier(); iter++)
34 {
35 fmpz_mat_t A, R, B, R2;
36 fmpz_t den, c, den2;
37 slong j, k, m, n, b, d, r, rank1, rank2;
38 slong *perm;
39 int equal;
40
41 m = n_randint(state, 10);
42 n = n_randint(state, 10);
43 r = n_randint(state, FLINT_MIN(m, n) + 1);
44
45 fmpz_mat_init(A, m, n);
46 fmpz_mat_init(R, m, n);
47 fmpz_mat_init(B, 2 * m, n);
48 fmpz_mat_init(R2, 2 * m, n);
49
50 fmpz_init(c);
51 fmpz_init(den);
52 fmpz_init(den2);
53
54 perm = _perm_init(2 * m);
55
56 /* sparse */
57 b = 1 + n_randint(state, 10) * n_randint(state, 10);
58 d = n_randint(state, 2*m*n + 1);
59 fmpz_mat_randrank(A, state, r, b);
60
61 /* dense */
62 if (n_randint(state, 2))
63 fmpz_mat_randops(A, state, d);
64
65 rank1 = fmpz_mat_rref(R, den, A);
66
67 if (r != rank1)
68 {
69 flint_printf("FAIL:\n");
70 flint_printf("wrong rank!\n");
71 abort();
72 }
73
74 if (!fmpz_mat_is_in_rref_with_rank(R, den, rank1))
75 {
76 flint_printf("FAIL matrix not in rref!\n");
77 fmpz_mat_print_pretty(A); flint_printf("\n\n");
78 fmpz_mat_print_pretty(R); flint_printf("\n\n");
79 abort();
80 }
81
82 /* Concatenate the original matrix with the rref, scramble the rows,
83 and check that the rref is the same */
84 _perm_randtest(perm, 2 * m, state);
85
86 for (j = 0; j < m; j++)
87 {
88 fmpz_randtest_not_zero(c, state, 5);
89 for (k = 0; k < n; k++)
90 fmpz_mul(fmpz_mat_entry(B, perm[j], k), fmpz_mat_entry(A, j, k), c);
91 }
92
93 for (j = 0; j < m; j++)
94 {
95 fmpz_randtest_not_zero(c, state, 5);
96 for (k = 0; k < n; k++)
97 fmpz_mul(fmpz_mat_entry(B, perm[m + j], k), fmpz_mat_entry(R, j, k), c);
98 }
99
100 rank2 = fmpz_mat_rref(R2, den2, B);
101 equal = (rank1 == rank2);
102
103 if (equal)
104 {
105 fmpz_mat_scalar_mul_fmpz(R, R, den2);
106 fmpz_mat_scalar_mul_fmpz(R2, R2, den);
107
108 for (j = 0; j < rank2; j++)
109 for (k = 0; k < n; k++)
110 equal = equal &&
111 fmpz_equal(fmpz_mat_entry(R, j, k), fmpz_mat_entry(R2, j, k));
112 for (j = rank2; j < 2 * rank2; j++)
113 for (k = 0; k < n; k++)
114 equal = equal && fmpz_is_zero(fmpz_mat_entry(R2, j, k));
115 }
116
117 if (!equal)
118 {
119 flint_printf("FAIL (rank1 = %wd, rank2 = %wd)!\n", rank1, rank2);
120 fmpz_mat_print_pretty(A); flint_printf("\n\n");
121 fmpz_mat_print_pretty(R); flint_printf("\n\n");
122 fmpz_mat_print_pretty(R2); flint_printf("\n\n");
123 abort();
124 }
125
126 fmpz_clear(c);
127 fmpz_clear(den);
128 fmpz_clear(den2);
129
130 _perm_clear(perm);
131
132 fmpz_mat_clear(A);
133 fmpz_mat_clear(R);
134 fmpz_mat_clear(B);
135 fmpz_mat_clear(R2);
136 }
137
138 FLINT_TEST_CLEANUP(state);
139
140 flint_printf("PASS\n");
141 return 0;
142 }
143
144