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 <https://www.gnu.org/licenses/>.
10 */
11 
12 #include "mpf_mat.h"
13 
14 void
mpf_mat_gso(mpf_mat_t B,const mpf_mat_t A)15 mpf_mat_gso(mpf_mat_t B, const mpf_mat_t A)
16 {
17     slong i, j, k;
18     int flag;
19     mpf_t t, s, tmp, eps;
20     flint_bitcnt_t exp;
21 
22     if (B->r != A->r || B->c != A->c)
23     {
24         flint_printf("Exception (mpf_mat_gso). Incompatible dimensions.\n");
25         flint_abort();
26     }
27 
28     if (B == A)
29     {
30         mpf_mat_t T;
31         mpf_mat_init(T, A->r, A->c, B->prec);
32         mpf_mat_gso(T, A);
33         mpf_mat_swap_entrywise(B, T);
34         mpf_mat_clear(T);
35         return;
36     }
37 
38     if (A->r == 0)
39     {
40         return;
41     }
42 
43     mpf_init2(t, B->prec);
44     mpf_init2(s, B->prec);
45     mpf_init2(tmp, B->prec);
46     mpf_init2(eps, B->prec);
47     exp = ceil((A->prec) / 64.0) * 64;
48     flint_mpf_set_ui(eps, 1);
49     mpf_div_2exp(eps, eps, exp);
50 
51     for (k = 0; k < A->c; k++)
52     {
53         for (j = 0; j < A->r; j++)
54         {
55             mpf_set(mpf_mat_entry(B, j, k), mpf_mat_entry(A, j, k));
56         }
57         flag = 1;
58         while (flag)
59         {
60             flint_mpf_set_ui(t, 0);
61             for (i = 0; i < k; i++)
62             {
63                 flint_mpf_set_ui(s, 0);
64                 for (j = 0; j < A->r; j++)
65                 {
66                     mpf_mul(tmp, mpf_mat_entry(B, j, i),
67                             mpf_mat_entry(B, j, k));
68                     mpf_add(s, s, tmp);
69                 }
70                 mpf_mul(tmp, s, s);
71                 mpf_add(t, t, tmp);
72                 for (j = 0; j < A->r; j++)
73                 {
74                     mpf_mul(tmp, s, mpf_mat_entry(B, j, i));
75                     mpf_sub(mpf_mat_entry(B, j, k), mpf_mat_entry(B, j, k),
76                             tmp);
77                 }
78             }
79             flint_mpf_set_ui(s, 0);
80             for (j = 0; j < A->r; j++)
81             {
82                 mpf_mul(tmp, mpf_mat_entry(B, j, k), mpf_mat_entry(B, j, k));
83                 mpf_add(s, s, tmp);
84             }
85             mpf_add(t, t, s);
86             flag = 0;
87             if (mpf_cmp(s, t) < 0)
88             {
89                 if (mpf_cmp(s, eps) < 0)
90                     flint_mpf_set_ui(s, 0);
91                 else
92                     flag = 1;
93             }
94         }
95         mpf_sqrt(s, s);
96         if (flint_mpf_cmp_ui(s, 0) != 0)
97             mpf_ui_div(s, 1, s);
98         for (j = 0; j < A->r; j++)
99         {
100             mpf_mul(mpf_mat_entry(B, j, k), mpf_mat_entry(B, j, k), s);
101         }
102     }
103     mpf_clears(t, s, tmp, eps, NULL);
104 }
105