1 /*
2 Copyright (C) 2012 Fredrik Johansson
3 Copyright (C) 2013 Mike Hansen
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <https://www.gnu.org/licenses/>.
11 */
12
13 #ifdef T
14
15 #include "templates.h"
16
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <gmp.h>
20 #include "perm.h"
21 #include "flint.h"
22
23 int
check_rref_form(slong * perm,TEMPLATE (T,mat_t)A,slong rank,const TEMPLATE (T,ctx_t)ctx)24 check_rref_form(slong * perm, TEMPLATE(T, mat_t) A, slong rank,
25 const TEMPLATE(T, ctx_t) ctx)
26 {
27 slong i, j, k, prev_pivot;
28
29 /* bottom should be zero */
30 for (i = rank; i < A->r; i++)
31 for (j = 0; j < A->c; j++)
32 if (!TEMPLATE(T, is_zero) (TEMPLATE(T, mat_entry) (A, i, j), ctx))
33 return 0;
34
35 prev_pivot = -1;
36
37 for (i = 0; i < rank; i++)
38 {
39 for (j = 0; j < A->c; j++)
40 {
41 if (!TEMPLATE(T, is_zero) (TEMPLATE(T, mat_entry) (A, i, j), ctx))
42 {
43 /* pivot should have a higher column index than previous */
44 if (j <= prev_pivot)
45 return 0;
46
47 /* column should be 0 ... 0 1 0 ... 0 */
48 for (k = 0; k < rank; k++)
49 {
50 if (i == k)
51 {
52 if (!TEMPLATE(T, is_one)
53 (TEMPLATE(T, mat_entry) (A, k, j), ctx))
54 return 0;
55 }
56 else
57 {
58 if (!TEMPLATE(T, is_zero)
59 (TEMPLATE(T, mat_entry) (A, k, j), ctx))
60 return 0;
61 }
62 }
63
64 prev_pivot = j;
65 break;
66 }
67 }
68 }
69
70 return 1;
71 }
72
73 int
main(void)74 main(void)
75 {
76 slong i;
77 FLINT_TEST_INIT(state);
78
79 printf("rref....");
80 fflush(stdout);
81
82 for (i = 0; i < 10 * flint_test_multiplier(); i++)
83 {
84 TEMPLATE(T, ctx_t) ctx;
85 TEMPLATE(T, mat_t) A, B, C, D;
86 TEMPLATE(T, t) c;
87 slong j, k, m, n, rank1, rank2;
88 slong *perm;
89 int equal;
90
91 TEMPLATE(T, ctx_randtest) (ctx, state);
92
93 TEMPLATE(T, init) (c, ctx);
94
95 m = n_randint(state, 20);
96 n = n_randint(state, 20);
97 perm = _perm_init(2 * m);
98
99 TEMPLATE(T, mat_init) (A, m, n, ctx);
100 TEMPLATE(T, mat_init) (D, 2 * m, n, ctx);
101
102 TEMPLATE(T, mat_randtest) (A, state, ctx);
103 TEMPLATE(T, mat_init_set) (B, A, ctx);
104 TEMPLATE(T, mat_init_set) (C, A, ctx);
105
106 rank1 = TEMPLATE(T, mat_rref) (B, ctx);
107
108 if (!check_rref_form(perm, B, rank1, ctx))
109 {
110 printf("FAIL (malformed rref)\n");
111 TEMPLATE(T, mat_print_pretty) (A, ctx);
112 printf("\n\n");
113 TEMPLATE(T, mat_print_pretty) (B, ctx);
114 printf("\n\n");
115 abort();
116 }
117
118 /* Concatenate the original matrix with the rref, scramble the rows,
119 and check that the rref is the same */
120 _perm_randtest(perm, 2 * m, state);
121
122 for (j = 0; j < m; j++)
123 {
124 TEMPLATE(T, randtest_not_zero) (c, state, ctx);
125 for (k = 0; k < n; k++)
126 {
127 TEMPLATE(T, mul) (TEMPLATE(T, mat_entry) (D, perm[j], k),
128 TEMPLATE(T, mat_entry) (A, j, k), c, ctx);
129 }
130 }
131
132 for (j = 0; j < m; j++)
133 {
134 TEMPLATE(T, randtest_not_zero) (c, state, ctx);
135 for (k = 0; k < n; k++)
136 {
137 TEMPLATE(T, mul) (TEMPLATE(T, mat_entry) (D, perm[m + j], k),
138 TEMPLATE(T, mat_entry) (B, j, k), c, ctx);
139
140 }
141 }
142
143 rank2 = TEMPLATE(T, mat_rref) (D, ctx);
144 equal = (rank1 == rank2);
145
146 if (equal)
147 {
148 for (j = 0; j < rank2; j++)
149 for (k = 0; k < n; k++)
150 {
151 equal = equal
152 && TEMPLATE(T,
153 equal) (TEMPLATE(T, mat_entry) (B, j, k),
154 TEMPLATE(T, mat_entry) (D, j, k),
155 ctx);
156 }
157 for (j = rank2; j < 2 * rank2; j++)
158 for (k = 0; k < n; k++)
159 {
160 equal = equal
161 && TEMPLATE(T,
162 is_zero) (TEMPLATE(T, mat_entry) (D, j, k),
163 ctx);
164 }
165 }
166
167 if (!equal)
168 {
169 flint_printf("FAIL (rank1 = %wd, rank2 = %wd)!\n", rank1, rank2);
170 TEMPLATE(T, mat_print_pretty) (A, ctx);
171 printf("\n\n");
172 TEMPLATE(T, mat_print_pretty) (B, ctx);
173 printf("\n\n");
174 TEMPLATE(T, mat_print_pretty) (D, ctx);
175 printf("\n\n");
176 abort();
177 }
178
179 _perm_clear(perm);
180 TEMPLATE(T, mat_clear) (A, ctx);
181 TEMPLATE(T, mat_clear) (B, ctx);
182 TEMPLATE(T, mat_clear) (C, ctx);
183 TEMPLATE(T, mat_clear) (D, ctx);
184
185 TEMPLATE(T, clear) (c, ctx);
186 TEMPLATE(T, ctx_clear) (ctx);
187 }
188
189 FLINT_TEST_CLEANUP(state);
190 printf("PASS\n");
191 return 0;
192 }
193
194
195 #endif
196