1 /*******************************************************************
2 *
3 * M4RI: Linear Algebra over GF(2)
4 *
5 * Copyright (C) 2010 Martin Albrecht <M.R.Albrecht@rhul.ac.uk>
6 *
7 * Distributed under the terms of the GNU General Public License (GPL)
8 * version 2 or higher.
9 *
10 * This code is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * The full text of the GPL is available at:
16 *
17 * http://www.gnu.org/licenses/
18 *
19 ********************************************************************/
20
21 #ifdef HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24
25 #include "echelonform.h"
26 #include "brilliantrussian.h"
27 #include "ple.h"
28 #include "triangular.h"
29
mzd_echelonize(mzd_t * A,int full)30 rci_t mzd_echelonize(mzd_t *A, int full) {
31 return _mzd_echelonize_m4ri(A, full, 0, 1, __M4RI_ECHELONFORM_CROSSOVER_DENSITY);
32 }
33
mzd_echelonize_m4ri(mzd_t * A,int full,int k)34 rci_t mzd_echelonize_m4ri(mzd_t *A, int full, int k) {
35 return _mzd_echelonize_m4ri(A, full, k, 0, 1.0);
36 }
37
mzd_echelonize_pluq(mzd_t * A,int full)38 rci_t mzd_echelonize_pluq(mzd_t *A, int full) {
39 mzp_t *P = mzp_init(A->nrows);
40 mzp_t *Q = mzp_init(A->ncols);
41
42 rci_t r;
43 if(full) {
44 #if 0
45 r = mzd_pluq(A, P, Q, 0);
46 mzd_t *U = mzd_init_window(A, 0, 0, r, r);
47 mzd_t *B = mzd_init_window(A, 0, r, r, A->ncols);
48 if(r!=A->ncols)
49 mzd_trsm_upper_left(U, B, 0);
50 if(r)
51 mzd_set_ui(U, 0);
52 for(rci_t i = 0; i < r; ++i)
53 mzd_write_bit(A, i, i, 1);
54 mzd_free_window(U);
55 mzd_free_window(B);
56
57 if(r) {
58 mzd_t *A0 = mzd_init_window(A, 0, 0, r, A->ncols);
59 mzd_apply_p_right(A0, Q);
60 mzd_free_window(A0);
61 } else {
62 mzd_apply_p_right(A, Q);
63 }
64 #else
65 r = mzd_pluq(A, P, Q, 0);
66
67 mzd_t *U = mzd_init_window(A, 0, 0, r, r);
68 const rci_t r_radix = m4ri_radix*(r/m4ri_radix);
69
70 if(r_radix == r && r!=A->ncols) {
71
72 mzd_t *B = mzd_init_window(A, 0, r, r, A->ncols);
73 if(r!=A->ncols)
74 mzd_trsm_upper_left(U, B, 0);
75 mzd_free_window(B);
76
77 } else if (r_radix != r && r!=A->ncols) {
78 assert(r_radix < r);
79
80 if(A->ncols > r_radix+m4ri_radix) {
81 mzd_t *B0 = mzd_submatrix(NULL, A, 0, r_radix, r, r_radix+m4ri_radix);
82 mzd_t *B0w = mzd_init_window( A, 0, r_radix, r, r_radix+m4ri_radix);
83 mzd_t *B1 = mzd_init_window(A, 0, r_radix+m4ri_radix, r, A->ncols);
84
85 mzd_trsm_upper_left(U, B0, 0);
86 mzd_trsm_upper_left(U, B1, 0);
87
88 mzd_copy(B0w, B0);
89 mzd_free(B0);
90 mzd_free_window(B0w);
91 mzd_free_window(B1);
92
93 } else {
94 mzd_t *B = mzd_submatrix(NULL, A, 0, r_radix, r, A->ncols);
95 mzd_t *Bw = mzd_init_window(A, 0, r_radix, r, A->ncols);
96
97 mzd_trsm_upper_left(U, B, 0);
98
99 mzd_copy(Bw, B);
100 mzd_free_window(Bw);
101 mzd_free(B);
102 }
103 }
104
105 mzd_set_ui(U, 1);
106 mzd_free_window(U);
107
108 if(r) {
109 mzd_t *A0 = mzd_init_window(A, 0, 0, r, A->ncols);
110 mzd_apply_p_right(A0, Q);
111 mzd_free_window(A0);
112 }
113 #endif
114
115 } else {
116 r = mzd_ple(A, P, Q, 0);
117 for(rci_t i = 0; i < r; ++i) {
118 for(rci_t j = 0; j <= i; j += m4ri_radix) {
119 int const length = MIN(m4ri_radix, i - j + 1);
120 mzd_clear_bits(A, i, j, length);
121 }
122 mzd_write_bit(A, i, Q->values[i], 1);
123 }
124 }
125
126 if(r != A->nrows) {
127 mzd_t *R = mzd_init_window(A, r, 0, A->nrows, A->ncols);
128 mzd_set_ui(R, 0);
129 mzd_free_window(R);
130 }
131
132 mzp_free(P);
133 mzp_free(Q);
134
135 __M4RI_DD_MZD(A);
136 __M4RI_DD_RCI(r);
137 return r;
138 }
139