1 /*
2 Copyright (C) 2014 Abhinav Baid
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 "mpf_mat.h"
13 #include "ulong_extras.h"
14
15 int
main(void)16 main(void)
17 {
18 int i;
19 FLINT_TEST_INIT(state);
20
21 flint_printf("qr....");
22 fflush(stdout);
23
24 for (i = 0; i < 100 * flint_test_multiplier(); i++)
25 {
26 mpf_t dot, tmp;
27 int j, k, l;
28 mpf_mat_t A, Q, R, B;
29 flint_bitcnt_t prec;
30
31 slong m, n;
32
33 m = n_randint(state, 10);
34 n = n_randint(state, 10);
35 prec = n_randint(state, 200) + 3;
36
37 mpf_mat_init(A, m, n, prec);
38 mpf_mat_init(Q, m, n, prec);
39 mpf_mat_init(R, n, n, prec);
40 mpf_mat_init(B, m, n, prec);
41 mpf_init2(dot, prec);
42 mpf_init2(tmp, prec);
43
44 mpf_mat_randtest(A, state, prec);
45
46 mpf_mat_qr(Q, R, A);
47
48 mpf_mat_mul(B, Q, R);
49
50 flint_mpf_set_ui(tmp, 1);
51 mpf_div_2exp(tmp, tmp, prec - 2);
52 for (j = 0; j < m; j++)
53 {
54 for (k = 0; k < n; k++)
55 {
56 mpf_sub(dot, mpf_mat_entry(A, j, k), mpf_mat_entry(B, j, k));
57 mpf_abs(dot, dot);
58 if (mpf_cmp(dot, tmp) > 0)
59 {
60 flint_printf("FAIL:\n");
61 flint_printf("A:\n");
62 mpf_mat_print(A);
63 flint_printf("Q:\n");
64 mpf_mat_print(Q);
65 flint_printf("R:\n");
66 mpf_mat_print(R);
67 flint_printf("B:\n");
68 mpf_mat_print(B);
69 mpf_out_str(stdout, 10, 0, dot);
70 flint_printf("\n");
71 flint_printf("%d\n", prec);
72 abort();
73 }
74 }
75 }
76
77 for (j = 0; j < n; j++)
78 {
79 mpf_t norm;
80 mpf_init2(norm, prec);
81 for (l = 0; l < m; l++)
82 {
83 mpf_mul(tmp, mpf_mat_entry(Q, l, j), mpf_mat_entry(Q, l, j));
84 mpf_add(norm, norm, tmp);
85 }
86
87 mpf_sub_ui(dot, norm, 1);
88 mpf_abs(dot, dot);
89 flint_mpf_set_ui(tmp, 1);
90 mpf_div_2exp(tmp, tmp, prec - 3);
91 if (flint_mpf_cmp_ui(norm, 0) != 0 && mpf_cmp(dot, tmp) > 0)
92 {
93 flint_printf("FAIL:\n");
94 flint_printf("Q:\n");
95 mpf_mat_print(Q);
96 mpf_out_str(stdout, 10, 0, norm);
97 flint_printf("\n");
98 flint_printf("%d\n", j);
99 abort();
100 }
101 mpf_clear(norm);
102 for (k = j + 1; k < n; k++)
103 {
104 flint_mpf_set_ui(dot, 0);
105 for (l = 0; l < m; l++)
106 {
107 mpf_mul(tmp, mpf_mat_entry(Q, l, j),
108 mpf_mat_entry(Q, l, k));
109 mpf_add(dot, dot, tmp);
110 }
111
112 flint_mpf_set_ui(tmp, 1);
113 mpf_div_2exp(tmp, tmp, prec);
114 mpf_abs(dot, dot);
115 if (mpf_cmp(dot, tmp) > 0)
116 {
117 flint_printf("FAIL:\n");
118 flint_printf("Q:\n");
119 mpf_mat_print(Q);
120 mpf_out_str(stdout, 10, 0, dot);
121 flint_printf("\n");
122 flint_printf("%d %d\n", j, k);
123 abort();
124 }
125 }
126 }
127
128 mpf_mat_clear(A);
129 mpf_mat_clear(Q);
130 mpf_mat_clear(R);
131 mpf_mat_clear(B);
132 mpf_clears(dot, tmp, NULL);
133 }
134
135 FLINT_TEST_CLEANUP(state);
136
137 flint_printf("PASS\n");
138 return EXIT_SUCCESS;
139 }
140