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 "templates.h"
13
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <limits.h>
17
18 int
main(void)19 main(void)
20 {
21 TEMPLATE(T, mat_t) A, B, C, I;
22 TEMPLATE(T, ctx_t) ctx;
23 slong i, j, m, r;
24 int result;
25 FLINT_TEST_INIT(state);
26
27
28 flint_printf("inv....");
29 fflush(stdout);
30
31 for (i = 0; i < 5 * flint_test_multiplier(); i++)
32 {
33 m = n_randint(state, 20);
34
35 TEMPLATE(T, ctx_randtest) (ctx, state);
36
37 TEMPLATE(T, mat_init)(A, m, m, ctx);
38 TEMPLATE(T, mat_init)(B, m, m, ctx);
39 TEMPLATE(T, mat_init)(C, m, m, ctx);
40 TEMPLATE(T, mat_init)(I, m, m, ctx);
41
42 for (j = 0; j < m; j++)
43 TEMPLATE(T, one)(TEMPLATE(T, mat_entry)(I, j, j), ctx);
44
45 /* Verify that A * A^-1 = I for random matrices */
46
47 TEMPLATE(T, mat_randrank)(A, state, m, ctx);
48 /* Dense or sparse? */
49 if (n_randint(state, 2))
50 TEMPLATE(T, mat_randops)(A, 1+n_randint(state, 1+m*m), state, ctx);
51
52 result = TEMPLATE(T, mat_inv)(B, A, ctx);
53 TEMPLATE(T, mat_mul)(C, A, B, ctx);
54
55 if (!TEMPLATE(T, mat_equal)(C, I, ctx) || !result)
56 {
57 flint_printf("FAIL:\n");
58 flint_printf("A * A^-1 != I!\n");
59 flint_printf("A:\n");
60 TEMPLATE(T, mat_print_pretty)(A, ctx);
61 flint_printf("A^-1:\n");
62 TEMPLATE(T, mat_print_pretty)(B, ctx);
63 flint_printf("A * A^-1:\n");
64 TEMPLATE(T, mat_print_pretty)(C, ctx);
65 flint_printf("\n");
66 abort();
67 }
68
69 /* Test aliasing */
70 TEMPLATE(T, mat_set)(C, A, ctx);
71 TEMPLATE(T, mat_inv)(A, A, ctx);
72 TEMPLATE(T, mat_mul)(B, A, C, ctx);
73
74 if (!TEMPLATE(T, mat_equal)(B, I, ctx))
75 {
76 flint_printf("FAIL:\n");
77 flint_printf("aliasing failed!\n");
78 TEMPLATE(T, mat_print_pretty)(C, ctx);
79 abort();
80 }
81
82 TEMPLATE(T, mat_clear)(A, ctx);
83 TEMPLATE(T, mat_clear)(B, ctx);
84 TEMPLATE(T, mat_clear)(C, ctx);
85 TEMPLATE(T, mat_clear)(I, ctx);
86
87 TEMPLATE(T, ctx_clear) (ctx);
88 }
89
90 /* Test singular systems */
91 for (i = 0; i < 1000 * flint_test_multiplier(); i++)
92 {
93 m = 1 + n_randint(state, 20);
94 r = n_randint(state, m);
95 TEMPLATE(T, ctx_randtest) (ctx, state);
96
97 TEMPLATE(T, mat_init)(A, m, m, ctx);
98 TEMPLATE(T, mat_init)(B, m, m, ctx);
99
100 TEMPLATE(T, mat_randrank)(A, state, r, ctx);
101
102 /* Dense */
103 if (n_randint(state, 2))
104 TEMPLATE(T, mat_randops)(A, 1+n_randint(state, 1+m*m), state, ctx);
105
106 result = TEMPLATE(T, mat_inv)(B, A, ctx);
107
108 if (result)
109 {
110 flint_printf("FAIL:\n");
111 flint_printf("singular matrix reported as invertible\n");
112 abort();
113 }
114
115 /* Aliasing */
116 result = TEMPLATE(T, mat_inv)(A, A, ctx);
117 if (result)
118 {
119 flint_printf("FAIL:\n");
120 flint_printf("singular matrix reported as invertible\n");
121 abort();
122 }
123
124 TEMPLATE(T, mat_clear)(A, ctx);
125 TEMPLATE(T, mat_clear)(B, ctx);
126
127 TEMPLATE(T, ctx_clear) (ctx);
128 }
129
130 FLINT_TEST_CLEANUP(state);
131
132 flint_printf("PASS\n");
133 return 0;
134 }
135