1 /*
2 Copyright (C) 2011 Fredrik Johansson
3 Copyright (C) 2013 Mike Hansen
4
5 This file is part of FLINT.
6
7 FLINT is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #ifdef T
14
15 #include "templates.h"
16
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include "flint.h"
20 #include "ulong_extras.h"
21
22 static __inline__ int
TEMPLATE(T,mat_pivot)23 TEMPLATE(T, mat_pivot) (TEMPLATE(T, mat_t) A, slong * P, slong start_row,
24 slong col, const TEMPLATE(T, ctx_t) ctx)
25 {
26 slong j, t;
27 TEMPLATE(T, struct) * u;
28
29 if (!TEMPLATE(T, is_zero)
30 (TEMPLATE(T, mat_entry) (A, start_row, col), ctx))
31 return 1;
32
33 for (j = start_row + 1; j < A->r; j++)
34 {
35 if (!TEMPLATE(T, is_zero) (TEMPLATE(T, mat_entry) (A, j, col), ctx))
36 {
37 u = A->rows[j];
38 A->rows[j] = A->rows[start_row];
39 A->rows[start_row] = u;
40
41 t = P[j];
42 P[j] = P[start_row];
43 P[start_row] = t;
44
45 return -1;
46 }
47 }
48 return 0;
49 }
50
51
52 slong
TEMPLATE(T,mat_lu_classical)53 TEMPLATE(T, mat_lu_classical) (slong * P,
54 TEMPLATE(T, mat_t) A,
55 int rank_check, const TEMPLATE(T, ctx_t) ctx)
56 {
57 TEMPLATE(T, t) d, e, neg_e;
58 TEMPLATE(T, struct) ** a;
59 slong i, m, n, rank, length, row, col;
60
61 m = A->r;
62 n = A->c;
63 a = A->rows;
64
65 rank = row = col = 0;
66
67 for (i = 0; i < m; i++)
68 P[i] = i;
69
70 TEMPLATE(T, init) (d, ctx);
71 TEMPLATE(T, init) (e, ctx);
72 TEMPLATE(T, init) (neg_e, ctx);
73
74 while (row < m && col < n)
75 {
76 if (TEMPLATE(T, mat_pivot) (A, P, row, col, ctx) == 0)
77 {
78 if (rank_check)
79 {
80 rank = 0;
81 goto cleanup;
82 }
83 col++;
84 continue;
85 }
86
87 rank++;
88
89 TEMPLATE(T, inv) (d, a[row] + col, ctx);
90
91 length = n - col - 1;
92
93 for (i = row + 1; i < m; i++)
94 {
95 TEMPLATE(T, mul) (e, a[i] + col, d, ctx);
96 if (length != 0)
97 {
98 TEMPLATE(T, neg) (neg_e, e, ctx);
99 _TEMPLATE3(T, vec_scalar_addmul, T) (a[i] + col + 1,
100 a[row] + col + 1,
101 length, neg_e, ctx);
102 }
103
104 TEMPLATE(T, zero) (a[i] + col, ctx);
105 TEMPLATE(T, set) (a[i] + rank - 1, e, ctx);
106 }
107 row++;
108 col++;
109 }
110
111 cleanup:
112 TEMPLATE(T, clear) (d, ctx);
113 TEMPLATE(T, clear) (e, ctx);
114 TEMPLATE(T, clear) (neg_e, ctx);
115
116 return rank;
117 }
118
119
120 #endif
121