1 /*
2 Copyright (C) 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 "flint.h"
15 #include "ulong_extras.h"
16 #include "nmod_vec.h"
17 #include "nmod_mat.h"
18
19
20 static __inline__ int
nmod_mat_pivot(nmod_mat_t A,slong * P,slong start_row,slong col)21 nmod_mat_pivot(nmod_mat_t A, slong * P, slong start_row, slong col)
22 {
23 slong j, t;
24 mp_ptr u;
25
26 if (nmod_mat_entry(A, start_row, col) != 0)
27 return 1;
28
29 for (j = start_row + 1; j < A->r; j++)
30 {
31 if (nmod_mat_entry(A, j, col) != 0)
32 {
33 u = A->rows[j];
34 A->rows[j] = A->rows[start_row];
35 A->rows[start_row] = u;
36
37 t = P[j];
38 P[j] = P[start_row];
39 P[start_row] = t;
40
41 return -1;
42 }
43 }
44 return 0;
45 }
46
47
48 slong
nmod_mat_lu_classical(slong * P,nmod_mat_t A,int rank_check)49 nmod_mat_lu_classical(slong * P, nmod_mat_t A, int rank_check)
50 {
51 mp_limb_t d, e, **a;
52 nmod_t mod;
53 slong i, m, n, rank, length, row, col;
54
55 m = A->r;
56 n = A->c;
57 a = A->rows;
58 mod = A->mod;
59
60 rank = row = col = 0;
61
62 for (i = 0; i < m; i++)
63 P[i] = i;
64
65 while (row < m && col < n)
66 {
67 if (nmod_mat_pivot(A, P, row, col) == 0)
68 {
69 if (rank_check)
70 return 0;
71 col++;
72 continue;
73 }
74
75 rank++;
76
77 d = a[row][col];
78 d = n_invmod(d, mod.n);
79 length = n - col - 1;
80
81 for (i = row + 1; i < m; i++)
82 {
83 e = n_mulmod2_preinv(a[i][col], d, mod.n, mod.ninv);
84 if (length != 0)
85 _nmod_vec_scalar_addmul_nmod(a[i] + col + 1,
86 a[row] + col + 1, length, nmod_neg(e, mod), mod);
87
88 a[i][col] = 0;
89 a[i][rank - 1] = e;
90 }
91 row++;
92 col++;
93 }
94
95 return rank;
96 }
97