1 /*
2 Copyright (C) 2010 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 <https://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 "ulong_extras.h"
20
21
22 int
main(void)23 main(void)
24 {
25 fmpz_mat_t A, B, C, I;
26 fmpz_t den;
27 slong i, j, m, r;
28
29 FLINT_TEST_INIT(state);
30
31 flint_printf("inv....");
32 fflush(stdout);
33
34 {
35 fmpz_t d;
36 fmpz_mat_t A, B, C;
37
38 fmpz_mat_init(A, 1, 1);
39 fmpz_one(fmpz_mat_entry(A, 0, 0));
40
41 fmpz_mat_window_init(B, A, 0, 0, 1, 1);
42
43 fmpz_mat_init(C, 1, 1);
44 fmpz_init(d);
45
46 fmpz_mat_inv(C, d, B);
47
48 fmpz_clear(d);
49 fmpz_mat_clear(C);
50 fmpz_mat_window_clear(B);
51 fmpz_mat_clear(A);
52 }
53
54 for (i = 0; i < 1000 * flint_test_multiplier(); i++)
55 {
56 m = n_randint(state, 10);
57
58 fmpz_mat_init(A, m, m);
59 fmpz_mat_init(B, m, m);
60 fmpz_mat_init(C, m, m);
61 fmpz_mat_init(I, m, m);
62 fmpz_init(den);
63
64 for (j = 0; j < m; j++)
65 fmpz_set_ui(&I->rows[j][j], UWORD(1));
66
67 /* Verify that A * A^-1 = I for random matrices */
68
69 fmpz_mat_randrank(A, state, m, 1+n_randint(state, 2)*n_randint(state, 100));
70 /* Dense or sparse? */
71 if (n_randint(state, 2))
72 fmpz_mat_randops(A, state, 1+n_randint(state, 1 + m*m));
73
74 fmpz_mat_inv(B, den, A);
75 fmpz_mat_mul(C, A, B);
76
77 _fmpz_vec_scalar_divexact_fmpz(C->entries, C->entries, m*m, den);
78
79 if (!fmpz_mat_equal(C, I))
80 {
81 flint_printf("FAIL:\n");
82 flint_printf("A * A^-1 != I!\n");
83 flint_printf("A:\n"), fmpz_mat_print_pretty(A), flint_printf("\n");
84 flint_printf("A^-1:\n"), fmpz_mat_print_pretty(B), flint_printf("\n");
85 flint_printf("den(A^-1) = "), fmpz_print(den), flint_printf("\n");
86 flint_printf("A * A^-1:\n"), fmpz_mat_print_pretty(C), flint_printf("\n");
87 abort();
88 }
89
90 /* Test aliasing */
91 fmpz_mat_set(C, A);
92 fmpz_mat_inv(A, den, A);
93 fmpz_mat_mul(B, A, C);
94 _fmpz_vec_scalar_divexact_fmpz(B->entries, B->entries, m*m, den);
95
96 if (!fmpz_mat_equal(B, I))
97 {
98 flint_printf("FAIL:\n");
99 flint_printf("aliasing failed!\n");
100 fmpz_mat_print(C); flint_printf("\n");
101 abort();
102 }
103
104 fmpz_mat_clear(A);
105 fmpz_mat_clear(B);
106 fmpz_mat_clear(C);
107 fmpz_mat_clear(I);
108 fmpz_clear(den);
109 }
110
111 /* Test singular matrices */
112 /* Test singular systems */
113 for (i = 0; i < 1000 * flint_test_multiplier(); i++)
114 {
115 m = 1 + n_randint(state, 10);
116 r = n_randint(state, m);
117
118 fmpz_mat_init(A, m, m);
119 fmpz_mat_init(B, m, m);
120 fmpz_init(den);
121
122 fmpz_mat_randrank(A, state, r, 1+n_randint(state, 2)*n_randint(state, 100));
123
124 /* Dense */
125 if (n_randint(state, 2))
126 fmpz_mat_randops(A, state, 1+n_randint(state, 1 + m*m));
127
128 fmpz_mat_inv(B, den, A);
129 if (!fmpz_is_zero(den))
130 {
131 flint_printf("FAIL:\n");
132 flint_printf("singular system gave nonzero denominator\n");
133 abort();
134 }
135
136 /* Aliasing */
137 fmpz_mat_inv(A, den, A);
138 if (!fmpz_is_zero(den))
139 {
140 flint_printf("FAIL:\n");
141 flint_printf("singular system gave nonzero denominator\n");
142 abort();
143 }
144
145 fmpz_mat_clear(A);
146 fmpz_mat_clear(B);
147 fmpz_clear(den);
148 }
149
150 FLINT_TEST_CLEANUP(state);
151
152 flint_printf("PASS\n");
153 return 0;
154 }
155