1 /*
2 Copyright (C) 2010,2011 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 <http://www.gnu.org/licenses/>.
10 */
11
12 #include <stdio.h>
13 #include <stdlib.h>
14 #include <gmp.h>
15 #include "flint.h"
16 #include "nmod_mat.h"
17 #include "nmod_vec.h"
18
19
20 void
nmod_mat_solve_tril_recursive(nmod_mat_t X,const nmod_mat_t L,const nmod_mat_t B,int unit)21 nmod_mat_solve_tril_recursive(nmod_mat_t X,
22 const nmod_mat_t L, const nmod_mat_t B,
23 int unit)
24 {
25 nmod_mat_t LA, LC, LD, XX, XY, BX, BY;
26 slong r, n, m;
27
28 n = L->r;
29 m = B->c;
30 r = n / 2;
31
32 if (n == 0 || m == 0)
33 return;
34
35 /*
36 Denoting inv(M) by M^, we have:
37
38 [A 0]^ [X] == [A^ 0 ] [X] == [A^ X]
39 [C D] [Y] == [-D^ C A^ D^] [Y] == [D^ (Y - C A^ X)]
40 */
41
42 nmod_mat_window_init(LA, L, 0, 0, r, r);
43 nmod_mat_window_init(LC, L, r, 0, n, r);
44 nmod_mat_window_init(LD, L, r, r, n, n);
45 nmod_mat_window_init(BX, B, 0, 0, r, m);
46 nmod_mat_window_init(BY, B, r, 0, n, m);
47 nmod_mat_window_init(XX, X, 0, 0, r, m);
48 nmod_mat_window_init(XY, X, r, 0, n, m);
49
50 nmod_mat_solve_tril(XX, LA, BX, unit);
51 nmod_mat_submul(XY, BY, LC, XX);
52 nmod_mat_solve_tril(XY, LD, XY, unit);
53
54 nmod_mat_window_clear(LA);
55 nmod_mat_window_clear(LC);
56 nmod_mat_window_clear(LD);
57 nmod_mat_window_clear(BX);
58 nmod_mat_window_clear(BY);
59 nmod_mat_window_clear(XX);
60 nmod_mat_window_clear(XY);
61 }
62