diff --git a/doc/source/fq_mat.rst b/doc/source/fq_mat.rst index 3f0bbec2b4..14e292f3b2 100644 --- a/doc/source/fq_mat.rst +++ b/doc/source/fq_mat.rst @@ -355,23 +355,6 @@ LU decomposition will abandon the output matrix in an undefined state and return 0 if `A` is detected to be rank-deficient. - This function calls ``fq_mat_lu_recursive``. - -.. function:: slong fq_mat_lu_classical(slong * P, fq_mat_t A, int rank_check, const fq_ctx_t ctx) - - Computes a generalised LU decomposition `LU = PA` of a given - matrix `A`, returning the rank of `A`. The behavior of this - function is identical to that of ``fq_mat_lu``. Uses Gaussian - elimination. - -.. function:: slong fq_mat_lu_recursive(slong * P, fq_mat_t A, int rank_check, const fq_ctx_t ctx) - - Computes a generalised LU decomposition `LU = PA` of a given - matrix `A`, returning the rank of `A`. The behavior of this - function is identical to that of ``fq_mat_lu``. Uses recursive - block decomposition, switching to classical Gaussian elimination - for sufficiently small blocks. - Reduced row echelon form -------------------------------------------------------------------------------- @@ -413,33 +396,6 @@ Triangular solving is allowed. Automatically chooses between the classical and recursive algorithms. -.. function:: void fq_mat_solve_tril_classical(fq_mat_t X, const fq_mat_t L, const fq_mat_t B, int unit, const fq_ctx_t ctx) - - Sets `X = L^{-1} B` where `L` is a full rank lower triangular - square matrix. If ``unit`` = 1, `L` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. Uses forward substitution. - -.. function:: void fq_mat_solve_tril_recursive(fq_mat_t X, const fq_mat_t L, const fq_mat_t B, int unit, const fq_ctx_t ctx) - - Sets `X = L^{-1} B` where `L` is a full rank lower triangular - square matrix. If ``unit`` = 1, `L` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. - - Uses the block inversion formula - - .. math:: - \begin{pmatrix} A & 0 \\ C & D \end{pmatrix}^{-1} - \begin{pmatrix} X \\ Y \end{pmatrix} = - \begin{pmatrix} A^{-1} X \\ D^{-1} ( Y - C A^{-1} X ) \end{pmatrix} - - - to reduce the problem to matrix multiplication and triangular - solving of smaller systems. - .. function:: void fq_mat_solve_triu(fq_mat_t X, const fq_mat_t U, const fq_mat_t B, int unit, const fq_ctx_t ctx) Sets `X = U^{-1} B` where `U` is a full rank upper triangular @@ -449,33 +405,6 @@ Triangular solving is allowed. Automatically chooses between the classical and recursive algorithms. -.. function:: void fq_mat_solve_triu_classical(fq_mat_t X, const fq_mat_t U, const fq_mat_t B, int unit, const fq_ctx_t ctx) - - Sets `X = U^{-1} B` where `U` is a full rank upper triangular - square matrix. If ``unit`` = 1, `U` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. Uses forward substitution. - -.. function:: void fq_mat_solve_triu_recursive(fq_mat_t X, const fq_mat_t U, const fq_mat_t B, int unit, const fq_ctx_t ctx) - - Sets `X = U^{-1} B` where `U` is a full rank upper triangular - square matrix. If ``unit`` = 1, `U` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. - - Uses the block inversion formula - - .. math:: - \begin{pmatrix} A & B \\ 0 & D \end{pmatrix}^{-1} - \begin{pmatrix} X \\ Y \end{pmatrix} = - \begin{pmatrix} A^{-1} (X - B D^{-1} Y) \\ D^{-1} Y \end{pmatrix} - - - to reduce the problem to matrix multiplication and triangular - solving of smaller systems. - Solving -------------------------------------------------------------------------------- diff --git a/doc/source/fq_nmod_mat.rst b/doc/source/fq_nmod_mat.rst index fbb9074be4..bb6c1273fe 100644 --- a/doc/source/fq_nmod_mat.rst +++ b/doc/source/fq_nmod_mat.rst @@ -354,23 +354,6 @@ LU decomposition will abandon the output matrix in an undefined state and return 0 if `A` is detected to be rank-deficient. - This function calls ``fq_nmod_mat_lu_recursive``. - -.. function:: slong fq_nmod_mat_lu_classical(slong * P, fq_nmod_mat_t A, int rank_check, const fq_nmod_ctx_t ctx) - - Computes a generalised LU decomposition `LU = PA` of a given - matrix `A`, returning the rank of `A`. The behavior of this - function is identical to that of ``fq_nmod_mat_lu``. Uses Gaussian - elimination. - -.. function:: slong fq_nmod_mat_lu_recursive(slong * P, fq_nmod_mat_t A, int rank_check, const fq_nmod_ctx_t ctx) - - Computes a generalised LU decomposition `LU = PA` of a given - matrix `A`, returning the rank of `A`. The behavior of this - function is identical to that of ``fq_nmod_mat_lu``. Uses recursive - block decomposition, switching to classical Gaussian elimination - for sufficiently small blocks. - Reduced row echelon form -------------------------------------------------------------------------------- @@ -412,33 +395,6 @@ Triangular solving is allowed. Automatically chooses between the classical and recursive algorithms. -.. function:: void fq_nmod_mat_solve_tril_classical(fq_nmod_mat_t X, const fq_nmod_mat_t L, const fq_nmod_mat_t B, int unit, const fq_nmod_ctx_t ctx) - - Sets `X = L^{-1} B` where `L` is a full rank lower triangular - square matrix. If ``unit`` = 1, `L` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. Uses forward substitution. - -.. function:: void fq_nmod_mat_solve_tril_recursive(fq_nmod_mat_t X, const fq_nmod_mat_t L, const fq_nmod_mat_t B, int unit, const fq_nmod_ctx_t ctx) - - Sets `X = L^{-1} B` where `L` is a full rank lower triangular - square matrix. If ``unit`` = 1, `L` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. - - Uses the block inversion formula - - .. math:: - \begin{pmatrix} A & 0 \\ C & D \end{pmatrix}^{-1} - \begin{pmatrix} X \\ Y \end{pmatrix} = - \begin{pmatrix} A^{-1} X \\ D^{-1} ( Y - C A^{-1} X ) \end{pmatrix} - - - to reduce the problem to matrix multiplication and triangular - solving of smaller systems. - .. function:: void fq_nmod_mat_solve_triu(fq_nmod_mat_t X, const fq_nmod_mat_t U, const fq_nmod_mat_t B, int unit, const fq_nmod_ctx_t ctx) Sets `X = U^{-1} B` where `U` is a full rank upper triangular @@ -448,33 +404,6 @@ Triangular solving is allowed. Automatically chooses between the classical and recursive algorithms. -.. function:: void fq_nmod_mat_solve_triu_classical(fq_nmod_mat_t X, const fq_nmod_mat_t U, const fq_nmod_mat_t B, int unit, const fq_nmod_ctx_t ctx) - - Sets `X = U^{-1} B` where `U` is a full rank upper triangular - square matrix. If ``unit`` = 1, `U` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. Uses forward substitution. - -.. function:: void fq_nmod_mat_solve_triu_recursive(fq_nmod_mat_t X, const fq_nmod_mat_t U, const fq_nmod_mat_t B, int unit, const fq_nmod_ctx_t ctx) - - Sets `X = U^{-1} B` where `U` is a full rank upper triangular - square matrix. If ``unit`` = 1, `U` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. - - Uses the block inversion formula - - .. math:: - \begin{pmatrix} A & B \\ 0 & D \end{pmatrix}^{-1} - \begin{pmatrix} X \\ Y \end{pmatrix} = - \begin{pmatrix} A^{-1} (X - B D^{-1} Y) \\ D^{-1} Y \end{pmatrix} - - - to reduce the problem to matrix multiplication and triangular - solving of smaller systems. - Solving -------------------------------------------------------------------------------- diff --git a/doc/source/fq_zech_mat.rst b/doc/source/fq_zech_mat.rst index 82b16ea9a4..afb6820b67 100644 --- a/doc/source/fq_zech_mat.rst +++ b/doc/source/fq_zech_mat.rst @@ -321,23 +321,6 @@ LU decomposition will abandon the output matrix in an undefined state and return 0 if `A` is detected to be rank-deficient. - This function calls ``fq_zech_mat_lu_recursive``. - -.. function:: slong fq_zech_mat_lu_classical(slong * P, fq_zech_mat_t A, int rank_check, const fq_zech_ctx_t ctx) - - Computes a generalised LU decomposition `LU = PA` of a given - matrix `A`, returning the rank of `A`. The behavior of this - function is identical to that of ``fq_zech_mat_lu``. Uses Gaussian - elimination. - -.. function:: slong fq_zech_mat_lu_recursive(slong * P, fq_zech_mat_t A, int rank_check, const fq_zech_ctx_t ctx) - - Computes a generalised LU decomposition `LU = PA` of a given - matrix `A`, returning the rank of `A`. The behavior of this - function is identical to that of ``fq_zech_mat_lu``. Uses recursive - block decomposition, switching to classical Gaussian elimination - for sufficiently small blocks. - Reduced row echelon form -------------------------------------------------------------------------------- @@ -365,33 +348,6 @@ Triangular solving is allowed. Automatically chooses between the classical and recursive algorithms. -.. function:: void fq_zech_mat_solve_tril_classical(fq_zech_mat_t X, const fq_zech_mat_t L, const fq_zech_mat_t B, int unit, const fq_zech_ctx_t ctx) - - Sets `X = L^{-1} B` where `L` is a full rank lower triangular - square matrix. If ``unit`` = 1, `L` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. Uses forward substitution. - -.. function:: void fq_zech_mat_solve_tril_recursive(fq_zech_mat_t X, const fq_zech_mat_t L, const fq_zech_mat_t B, int unit, const fq_zech_ctx_t ctx) - - Sets `X = L^{-1} B` where `L` is a full rank lower triangular - square matrix. If ``unit`` = 1, `L` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. - - Uses the block inversion formula - - .. math:: - \begin{pmatrix} A & 0 \\ C & D \end{pmatrix}^{-1} - \begin{pmatrix} X \\ Y \end{pmatrix} = - \begin{pmatrix} A^{-1} X \\ D^{-1} ( Y - C A^{-1} X ) \end{pmatrix} - - - to reduce the problem to matrix multiplication and triangular - solving of smaller systems. - .. function:: void fq_zech_mat_solve_triu(fq_zech_mat_t X, const fq_zech_mat_t U, const fq_zech_mat_t B, int unit, const fq_zech_ctx_t ctx) Sets `X = U^{-1} B` where `U` is a full rank upper triangular @@ -401,32 +357,6 @@ Triangular solving is allowed. Automatically chooses between the classical and recursive algorithms. -.. function:: void fq_zech_mat_solve_triu_classical(fq_zech_mat_t X, const fq_zech_mat_t U, const fq_zech_mat_t B, int unit, const fq_zech_ctx_t ctx) - - Sets `X = U^{-1} B` where `U` is a full rank upper triangular - square matrix. If ``unit`` = 1, `U` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. Uses forward substitution. - -.. function:: void fq_zech_mat_solve_triu_recursive(fq_zech_mat_t X, const fq_zech_mat_t U, const fq_zech_mat_t B, int unit, const fq_zech_ctx_t ctx) - - Sets `X = U^{-1} B` where `U` is a full rank upper triangular - square matrix. If ``unit`` = 1, `U` is assumed to have ones on - its main diagonal, and the main diagonal will not be read. `X` - and `B` are allowed to be the same matrix, but no other aliasing - is allowed. - - Uses the block inversion formula - - .. math:: - \begin{pmatrix} A & B \\ 0 & D \end{pmatrix}^{-1} - \begin{pmatrix} X \\ Y \end{pmatrix} = - \begin{pmatrix} A^{-1} (X - B D^{-1} Y) \\ D^{-1} Y \end{pmatrix} - - - to reduce the problem to matrix multiplication and triangular - solving of smaller systems. Solving -------------------------------------------------------------------------------- diff --git a/src/fq/mat_templates.c b/src/fq/mat_templates.c index f1af4ebd6a..54521071e2 100644 --- a/src/fq/mat_templates.c +++ b/src/fq/mat_templates.c @@ -42,8 +42,6 @@ #include "fq_mat_templates/is_one.c" #include "fq_mat_templates/is_zero.c" #include "fq_mat_templates/lu.c" -#include "fq_mat_templates/lu_classical.c" -#include "fq_mat_templates/lu_recursive.c" #include "fq_mat_templates/mat_entry_set.c" #include "fq_mat_templates/mat_invert_cols.c" #include "fq_mat_templates/mat_swap_cols.c" @@ -70,11 +68,7 @@ #include "fq_mat_templates/similarity.c" #include "fq_mat_templates/solve.c" #include "fq_mat_templates/solve_tril.c" -#include "fq_mat_templates/solve_tril_classical.c" -#include "fq_mat_templates/solve_tril_recursive.c" #include "fq_mat_templates/solve_triu.c" -#include "fq_mat_templates/solve_triu_classical.c" -#include "fq_mat_templates/solve_triu_recursive.c" #include "fq_mat_templates/sub.c" #include "fq_mat_templates/submul.c" #include "fq_mat_templates/swap.c" diff --git a/src/fq_mat/test/main.c b/src/fq_mat/test/main.c index 9b3ba2d59b..a637320533 100644 --- a/src/fq_mat/test/main.c +++ b/src/fq_mat/test/main.c @@ -20,8 +20,7 @@ #include "t-inv.c" #include "t-invert_rows_cols.c" #include "t-is_zero.c" -#include "t-lu_classical.c" -#include "t-lu_recursive.c" +#include "t-lu.c" #include "t-minpoly.c" #include "t-mul.c" #include "t-mul_KS.c" @@ -34,11 +33,7 @@ #include "t-set_nmod_mat.c" #include "t-solve.c" #include "t-solve_tril.c" -#include "t-solve_tril_classical.c" -#include "t-solve_tril_recursive.c" #include "t-solve_triu.c" -#include "t-solve_triu_classical.c" -#include "t-solve_triu_recursive.c" #include "t-submul.c" #include "t-vec_mul.c" #include "t-window_init_clear.c" @@ -57,8 +52,7 @@ test_struct tests[] = TEST_FUNCTION(fq_mat_inv), TEST_FUNCTION(fq_mat_invert_rows_cols), TEST_FUNCTION(fq_mat_is_zero), - TEST_FUNCTION(fq_mat_lu_classical), - TEST_FUNCTION(fq_mat_lu_recursive), + TEST_FUNCTION(fq_mat_lu), TEST_FUNCTION(fq_mat_minpoly), TEST_FUNCTION(fq_mat_mul), TEST_FUNCTION(fq_mat_mul_KS), @@ -71,11 +65,7 @@ test_struct tests[] = TEST_FUNCTION(fq_mat_set_nmod_mat), TEST_FUNCTION(fq_mat_solve), TEST_FUNCTION(fq_mat_solve_tril), - TEST_FUNCTION(fq_mat_solve_tril_classical), - TEST_FUNCTION(fq_mat_solve_tril_recursive), TEST_FUNCTION(fq_mat_solve_triu), - TEST_FUNCTION(fq_mat_solve_triu_classical), - TEST_FUNCTION(fq_mat_solve_triu_recursive), TEST_FUNCTION(fq_mat_submul), TEST_FUNCTION(fq_mat_vec_mul), TEST_FUNCTION(fq_mat_window_init_clear), diff --git a/src/fq_mat/test/t-lu_classical.c b/src/fq_mat/test/t-lu.c similarity index 90% rename from src/fq_mat/test/t-lu_classical.c rename to src/fq_mat/test/t-lu.c index e60d80392f..5c88d8ff4f 100644 --- a/src/fq_mat/test/t-lu_classical.c +++ b/src/fq_mat/test/t-lu.c @@ -18,6 +18,6 @@ #define T fq #define CAP_T FQ -#include "fq_mat_templates/test/t-lu_classical.c" +#include "fq_mat_templates/test/t-lu.c" #undef CAP_T #undef T diff --git a/src/fq_mat/test/t-lu_recursive.c b/src/fq_mat/test/t-lu_recursive.c deleted file mode 100644 index b30a9b9e11..0000000000 --- a/src/fq_mat/test/t-lu_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq.h" -#include "fq_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq -#define CAP_T FQ -#include "fq_mat_templates/test/t-lu_recursive.c" -#undef CAP_T -#undef T diff --git a/src/fq_mat/test/t-solve_tril_classical.c b/src/fq_mat/test/t-solve_tril_classical.c deleted file mode 100644 index 38f6fcf4f4..0000000000 --- a/src/fq_mat/test/t-solve_tril_classical.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq.h" -#include "fq_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq -#define CAP_T FQ -#include "fq_mat_templates/test/t-solve_tril_classical.c" -#undef CAP_T -#undef T diff --git a/src/fq_mat/test/t-solve_tril_recursive.c b/src/fq_mat/test/t-solve_tril_recursive.c deleted file mode 100644 index 2d08d66f72..0000000000 --- a/src/fq_mat/test/t-solve_tril_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq.h" -#include "fq_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq -#define CAP_T FQ -#include "fq_mat_templates/test/t-solve_tril_recursive.c" -#undef CAP_T -#undef T diff --git a/src/fq_mat/test/t-solve_triu_classical.c b/src/fq_mat/test/t-solve_triu_classical.c deleted file mode 100644 index 0c86ab1d6b..0000000000 --- a/src/fq_mat/test/t-solve_triu_classical.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq.h" -#include "fq_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq -#define CAP_T FQ -#include "fq_mat_templates/test/t-solve_triu_classical.c" -#undef CAP_T -#undef T diff --git a/src/fq_mat/test/t-solve_triu_recursive.c b/src/fq_mat/test/t-solve_triu_recursive.c deleted file mode 100644 index 2a3b0ed911..0000000000 --- a/src/fq_mat/test/t-solve_triu_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq.h" -#include "fq_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq -#define CAP_T FQ -#include "fq_mat_templates/test/t-solve_triu_recursive.c" -#undef CAP_T -#undef T diff --git a/src/fq_mat_templates.h b/src/fq_mat_templates.h index c48123747f..6f320f1bb7 100644 --- a/src/fq_mat_templates.h +++ b/src/fq_mat_templates.h @@ -236,14 +236,6 @@ slong TEMPLATE(T, mat_lu)(slong * P, int rank_check, const TEMPLATE(T, ctx_t) ctx); -slong TEMPLATE(T, mat_lu_recursive)(slong * P, - TEMPLATE(T, mat_t) A, - int rank_check, - const TEMPLATE(T, ctx_t) ctx); - -slong TEMPLATE(T, mat_lu_classical)(slong * P, TEMPLATE(T, mat_t) A, int rank_check, - const TEMPLATE(T, ctx_t) ctx); - /* Inverse *******************************************************************/ int TEMPLATE(T, mat_inv)(TEMPLATE(T, mat_t) B, TEMPLATE(T, mat_t) A, @@ -267,33 +259,10 @@ void TEMPLATE(T, mat_solve_tril)(TEMPLATE(T, mat_t) X, const TEMPLATE(T, mat_t) const TEMPLATE(T, mat_t) B, int unit, const TEMPLATE(T, ctx_t) ctx); -void TEMPLATE(T, mat_solve_tril_classical)(TEMPLATE(T, mat_t) X, - const TEMPLATE(T, mat_t) L, - const TEMPLATE(T, mat_t) B, - int unit, - const TEMPLATE(T, ctx_t) ctx); - -void TEMPLATE(T, mat_solve_tril_recursive)(TEMPLATE(T, mat_t) X, - const TEMPLATE(T, mat_t) L, - const TEMPLATE(T, mat_t) B, - int unit, - const TEMPLATE(T, ctx_t) ctx); - void TEMPLATE(T, mat_solve_triu)(TEMPLATE(T, mat_t) X, const TEMPLATE(T, mat_t) U, const TEMPLATE(T, mat_t) B, int unit, const TEMPLATE(T, ctx_t) ctx); -void TEMPLATE(T, mat_solve_triu_classical)(TEMPLATE(T, mat_t) X, - const TEMPLATE(T, mat_t) U, - const TEMPLATE(T, mat_t) B, - int unit, - const TEMPLATE(T, ctx_t) ctx); -void TEMPLATE(T, mat_solve_triu_recursive)(TEMPLATE(T, mat_t) X, - const TEMPLATE(T, mat_t) U, - const TEMPLATE(T, mat_t) B, - int unit, - const TEMPLATE(T, ctx_t) ctx); - void TEMPLATE(T, mat_mul_vec)(TEMPLATE(T, struct) * c, const TEMPLATE(T, mat_t) A, const TEMPLATE(T, struct) * b, slong blen, diff --git a/src/fq_mat_templates/can_solve.c b/src/fq_mat_templates/can_solve.c index 9b0a2b84f4..f07bd0d0ab 100644 --- a/src/fq_mat_templates/can_solve.c +++ b/src/fq_mat_templates/can_solve.c @@ -11,132 +11,29 @@ */ #ifdef T + +#include "gr.h" +#include "gr_mat.h" #include "templates.h" int TEMPLATE(T, mat_can_solve)(TEMPLATE(T, mat_t) X, const TEMPLATE(T, mat_t) A, const TEMPLATE(T, mat_t) B, const TEMPLATE(T, ctx_t) ctx) { - slong i, j, k, col, *pivots, rank, *perm; - TEMPLATE(T, mat_t) LU, LU2, PB; - int result = 1; - - if (X->r != A->c || X->c != B->c) - return 0; - - if (A->r == 0 || B->c == 0) - { - TEMPLATE(T, mat_zero)(X, ctx); - - return 1; - } - - if (A->c == 0) - { - TEMPLATE(T, mat_zero)(X, ctx); - - return TEMPLATE(T, mat_is_zero)(B, ctx); - } - - TEMPLATE(T, mat_init_set)(LU, A, ctx); - perm = flint_malloc(sizeof(slong) * A->r); - for (i = 0; i < A->r; i++) - perm[i] = i; - - rank = TEMPLATE(T, mat_lu)(perm, LU, 0, ctx); - - TEMPLATE(T, mat_init)(PB, B->r, B->c, ctx); - for (i = 0; i < B->r; i++) - _TEMPLATE(T, vec_set)(TEMPLATE(T, mat_entry)(PB, i, 0), TEMPLATE(T, mat_entry)(B, perm[i], 0), B->c, ctx); - - TEMPLATE(T, mat_init)(LU2, rank, rank, ctx); - - pivots = flint_malloc(sizeof(slong)*rank); - - col = 0; - for (i = 0; i < rank; i++) - { - while (TEMPLATE(T, is_zero)(TEMPLATE(T, mat_entry)(LU, i, col), ctx)) - col++; - - pivots[i] = col; - - for (j = 0; j < rank; j++) - TEMPLATE(T, mat_entry_set)(LU2, j, i, TEMPLATE(T, mat_entry)(LU, j, col), ctx); - - col++; - } - - X->r = rank; - PB->r = rank; - LU->r = rank; - TEMPLATE(T, mat_solve_tril)(X, LU, PB, 1, ctx); - LU->r = A->r; - - if (A->r > rank) - { - TEMPLATE(T, mat_t) P; - - LU->entries += rank * LU->stride; - LU->r = A->r - rank; - X->r = LU->c; - - TEMPLATE(T, mat_init)(P, LU->r, B->c, ctx); - - TEMPLATE(T, mat_mul)(P, LU, X, ctx); - - PB->r = LU->r; - PB->entries += rank * PB->stride; - - result = TEMPLATE(T, mat_equal)(P, PB, ctx); - - PB->entries -= rank * PB->stride; - TEMPLATE(T, mat_clear)(P, ctx); - - LU->entries -= rank * LU->stride; - - if (!result) - { - X->r = A->c; - TEMPLATE(T, mat_zero)(X, ctx); - goto cleanup; - } - } - - TEMPLATE(T, mat_solve_triu)(X, LU2, X, 0, ctx); - - X->r = A->c; - - k = rank - 1; - for (i = A->c - 1; i >= 0; i--) - { - if (k < 0 || i != pivots[k]) - { - for (j = 0; j < B->c; j++) - TEMPLATE(T, zero)(TEMPLATE(T, mat_entry)(X, i, j), ctx); - } else - { - for (j = 0; j < B->c; j++) - TEMPLATE(T, mat_entry_set)(X, i, j, TEMPLATE(T, mat_entry)(X, k, j), ctx); - - k--; - } - } - -cleanup: - - TEMPLATE(T, mat_clear)(LU2, ctx); + gr_ctx_t gr_ctx; + int status; - PB->r = B->r; - TEMPLATE(T, mat_clear)(PB, ctx); + TEMPLATE3(_gr_ctx_init, T, from_ref)(gr_ctx, ctx); - LU->r = A->r; - TEMPLATE(T, mat_clear)(LU, ctx); - flint_free(perm); + status = gr_mat_solve_field((gr_mat_struct *) X, + (const gr_mat_struct *) A, + (const gr_mat_struct *) B, gr_ctx); - flint_free(pivots); + /* should not happen */ + if (status == GR_UNABLE) + flint_abort(); - return result; + return status == GR_SUCCESS; } #endif diff --git a/src/fq_mat_templates/lu.c b/src/fq_mat_templates/lu.c index f354ba0456..4fb4db353f 100644 --- a/src/fq_mat_templates/lu.c +++ b/src/fq_mat_templates/lu.c @@ -12,13 +12,19 @@ #ifdef T +#include "gr.h" +#include "gr_mat.h" #include "templates.h" slong TEMPLATE(T, mat_lu) (slong * P, TEMPLATE(T, mat_t) A, int rank_check, const TEMPLATE(T, ctx_t) ctx) { - return TEMPLATE(T, mat_lu_recursive) (P, A, rank_check, ctx); + gr_ctx_t gr_ctx; + slong rank; + TEMPLATE3(_gr_ctx_init, T, from_ref)(gr_ctx, ctx); + GR_MUST_SUCCEED(gr_mat_lu(&rank, P, (gr_mat_struct *) A, (const gr_mat_struct *) A, rank_check, gr_ctx)); + return rank; } diff --git a/src/fq_mat_templates/lu_classical.c b/src/fq_mat_templates/lu_classical.c deleted file mode 100644 index 8e04018199..0000000000 --- a/src/fq_mat_templates/lu_classical.c +++ /dev/null @@ -1,106 +0,0 @@ -/* - Copyright (C) 2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "templates.h" - - -static inline int -TEMPLATE(T, mat_pivot) (TEMPLATE(T, mat_t) A, slong * P, slong start_row, - slong col, const TEMPLATE(T, ctx_t) ctx) -{ - slong j; - - if (!TEMPLATE(T, is_zero) - (TEMPLATE(T, mat_entry) (A, start_row, col), ctx)) - return 1; - - for (j = start_row + 1; j < A->r; j++) - { - if (!TEMPLATE(T, is_zero) (TEMPLATE(T, mat_entry) (A, j, col), ctx)) - { - TEMPLATE(T, mat_swap_rows) (A, P, j, start_row, ctx); - return -1; - } - } - return 0; -} - - -slong -TEMPLATE(T, mat_lu_classical) (slong * P, - TEMPLATE(T, mat_t) A, - int rank_check, const TEMPLATE(T, ctx_t) ctx) -{ - TEMPLATE(T, t) d, e, neg_e; - slong i, m, n, rank, length, row, col; - - m = A->r; - n = A->c; - - rank = row = col = 0; - - for (i = 0; i < m; i++) - P[i] = i; - - TEMPLATE(T, init) (d, ctx); - TEMPLATE(T, init) (e, ctx); - TEMPLATE(T, init) (neg_e, ctx); - - while (row < m && col < n) - { - if (TEMPLATE(T, mat_pivot) (A, P, row, col, ctx) == 0) - { - if (rank_check) - { - rank = 0; - goto cleanup; - } - col++; - continue; - } - - rank++; - - TEMPLATE(T, inv) (d, TEMPLATE(T, mat_entry) (A, row, col), ctx); - - length = n - col - 1; - - for (i = row + 1; i < m; i++) - { - TEMPLATE(T, mul) (e, TEMPLATE(T, mat_entry) (A, i, col), d, ctx); - if (length != 0) - { - TEMPLATE(T, neg) (neg_e, e, ctx); - _TEMPLATE3(T, vec_scalar_addmul, T) (TEMPLATE(T, mat_entry) (A, i, col + 1), - TEMPLATE(T, mat_entry) (A, row, col + 1), - length, neg_e, ctx); - } - - TEMPLATE(T, zero) (TEMPLATE(T, mat_entry) (A, i, col), ctx); - TEMPLATE(T, set) (TEMPLATE(T, mat_entry) (A, i, rank - 1), e, ctx); - } - row++; - col++; - } - -cleanup: - TEMPLATE(T, clear) (d, ctx); - TEMPLATE(T, clear) (e, ctx); - TEMPLATE(T, clear) (neg_e, ctx); - - return rank; -} - - -#endif diff --git a/src/fq_mat_templates/lu_recursive.c b/src/fq_mat_templates/lu_recursive.c deleted file mode 100644 index 1b1e2aa814..0000000000 --- a/src/fq_mat_templates/lu_recursive.c +++ /dev/null @@ -1,139 +0,0 @@ -/* - Copyright (C) 2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include -#include "templates.h" - - -static void -_apply_permutation(slong * AP, TEMPLATE(T, mat_t) A, const slong * P, - slong num_rows, slong row_offset, slong num_cols, slong col_offset) -{ - if (num_rows != 0) - { - TEMPLATE(T, struct) * Atmp; - slong * APtmp; - slong i; - - /* todo: reduce memory allocation */ - Atmp = flint_malloc(sizeof(TEMPLATE(T, struct)) * num_rows * num_cols); - APtmp = flint_malloc(sizeof(slong) * num_rows); - - for (i = 0; i < num_rows; i++) - memcpy(Atmp + i * num_cols, TEMPLATE(T, mat_entry) (A, P[i] + row_offset, col_offset), num_cols * sizeof(TEMPLATE(T, struct))); - for (i = 0; i < num_rows; i++) - memcpy(TEMPLATE(T, mat_entry) (A, i + row_offset, col_offset), Atmp + i * num_cols, num_cols * sizeof(TEMPLATE(T, struct))); - - for (i = 0; i < num_rows; i++) APtmp[i] = AP[P[i] + row_offset]; - for (i = 0; i < num_rows; i++) AP[i + row_offset] = APtmp[i]; - - flint_free(Atmp); - flint_free(APtmp); - } -} - - -slong -TEMPLATE(T, mat_lu_recursive) (slong * P, - TEMPLATE(T, mat_t) A, - int rank_check, const TEMPLATE(T, ctx_t) ctx) -{ - - slong i, j, m, n, r1, r2, n1; - TEMPLATE(T, mat_t) A0, A1, A00, A01, A10, A11; - slong *P1; - - m = A->r; - n = A->c; - - if (m < TEMPLATE(CAP_T, MAT_LU_RECURSIVE_CUTOFF) - || n < TEMPLATE(CAP_T, MAT_LU_RECURSIVE_CUTOFF)) - { - r1 = TEMPLATE(T, mat_lu_classical) (P, A, rank_check, ctx); - return r1; - } - - n1 = n / 2; - - for (i = 0; i < m; i++) - P[i] = i; - - P1 = flint_malloc(sizeof(slong) * m); - TEMPLATE(T, mat_window_init) (A0, A, 0, 0, m, n1, ctx); - TEMPLATE(T, mat_window_init) (A1, A, 0, n1, m, n, ctx); - - r1 = TEMPLATE(T, mat_lu) (P1, A0, rank_check, ctx); - - if (rank_check && (r1 != n1)) - { - flint_free(P1); - TEMPLATE(T, mat_window_clear) (A0, ctx); - TEMPLATE(T, mat_window_clear) (A1, ctx); - return 0; - } - - if (r1 != 0) - { - _apply_permutation(P, A, P1, m, 0, n - n1, n1); - } - - TEMPLATE(T, mat_window_init) (A00, A, 0, 0, r1, r1, ctx); - TEMPLATE(T, mat_window_init) (A10, A, r1, 0, m, r1, ctx); - TEMPLATE(T, mat_window_init) (A01, A, 0, n1, r1, n, ctx); - TEMPLATE(T, mat_window_init) (A11, A, r1, n1, m, n, ctx); - - if (r1 != 0) - { - TEMPLATE(T, mat_solve_tril) (A01, A00, A01, 1, ctx); - TEMPLATE(T, mat_submul) (A11, A11, A10, A01, ctx); - } - - r2 = TEMPLATE(T, mat_lu) (P1, A11, rank_check, ctx); - - if (rank_check && (r1 + r2 < FLINT_MIN(m, n))) - { - r1 = r2 = 0; - } - else - { - _apply_permutation(P, A, P1, m - r1, r1, n1, 0); - - /* Compress L */ - if (r1 != n1) - { - for (i = 0; i < m - r1; i++) - { - TEMPLATE(T, struct) * row = TEMPLATE(T, mat_entry) (A, r1 + i, 0); - for (j = 0; j < FLINT_MIN(i, r2); j++) - { - TEMPLATE(T, set) (row + r1 + j, row + n1 + j, ctx); - TEMPLATE(T, zero) (row + n1 + j, ctx); - } - } - } - } - - flint_free(P1); - TEMPLATE(T, mat_window_clear) (A00, ctx); - TEMPLATE(T, mat_window_clear) (A01, ctx); - TEMPLATE(T, mat_window_clear) (A10, ctx); - TEMPLATE(T, mat_window_clear) (A11, ctx); - TEMPLATE(T, mat_window_clear) (A0, ctx); - TEMPLATE(T, mat_window_clear) (A1, ctx); - - return r1 + r2; -} - - -#endif diff --git a/src/fq_mat_templates/solve.c b/src/fq_mat_templates/solve.c index c8b0b60a82..46aaa64a5c 100644 --- a/src/fq_mat_templates/solve.c +++ b/src/fq_mat_templates/solve.c @@ -10,49 +10,22 @@ */ #ifdef T + +#include "gr.h" +#include "gr_mat.h" #include "templates.h" int TEMPLATE(T, mat_solve)(TEMPLATE(T, mat_t) X, const TEMPLATE(T, mat_t) A, const TEMPLATE(T, mat_t) B, const TEMPLATE(T, ctx_t) ctx) { - slong i, rank, *perm; - TEMPLATE(T, mat_t) LU; - int result; - - if (A->r == 0 || B->c == 0) - return 1; - - TEMPLATE(T, mat_init_set)(LU, A, ctx); - perm = flint_malloc(sizeof(slong) * A->r); - for (i = 0; i < A->r; i++) - perm[i] = i; - - rank = TEMPLATE(T, mat_lu)(perm, LU, 1, ctx); - - if (rank == A->r) - { - TEMPLATE(T, mat_t) PB; - TEMPLATE(T, mat_init)(PB, B->r, B->c, ctx); - for (i = 0; i < A->r; i++) - _TEMPLATE(T, vec_set)(TEMPLATE(T, mat_entry)(PB, i, 0), - TEMPLATE(T, mat_entry)(B, perm[i], 0), B->c, ctx); - - TEMPLATE(T, mat_solve_tril)(X, LU, PB, 1, ctx); - TEMPLATE(T, mat_solve_triu)(X, LU, X, 0, ctx); - - TEMPLATE(T, mat_clear)(PB, ctx); - result = 1; - } - else - { - result = 0; - } - - TEMPLATE(T, mat_clear)(LU, ctx); - flint_free(perm); - - return result; + gr_ctx_t gr_ctx; + int status; + TEMPLATE3(_gr_ctx_init, T, from_ref)(gr_ctx, ctx); + status = gr_mat_nonsingular_solve_lu((gr_mat_struct *) X, + (const gr_mat_struct *) A, + (const gr_mat_struct *) B, gr_ctx); + return (status == GR_SUCCESS); } #endif diff --git a/src/fq_mat_templates/solve_tril.c b/src/fq_mat_templates/solve_tril.c index 00a007489e..afa6976565 100644 --- a/src/fq_mat_templates/solve_tril.c +++ b/src/fq_mat_templates/solve_tril.c @@ -12,6 +12,8 @@ #ifdef T +#include "gr.h" +#include "gr_mat.h" #include "templates.h" void @@ -19,16 +21,11 @@ TEMPLATE(T, mat_solve_tril) (TEMPLATE(T, mat_t) X, const TEMPLATE(T, mat_t) L, const TEMPLATE(T, mat_t) B, int unit, const TEMPLATE(T, ctx_t) ctx) { - if (B->r < TEMPLATE(CAP_T, MAT_SOLVE_TRI_ROWS_CUTOFF) || - B->c < TEMPLATE(CAP_T, MAT_SOLVE_TRI_COLS_CUTOFF)) - { - TEMPLATE(T, mat_solve_tril_classical) (X, L, B, unit, ctx); - } - else - { - TEMPLATE(T, mat_solve_tril_recursive) (X, L, B, unit, ctx); - } + gr_ctx_t gr_ctx; + TEMPLATE3(_gr_ctx_init, T, from_ref)(gr_ctx, ctx); + GR_MUST_SUCCEED(gr_mat_nonsingular_solve_tril((gr_mat_struct *) X, + (const gr_mat_struct *) L, + (const gr_mat_struct *) B, unit, gr_ctx)); } - #endif diff --git a/src/fq_mat_templates/solve_tril_classical.c b/src/fq_mat_templates/solve_tril_classical.c deleted file mode 100644 index 58258bd648..0000000000 --- a/src/fq_mat_templates/solve_tril_classical.c +++ /dev/null @@ -1,67 +0,0 @@ -/* - Copyright (C) 2010,2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "templates.h" - -void -TEMPLATE(T, mat_solve_tril_classical) (TEMPLATE(T, mat_t) X, - const TEMPLATE(T, mat_t) L, - const TEMPLATE(T, mat_t) B, - int unit, const TEMPLATE(T, ctx_t) ctx) -{ - slong i, j, n, m; - TEMPLATE(T, struct) * inv, *tmp; - - n = L->r; - m = B->c; - - if (!unit) - { - inv = _TEMPLATE(T, vec_init) (n, ctx); - for (i = 0; i < n; i++) - TEMPLATE(T, inv) (inv + i, TEMPLATE(T, mat_entry) (L, i, i), ctx); - } - else - inv = NULL; - - tmp = _TEMPLATE(T, vec_init) (n, ctx); - - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - TEMPLATE(T, set) (tmp + j, TEMPLATE(T, mat_entry) (X, j, i), ctx); - - for (j = 0; j < n; j++) - { - TEMPLATE(T, t) s; - TEMPLATE(T, init) (s, ctx); - _TEMPLATE(T, vec_dot) (s, TEMPLATE(T, mat_entry) (L, j, 0), tmp, j, ctx); - TEMPLATE(T, sub) (s, TEMPLATE(T, mat_entry) (B, j, i), s, ctx); - if (!unit) - TEMPLATE(T, mul) (s, s, inv + j, ctx); - TEMPLATE(T, set) (tmp + j, s, ctx); - TEMPLATE(T, clear) (s, ctx); - } - - for (j = 0; j < n; j++) - TEMPLATE(T, mat_entry_set) (X, j, i, tmp + j, ctx); - } - - _TEMPLATE(T, vec_clear) (tmp, n, ctx); - if (!unit) - _TEMPLATE(T, vec_clear) (inv, n, ctx); -} - - -#endif diff --git a/src/fq_mat_templates/solve_tril_recursive.c b/src/fq_mat_templates/solve_tril_recursive.c deleted file mode 100644 index a7ef5bfd46..0000000000 --- a/src/fq_mat_templates/solve_tril_recursive.c +++ /dev/null @@ -1,62 +0,0 @@ -/* - Copyright (C) 2010,2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "templates.h" - -void -TEMPLATE(T, mat_solve_tril_recursive) (TEMPLATE(T, mat_t) X, - const TEMPLATE(T, mat_t) L, - const TEMPLATE(T, mat_t) B, - int unit, const TEMPLATE(T, ctx_t) ctx) -{ - TEMPLATE(T, mat_t) LA, LC, LD, XX, XY, BX, BY; - slong r, n, m; - - n = L->r; - m = B->c; - r = n / 2; - - if (n == 0 || m == 0) - return; - - /* - Denoting inv(M) by M^, we have: - - [A 0]^ [X] == [A^ 0 ] [X] == [A^ X] - [C D] [Y] == [-D^ C A^ D^] [Y] == [D^ (Y - C A^ X)] - */ - - TEMPLATE(T, mat_window_init) (LA, L, 0, 0, r, r, ctx); - TEMPLATE(T, mat_window_init) (LC, L, r, 0, n, r, ctx); - TEMPLATE(T, mat_window_init) (LD, L, r, r, n, n, ctx); - TEMPLATE(T, mat_window_init) (BX, B, 0, 0, r, m, ctx); - TEMPLATE(T, mat_window_init) (BY, B, r, 0, n, m, ctx); - TEMPLATE(T, mat_window_init) (XX, X, 0, 0, r, m, ctx); - TEMPLATE(T, mat_window_init) (XY, X, r, 0, n, m, ctx); - - TEMPLATE(T, mat_solve_tril) (XX, LA, BX, unit, ctx); - TEMPLATE(T, mat_submul) (XY, BY, LC, XX, ctx); - TEMPLATE(T, mat_solve_tril) (XY, LD, XY, unit, ctx); - - TEMPLATE(T, mat_window_clear) (LA, ctx); - TEMPLATE(T, mat_window_clear) (LC, ctx); - TEMPLATE(T, mat_window_clear) (LD, ctx); - TEMPLATE(T, mat_window_clear) (BX, ctx); - TEMPLATE(T, mat_window_clear) (BY, ctx); - TEMPLATE(T, mat_window_clear) (XX, ctx); - TEMPLATE(T, mat_window_clear) (XY, ctx); -} - - -#endif diff --git a/src/fq_mat_templates/solve_triu.c b/src/fq_mat_templates/solve_triu.c index 0e767e34c1..2403de21fd 100644 --- a/src/fq_mat_templates/solve_triu.c +++ b/src/fq_mat_templates/solve_triu.c @@ -12,6 +12,8 @@ #ifdef T +#include "gr.h" +#include "gr_mat.h" #include "templates.h" void @@ -19,16 +21,11 @@ TEMPLATE(T, mat_solve_triu) (TEMPLATE(T, mat_t) X, const TEMPLATE(T, mat_t) U, const TEMPLATE(T, mat_t) B, int unit, const TEMPLATE(T, ctx_t) ctx) { - if (B->r < TEMPLATE(CAP_T, MAT_SOLVE_TRI_ROWS_CUTOFF) || - B->c < TEMPLATE(CAP_T, MAT_SOLVE_TRI_COLS_CUTOFF)) - { - TEMPLATE(T, mat_solve_triu_classical) (X, U, B, unit, ctx); - } - else - { - TEMPLATE(T, mat_solve_triu_recursive) (X, U, B, unit, ctx); - } + gr_ctx_t gr_ctx; + TEMPLATE3(_gr_ctx_init, T, from_ref)(gr_ctx, ctx); + GR_MUST_SUCCEED(gr_mat_nonsingular_solve_triu((gr_mat_struct *) X, + (const gr_mat_struct *) U, + (const gr_mat_struct *) B, unit, gr_ctx)); } - #endif diff --git a/src/fq_mat_templates/solve_triu_classical.c b/src/fq_mat_templates/solve_triu_classical.c deleted file mode 100644 index eb405ea3d2..0000000000 --- a/src/fq_mat_templates/solve_triu_classical.c +++ /dev/null @@ -1,68 +0,0 @@ -/* - Copyright (C) 2010,2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "templates.h" - -void -TEMPLATE(T, mat_solve_triu_classical) (TEMPLATE(T, mat_t) X, - const TEMPLATE(T, mat_t) U, - const TEMPLATE(T, mat_t) B, - int unit, const TEMPLATE(T, ctx_t) ctx) -{ - slong i, j, n, m; - TEMPLATE(T, struct) * inv, *tmp; - - n = U->r; - m = B->c; - - if (!unit) - { - inv = _TEMPLATE(T, vec_init) (n, ctx); - for (i = 0; i < n; i++) - TEMPLATE(T, inv) (inv + i, TEMPLATE(T, mat_entry) (U, i, i), ctx); - } - else - inv = NULL; - - tmp = _TEMPLATE(T, vec_init) (n, ctx); - - for (i = 0; i < m; i++) - { - for (j = 0; j < n; j++) - TEMPLATE(T, set) (tmp + j, TEMPLATE(T, mat_entry) (X, j, i), ctx); - - for (j = n - 1; j >= 0; j--) - { - TEMPLATE(T, t) s; - TEMPLATE(T, init) (s, ctx); - _TEMPLATE(T, vec_dot) (s, TEMPLATE(T, mat_entry) (U, j, j + 1), tmp + j + 1, - n - j - 1, ctx); - TEMPLATE(T, sub) (s, TEMPLATE(T, mat_entry) (B, j, i), s, ctx); - if (!unit) - TEMPLATE(T, mul) (s, s, inv + j, ctx); - TEMPLATE(T, set) (tmp + j, s, ctx); - TEMPLATE(T, clear) (s, ctx); - } - - for (j = 0; j < n; j++) - TEMPLATE(T, mat_entry_set) (X, j, i, tmp + j, ctx); - } - - _TEMPLATE(T, vec_clear) (tmp, n, ctx); - if (!unit) - _TEMPLATE(T, vec_clear) (inv, n, ctx); -} - - -#endif diff --git a/src/fq_mat_templates/solve_triu_recursive.c b/src/fq_mat_templates/solve_triu_recursive.c deleted file mode 100644 index d5937a7518..0000000000 --- a/src/fq_mat_templates/solve_triu_recursive.c +++ /dev/null @@ -1,62 +0,0 @@ -/* - Copyright (C) 2010,2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "templates.h" - -void -TEMPLATE(T, mat_solve_triu_recursive) (TEMPLATE(T, mat_t) X, - const TEMPLATE(T, mat_t) U, - const TEMPLATE(T, mat_t) B, - int unit, const TEMPLATE(T, ctx_t) ctx) -{ - TEMPLATE(T, mat_t) UA, UB, UD, XX, XY, BX, BY; - slong r, n, m; - - n = U->r; - m = B->c; - r = n / 2; - - if (n == 0 || m == 0) - return; - - /* - Denoting inv(M) by M^, we have: - - [A B]^ [X] == [A^ (X - B D^ Y)] - [0 D] [Y] == [ D^ Y ] - */ - - TEMPLATE(T, mat_window_init) (UA, U, 0, 0, r, r, ctx); - TEMPLATE(T, mat_window_init) (UB, U, 0, r, r, n, ctx); - TEMPLATE(T, mat_window_init) (UD, U, r, r, n, n, ctx); - TEMPLATE(T, mat_window_init) (BX, B, 0, 0, r, m, ctx); - TEMPLATE(T, mat_window_init) (BY, B, r, 0, n, m, ctx); - TEMPLATE(T, mat_window_init) (XX, X, 0, 0, r, m, ctx); - TEMPLATE(T, mat_window_init) (XY, X, r, 0, n, m, ctx); - - TEMPLATE(T, mat_solve_triu) (XY, UD, BY, unit, ctx); - TEMPLATE(T, mat_submul) (XX, BX, UB, XY, ctx); - TEMPLATE(T, mat_solve_triu) (XX, UA, XX, unit, ctx); - - TEMPLATE(T, mat_window_clear) (UA, ctx); - TEMPLATE(T, mat_window_clear) (UB, ctx); - TEMPLATE(T, mat_window_clear) (UD, ctx); - TEMPLATE(T, mat_window_clear) (BX, ctx); - TEMPLATE(T, mat_window_clear) (BY, ctx); - TEMPLATE(T, mat_window_clear) (XX, ctx); - TEMPLATE(T, mat_window_clear) (XY, ctx); -} - - -#endif diff --git a/src/fq_mat_templates/test/t-lu_classical.c b/src/fq_mat_templates/test/t-lu.c similarity index 94% rename from src/fq_mat_templates/test/t-lu_classical.c rename to src/fq_mat_templates/test/t-lu.c index 6c53a8dfc3..661b6f322a 100644 --- a/src/fq_mat_templates/test/t-lu_classical.c +++ b/src/fq_mat_templates/test/t-lu.c @@ -16,7 +16,6 @@ #include "test_helpers.h" #include "templates.h" -/* Defined in t-lu_classical.c and t-lu_recursive.c */ #ifndef perm #define perm perm void @@ -37,7 +36,6 @@ perm(TEMPLATE(T, mat_t) A, slong * P) } #endif -/* Defined in t-lu_classical.c and t-lu_recursive.c */ #ifndef check #define check check void @@ -105,7 +103,7 @@ check(slong * P, TEMPLATE(T, mat_t) LU, const TEMPLATE(T, mat_t) A, slong rank, } #endif -TEST_TEMPLATE_FUNCTION_START(T, mat_lu_classical, state) +TEST_TEMPLATE_FUNCTION_START(T, mat_lu, state) { slong i; @@ -136,7 +134,7 @@ TEST_TEMPLATE_FUNCTION_START(T, mat_lu_classical, state) TEMPLATE(T, mat_init_set) (LU, A, ctx); P = flint_malloc(sizeof(slong) * m); - rank = TEMPLATE(T, mat_lu_classical) (P, LU, 0, ctx); + rank = TEMPLATE(T, mat_lu) (P, LU, 0, ctx); if (r != rank) { diff --git a/src/fq_mat_templates/test/t-lu_recursive.c b/src/fq_mat_templates/test/t-lu_recursive.c deleted file mode 100644 index 0482000636..0000000000 --- a/src/fq_mat_templates/test/t-lu_recursive.c +++ /dev/null @@ -1,165 +0,0 @@ -/* - Copyright (C) 2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include -#include "test_helpers.h" -#include "templates.h" - -/* Defined in t-lu_classical.c and t-lu_recursive.c */ -#ifndef perm -#define perm perm -void -perm(TEMPLATE(T, mat_t) A, slong * P) -{ - slong i; - TEMPLATE(T, struct) * tmp; - - if (A->c == 0 || A->r == 0) - return; - - tmp = flint_malloc(sizeof(TEMPLATE(T, struct)) * A->r * A->c); - - for (i = 0; i < A->r; i++) memcpy(tmp + P[i] * A->c, TEMPLATE(T, mat_entry)(A, i, 0), A->c * sizeof(TEMPLATE(T, struct))); - for (i = 0; i < A->r; i++) memcpy(TEMPLATE(T, mat_entry)(A, i, 0), tmp + i * A->c, A->c * sizeof(TEMPLATE(T, struct))); - - flint_free(tmp); -} -#endif - -/* Defined in t-lu_classical.c and t-lu_recursive.c */ -#ifndef check -#define check check -void -check(slong * P, TEMPLATE(T, mat_t) LU, const TEMPLATE(T, mat_t) A, slong rank, - const TEMPLATE(T, ctx_t) ctx) -{ - TEMPLATE(T, mat_t) B, L, U; - slong m, n, i, j; - - m = A->r; - n = A->c; - - TEMPLATE(T, mat_init) (B, m, n, ctx); - TEMPLATE(T, mat_init) (L, m, m, ctx); - TEMPLATE(T, mat_init) (U, m, n, ctx); - - rank = FLINT_ABS(rank); - - for (i = rank; i < FLINT_MIN(m, n); i++) - { - for (j = i; j < n; j++) - { - if (!TEMPLATE(T, is_zero) (TEMPLATE(T, mat_entry) (LU, i, j), ctx)) - { - printf("FAIL: wrong shape!\n"); - fflush(stdout); - flint_abort(); - } - } - } - - for (i = 0; i < m; i++) - { - for (j = 0; j < FLINT_MIN(i, n); j++) - TEMPLATE(T, mat_entry_set) (L, i, j, - TEMPLATE(T, mat_entry) (LU, i, j), - ctx); - if (i < rank) - TEMPLATE(T, one) (TEMPLATE(T, mat_entry) (L, i, i), ctx); - for (j = i; j < n; j++) - TEMPLATE(T, mat_entry_set) (U, i, j, - TEMPLATE(T, mat_entry) (LU, i, j), - ctx); - } - - TEMPLATE(T, mat_mul) (B, L, U, ctx); - perm(B, P); - - if (!TEMPLATE(T, mat_equal) (A, B, ctx)) - { - printf("FAIL\n"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("LU:\n"); - TEMPLATE(T, mat_print_pretty) (LU, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - fflush(stdout); - flint_abort(); - } - - TEMPLATE(T, mat_clear) (B, ctx); - TEMPLATE(T, mat_clear) (L, ctx); - TEMPLATE(T, mat_clear) (U, ctx); -} -#endif - -TEST_TEMPLATE_FUNCTION_START(T, mat_lu_recursive, state) -{ - slong i; - - for (i = 0; i < 10 * flint_test_multiplier(); i++) - { - TEMPLATE(T, ctx_t) ctx; - TEMPLATE(T, mat_t) A, LU; - - slong m, n, r, d, rank; - slong *P; - - TEMPLATE(T, ctx_init_randtest)(ctx, state, 3); - - m = n_randint(state, 20); - n = n_randint(state, 20); - - for (r = 0; r <= FLINT_MIN(m, n); r++) - { - TEMPLATE(T, mat_init) (A, m, n, ctx); - TEMPLATE(T, mat_randrank) (A, state, r, ctx); - - if (n_randint(state, 2)) - { - d = n_randint(state, 2 * m * n + 1); - TEMPLATE(T, mat_randops) (A, state, d, ctx); - } - - TEMPLATE(T, mat_init_set) (LU, A, ctx); - P = flint_malloc(sizeof(slong) * m); - - rank = TEMPLATE(T, mat_lu_recursive) (P, LU, 0, ctx); - - if (r != rank) - { - printf("FAIL:\n"); - printf("wrong rank!\n"); - printf("A:"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("LU:"); - TEMPLATE(T, mat_print_pretty) (LU, ctx); - fflush(stdout); - flint_abort(); - } - - check(P, LU, A, rank, ctx); - - TEMPLATE(T, mat_clear) (A, ctx); - TEMPLATE(T, mat_clear) (LU, ctx); - flint_free(P); - } - - TEMPLATE(T, ctx_clear) (ctx); - } - - TEST_FUNCTION_END(state); -} -#endif diff --git a/src/fq_mat_templates/test/t-solve_tril_classical.c b/src/fq_mat_templates/test/t-solve_tril_classical.c deleted file mode 100644 index dea03197d9..0000000000 --- a/src/fq_mat_templates/test/t-solve_tril_classical.c +++ /dev/null @@ -1,86 +0,0 @@ -/* - Copyright (C) 2010, 2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "test_helpers.h" -#include "templates.h" - -TEST_TEMPLATE_FUNCTION_START(T, mat_solve_tril_classical, state) -{ - slong i; - - for (i = 0; i < 5 * flint_test_multiplier(); i++) - { - TEMPLATE(T, ctx_t) ctx; - TEMPLATE(T, mat_t) A, X, B, Y; - - slong rows, cols; - int unit; - - TEMPLATE(T, ctx_init_randtest)(ctx, state, 3); - - rows = n_randint(state, 50); - cols = n_randint(state, 50); - unit = n_randint(state, 2); - - TEMPLATE(T, mat_init) (A, rows, rows, ctx); - TEMPLATE(T, mat_init) (B, rows, cols, ctx); - TEMPLATE(T, mat_init) (X, rows, cols, ctx); - TEMPLATE(T, mat_init) (Y, rows, cols, ctx); - - TEMPLATE(T, mat_randtril) (A, state, unit, ctx); - TEMPLATE(T, mat_randtest) (X, state, ctx); - TEMPLATE(T, mat_mul) (B, A, X, ctx); - - /* Check Y = A^(-1) * (A * X) = X */ - TEMPLATE(T, mat_solve_tril_classical) (Y, A, B, unit, ctx); - if (!TEMPLATE(T, mat_equal) (Y, X, ctx)) - { - printf("FAIL!\n"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("X:\n"); - TEMPLATE(T, mat_print_pretty) (X, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - printf("Y:\n"); - TEMPLATE(T, mat_print_pretty) (Y, ctx); - fflush(stdout); - flint_abort(); - } - - /* Check aliasing */ - TEMPLATE(T, mat_solve_tril_classical) (B, A, B, unit, ctx); - if (!TEMPLATE(T, mat_equal) (B, X, ctx)) - { - printf("FAIL!\n"); - printf("aliasing test failed"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - fflush(stdout); - flint_abort(); - } - - TEMPLATE(T, mat_clear) (A, ctx); - TEMPLATE(T, mat_clear) (B, ctx); - TEMPLATE(T, mat_clear) (X, ctx); - TEMPLATE(T, mat_clear) (Y, ctx); - - TEMPLATE(T, ctx_clear) (ctx); - } - - TEST_FUNCTION_END(state); -} -#endif diff --git a/src/fq_mat_templates/test/t-solve_tril_recursive.c b/src/fq_mat_templates/test/t-solve_tril_recursive.c deleted file mode 100644 index 180d577017..0000000000 --- a/src/fq_mat_templates/test/t-solve_tril_recursive.c +++ /dev/null @@ -1,86 +0,0 @@ -/* - Copyright (C) 2010, 2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "test_helpers.h" -#include "templates.h" - -TEST_TEMPLATE_FUNCTION_START(T, mat_solve_tril_recursive, state) -{ - slong i; - - for (i = 0; i < 5 * flint_test_multiplier(); i++) - { - TEMPLATE(T, ctx_t) ctx; - TEMPLATE(T, mat_t) A, X, B, Y; - - slong rows, cols; - int unit; - - TEMPLATE(T, ctx_init_randtest)(ctx, state, 3); - - rows = n_randint(state, 50); - cols = n_randint(state, 50); - unit = n_randint(state, 2); - - TEMPLATE(T, mat_init) (A, rows, rows, ctx); - TEMPLATE(T, mat_init) (B, rows, cols, ctx); - TEMPLATE(T, mat_init) (X, rows, cols, ctx); - TEMPLATE(T, mat_init) (Y, rows, cols, ctx); - - TEMPLATE(T, mat_randtril) (A, state, unit, ctx); - TEMPLATE(T, mat_randtest) (X, state, ctx); - TEMPLATE(T, mat_mul) (B, A, X, ctx); - - /* Check Y = A^(-1) * (A * X) = X */ - TEMPLATE(T, mat_solve_tril_recursive) (Y, A, B, unit, ctx); - if (!TEMPLATE(T, mat_equal) (Y, X, ctx)) - { - printf("FAIL!\n"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("X:\n"); - TEMPLATE(T, mat_print_pretty) (X, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - printf("Y:\n"); - TEMPLATE(T, mat_print_pretty) (Y, ctx); - fflush(stdout); - flint_abort(); - } - - /* Check aliasing */ - TEMPLATE(T, mat_solve_tril_recursive) (B, A, B, unit, ctx); - if (!TEMPLATE(T, mat_equal) (B, X, ctx)) - { - printf("FAIL!\n"); - printf("aliasing test failed"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - fflush(stdout); - flint_abort(); - } - - TEMPLATE(T, mat_clear) (A, ctx); - TEMPLATE(T, mat_clear) (B, ctx); - TEMPLATE(T, mat_clear) (X, ctx); - TEMPLATE(T, mat_clear) (Y, ctx); - - TEMPLATE(T, ctx_clear) (ctx); - } - - TEST_FUNCTION_END(state); -} -#endif diff --git a/src/fq_mat_templates/test/t-solve_triu_classical.c b/src/fq_mat_templates/test/t-solve_triu_classical.c deleted file mode 100644 index 811aebc60e..0000000000 --- a/src/fq_mat_templates/test/t-solve_triu_classical.c +++ /dev/null @@ -1,85 +0,0 @@ -/* - Copyright (C) 2010, 2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "test_helpers.h" -#include "templates.h" - -TEST_TEMPLATE_FUNCTION_START(T, mat_solve_triu_classical, state) -{ - slong i; - - for (i = 0; i < 5 * flint_test_multiplier(); i++) - { - TEMPLATE(T, ctx_t) ctx; - TEMPLATE(T, mat_t) A, X, B, Y; - slong rows, cols; - int unit; - - TEMPLATE(T, ctx_init_randtest)(ctx, state, 3); - - rows = n_randint(state, 50); - cols = n_randint(state, 50); - unit = n_randint(state, 2); - - TEMPLATE(T, mat_init) (A, rows, rows, ctx); - TEMPLATE(T, mat_init) (B, rows, cols, ctx); - TEMPLATE(T, mat_init) (X, rows, cols, ctx); - TEMPLATE(T, mat_init) (Y, rows, cols, ctx); - - TEMPLATE(T, mat_randtriu) (A, state, unit, ctx); - TEMPLATE(T, mat_randtest) (X, state, ctx); - TEMPLATE(T, mat_mul) (B, A, X, ctx); - - /* Check Y = A^(-1) * (A * X) = X */ - TEMPLATE(T, mat_solve_triu_classical) (Y, A, B, unit, ctx); - if (!TEMPLATE(T, mat_equal) (Y, X, ctx)) - { - printf("FAIL!\n"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("X:\n"); - TEMPLATE(T, mat_print_pretty) (X, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - printf("Y:\n"); - TEMPLATE(T, mat_print_pretty) (Y, ctx); - fflush(stdout); - flint_abort(); - } - - /* Check aliasing */ - TEMPLATE(T, mat_solve_triu_classical) (B, A, B, unit, ctx); - if (!TEMPLATE(T, mat_equal) (B, X, ctx)) - { - printf("FAIL!\n"); - printf("aliasing test failed"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - fflush(stdout); - flint_abort(); - } - - TEMPLATE(T, mat_clear) (A, ctx); - TEMPLATE(T, mat_clear) (B, ctx); - TEMPLATE(T, mat_clear) (X, ctx); - TEMPLATE(T, mat_clear) (Y, ctx); - - TEMPLATE(T, ctx_clear) (ctx); - } - - TEST_FUNCTION_END(state); -} -#endif diff --git a/src/fq_mat_templates/test/t-solve_triu_recursive.c b/src/fq_mat_templates/test/t-solve_triu_recursive.c deleted file mode 100644 index 49b8b9a74e..0000000000 --- a/src/fq_mat_templates/test/t-solve_triu_recursive.c +++ /dev/null @@ -1,85 +0,0 @@ -/* - Copyright (C) 2010, 2011 Fredrik Johansson - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#ifdef T - -#include "test_helpers.h" -#include "templates.h" - -TEST_TEMPLATE_FUNCTION_START(T, mat_solve_triu_recursive, state) -{ - slong i; - - for (i = 0; i < 5 * flint_test_multiplier(); i++) - { - TEMPLATE(T, ctx_t) ctx; - TEMPLATE(T, mat_t) A, X, B, Y; - slong rows, cols; - int unit; - - TEMPLATE(T, ctx_init_randtest)(ctx, state, 3); - - rows = n_randint(state, 50); - cols = n_randint(state, 50); - unit = n_randint(state, 2); - - TEMPLATE(T, mat_init) (A, rows, rows, ctx); - TEMPLATE(T, mat_init) (B, rows, cols, ctx); - TEMPLATE(T, mat_init) (X, rows, cols, ctx); - TEMPLATE(T, mat_init) (Y, rows, cols, ctx); - - TEMPLATE(T, mat_randtriu) (A, state, unit, ctx); - TEMPLATE(T, mat_randtest) (X, state, ctx); - TEMPLATE(T, mat_mul) (B, A, X, ctx); - - /* Check Y = A^(-1) * (A * X) = X */ - TEMPLATE(T, mat_solve_triu_recursive) (Y, A, B, unit, ctx); - if (!TEMPLATE(T, mat_equal) (Y, X, ctx)) - { - printf("FAIL!\n"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("X:\n"); - TEMPLATE(T, mat_print_pretty) (X, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - printf("Y:\n"); - TEMPLATE(T, mat_print_pretty) (Y, ctx); - fflush(stdout); - flint_abort(); - } - - /* Check aliasing */ - TEMPLATE(T, mat_solve_triu_recursive) (B, A, B, unit, ctx); - if (!TEMPLATE(T, mat_equal) (B, X, ctx)) - { - printf("FAIL!\n"); - printf("aliasing test failed"); - printf("A:\n"); - TEMPLATE(T, mat_print_pretty) (A, ctx); - printf("B:\n"); - TEMPLATE(T, mat_print_pretty) (B, ctx); - fflush(stdout); - flint_abort(); - } - - TEMPLATE(T, mat_clear) (A, ctx); - TEMPLATE(T, mat_clear) (B, ctx); - TEMPLATE(T, mat_clear) (X, ctx); - TEMPLATE(T, mat_clear) (Y, ctx); - - TEMPLATE(T, ctx_clear) (ctx); - } - - TEST_FUNCTION_END(state); -} -#endif diff --git a/src/fq_nmod/mat_templates.c b/src/fq_nmod/mat_templates.c index d7ec0e2daf..dfaaae6237 100644 --- a/src/fq_nmod/mat_templates.c +++ b/src/fq_nmod/mat_templates.c @@ -43,8 +43,6 @@ #include "fq_mat_templates/is_one.c" #include "fq_mat_templates/is_zero.c" #include "fq_mat_templates/lu.c" -#include "fq_mat_templates/lu_classical.c" -#include "fq_mat_templates/lu_recursive.c" #include "fq_mat_templates/mat_entry_set.c" #include "fq_mat_templates/mat_invert_cols.c" #include "fq_mat_templates/mat_swap_cols.c" @@ -71,11 +69,7 @@ #include "fq_mat_templates/similarity.c" #include "fq_mat_templates/solve.c" #include "fq_mat_templates/solve_tril.c" -#include "fq_mat_templates/solve_tril_classical.c" -#include "fq_mat_templates/solve_tril_recursive.c" #include "fq_mat_templates/solve_triu.c" -#include "fq_mat_templates/solve_triu_classical.c" -#include "fq_mat_templates/solve_triu_recursive.c" #include "fq_mat_templates/sub.c" #include "fq_mat_templates/submul.c" #include "fq_mat_templates/swap.c" diff --git a/src/fq_nmod_mat/test/main.c b/src/fq_nmod_mat/test/main.c index 941fd74ad1..d1039167f8 100644 --- a/src/fq_nmod_mat/test/main.c +++ b/src/fq_nmod_mat/test/main.c @@ -20,8 +20,7 @@ #include "t-inv.c" #include "t-invert_rows_cols.c" #include "t-is_zero.c" -#include "t-lu_classical.c" -#include "t-lu_recursive.c" +#include "t-lu.c" #include "t-minpoly.c" #include "t-mul.c" #include "t-mul_KS.c" @@ -34,11 +33,7 @@ #include "t-set_nmod_mat.c" #include "t-solve.c" #include "t-solve_tril.c" -#include "t-solve_tril_classical.c" -#include "t-solve_tril_recursive.c" #include "t-solve_triu.c" -#include "t-solve_triu_classical.c" -#include "t-solve_triu_recursive.c" #include "t-submul.c" #include "t-vec_mul.c" #include "t-window_init_clear.c" @@ -57,8 +52,7 @@ test_struct tests[] = TEST_FUNCTION(fq_nmod_mat_inv), TEST_FUNCTION(fq_nmod_mat_invert_rows_cols), TEST_FUNCTION(fq_nmod_mat_is_zero), - TEST_FUNCTION(fq_nmod_mat_lu_classical), - TEST_FUNCTION(fq_nmod_mat_lu_recursive), + TEST_FUNCTION(fq_nmod_mat_lu), TEST_FUNCTION(fq_nmod_mat_minpoly), TEST_FUNCTION(fq_nmod_mat_mul), TEST_FUNCTION(fq_nmod_mat_mul_KS), @@ -71,11 +65,7 @@ test_struct tests[] = TEST_FUNCTION(fq_nmod_mat_set_nmod_mat), TEST_FUNCTION(fq_nmod_mat_solve), TEST_FUNCTION(fq_nmod_mat_solve_tril), - TEST_FUNCTION(fq_nmod_mat_solve_tril_classical), - TEST_FUNCTION(fq_nmod_mat_solve_tril_recursive), TEST_FUNCTION(fq_nmod_mat_solve_triu), - TEST_FUNCTION(fq_nmod_mat_solve_triu_classical), - TEST_FUNCTION(fq_nmod_mat_solve_triu_recursive), TEST_FUNCTION(fq_nmod_mat_submul), TEST_FUNCTION(fq_nmod_mat_vec_mul), TEST_FUNCTION(fq_nmod_mat_window_init_clear), diff --git a/src/fq_nmod_mat/test/t-lu_classical.c b/src/fq_nmod_mat/test/t-lu.c similarity index 91% rename from src/fq_nmod_mat/test/t-lu_classical.c rename to src/fq_nmod_mat/test/t-lu.c index 59844f6507..1dc2ac614b 100644 --- a/src/fq_nmod_mat/test/t-lu_classical.c +++ b/src/fq_nmod_mat/test/t-lu.c @@ -18,6 +18,6 @@ #define T fq_nmod #define CAP_T FQ_NMOD -#include "fq_mat_templates/test/t-lu_classical.c" +#include "fq_mat_templates/test/t-lu.c" #undef CAP_T #undef T diff --git a/src/fq_nmod_mat/test/t-lu_recursive.c b/src/fq_nmod_mat/test/t-lu_recursive.c deleted file mode 100644 index 6f8604e9b3..0000000000 --- a/src/fq_nmod_mat/test/t-lu_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_nmod.h" -#include "fq_nmod_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_nmod -#define CAP_T FQ_NMOD -#include "fq_mat_templates/test/t-lu_recursive.c" -#undef CAP_T -#undef T diff --git a/src/fq_nmod_mat/test/t-solve_tril_classical.c b/src/fq_nmod_mat/test/t-solve_tril_classical.c deleted file mode 100644 index 80b78feb48..0000000000 --- a/src/fq_nmod_mat/test/t-solve_tril_classical.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_nmod.h" -#include "fq_nmod_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_nmod -#define CAP_T FQ_NMOD -#include "fq_mat_templates/test/t-solve_tril_classical.c" -#undef CAP_T -#undef T diff --git a/src/fq_nmod_mat/test/t-solve_tril_recursive.c b/src/fq_nmod_mat/test/t-solve_tril_recursive.c deleted file mode 100644 index 6fd834209e..0000000000 --- a/src/fq_nmod_mat/test/t-solve_tril_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_nmod.h" -#include "fq_nmod_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_nmod -#define CAP_T FQ_NMOD -#include "fq_mat_templates/test/t-solve_tril_recursive.c" -#undef CAP_T -#undef T diff --git a/src/fq_nmod_mat/test/t-solve_triu_classical.c b/src/fq_nmod_mat/test/t-solve_triu_classical.c deleted file mode 100644 index 6b4c60749f..0000000000 --- a/src/fq_nmod_mat/test/t-solve_triu_classical.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_nmod.h" -#include "fq_nmod_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_nmod -#define CAP_T FQ_NMOD -#include "fq_mat_templates/test/t-solve_triu_classical.c" -#undef CAP_T -#undef T diff --git a/src/fq_nmod_mat/test/t-solve_triu_recursive.c b/src/fq_nmod_mat/test/t-solve_triu_recursive.c deleted file mode 100644 index 4336d3a440..0000000000 --- a/src/fq_nmod_mat/test/t-solve_triu_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_nmod.h" -#include "fq_nmod_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_nmod -#define CAP_T FQ_NMOD -#include "fq_mat_templates/test/t-solve_triu_recursive.c" -#undef CAP_T -#undef T diff --git a/src/fq_zech/mat_templates.c b/src/fq_zech/mat_templates.c index 0f76cc8839..27cd851dca 100644 --- a/src/fq_zech/mat_templates.c +++ b/src/fq_zech/mat_templates.c @@ -42,8 +42,6 @@ #include "fq_mat_templates/is_one.c" #include "fq_mat_templates/is_zero.c" #include "fq_mat_templates/lu.c" -#include "fq_mat_templates/lu_classical.c" -#include "fq_mat_templates/lu_recursive.c" #include "fq_mat_templates/mat_entry_set.c" #include "fq_mat_templates/mat_invert_cols.c" #include "fq_mat_templates/mat_swap_cols.c" @@ -70,11 +68,7 @@ #include "fq_mat_templates/similarity.c" #include "fq_mat_templates/solve.c" #include "fq_mat_templates/solve_tril.c" -#include "fq_mat_templates/solve_tril_classical.c" -#include "fq_mat_templates/solve_tril_recursive.c" #include "fq_mat_templates/solve_triu.c" -#include "fq_mat_templates/solve_triu_classical.c" -#include "fq_mat_templates/solve_triu_recursive.c" #include "fq_mat_templates/sub.c" #include "fq_mat_templates/submul.c" #include "fq_mat_templates/swap.c" diff --git a/src/fq_zech_mat/test/main.c b/src/fq_zech_mat/test/main.c index 9c0488bcc7..217a2c9cb3 100644 --- a/src/fq_zech_mat/test/main.c +++ b/src/fq_zech_mat/test/main.c @@ -19,8 +19,7 @@ #include "t-equal.c" #include "t-inv.c" #include "t-is_zero.c" -#include "t-lu_classical.c" -#include "t-lu_recursive.c" +#include "t-lu.c" #include "t-minpoly.c" #include "t-mul.c" #include "t-mul_KS.c" @@ -33,11 +32,7 @@ #include "t-set_nmod_mat.c" #include "t-solve.c" #include "t-solve_tril.c" -#include "t-solve_tril_classical.c" -#include "t-solve_tril_recursive.c" #include "t-solve_triu.c" -#include "t-solve_triu_classical.c" -#include "t-solve_triu_recursive.c" #include "t-submul.c" #include "t-vec_mul.c" #include "t-window_init_clear.c" @@ -55,8 +50,7 @@ test_struct tests[] = TEST_FUNCTION(fq_zech_mat_equal), TEST_FUNCTION(fq_zech_mat_inv), TEST_FUNCTION(fq_zech_mat_is_zero), - TEST_FUNCTION(fq_zech_mat_lu_classical), - TEST_FUNCTION(fq_zech_mat_lu_recursive), + TEST_FUNCTION(fq_zech_mat_lu), TEST_FUNCTION(fq_zech_mat_minpoly), TEST_FUNCTION(fq_zech_mat_mul), TEST_FUNCTION(fq_zech_mat_mul_KS), @@ -69,11 +63,7 @@ test_struct tests[] = TEST_FUNCTION(fq_zech_mat_set_nmod_mat), TEST_FUNCTION(fq_zech_mat_solve), TEST_FUNCTION(fq_zech_mat_solve_tril), - TEST_FUNCTION(fq_zech_mat_solve_tril_classical), - TEST_FUNCTION(fq_zech_mat_solve_tril_recursive), TEST_FUNCTION(fq_zech_mat_solve_triu), - TEST_FUNCTION(fq_zech_mat_solve_triu_classical), - TEST_FUNCTION(fq_zech_mat_solve_triu_recursive), TEST_FUNCTION(fq_zech_mat_submul), TEST_FUNCTION(fq_zech_mat_vec_mul), TEST_FUNCTION(fq_zech_mat_window_init_clear), diff --git a/src/fq_zech_mat/test/t-lu_classical.c b/src/fq_zech_mat/test/t-lu.c similarity index 91% rename from src/fq_zech_mat/test/t-lu_classical.c rename to src/fq_zech_mat/test/t-lu.c index 6c2df4c498..50f73cc418 100644 --- a/src/fq_zech_mat/test/t-lu_classical.c +++ b/src/fq_zech_mat/test/t-lu.c @@ -18,6 +18,6 @@ #define T fq_zech #define CAP_T FQ_ZECH -#include "fq_mat_templates/test/t-lu_classical.c" +#include "fq_mat_templates/test/t-lu.c" #undef CAP_T #undef T diff --git a/src/fq_zech_mat/test/t-lu_recursive.c b/src/fq_zech_mat/test/t-lu_recursive.c deleted file mode 100644 index 67485fa75b..0000000000 --- a/src/fq_zech_mat/test/t-lu_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_zech.h" -#include "fq_zech_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_zech -#define CAP_T FQ_ZECH -#include "fq_mat_templates/test/t-lu_recursive.c" -#undef CAP_T -#undef T diff --git a/src/fq_zech_mat/test/t-solve_tril_classical.c b/src/fq_zech_mat/test/t-solve_tril_classical.c deleted file mode 100644 index 250fa8d91a..0000000000 --- a/src/fq_zech_mat/test/t-solve_tril_classical.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_zech.h" -#include "fq_zech_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_zech -#define CAP_T FQ_ZECH -#include "fq_mat_templates/test/t-solve_tril_classical.c" -#undef CAP_T -#undef T diff --git a/src/fq_zech_mat/test/t-solve_tril_recursive.c b/src/fq_zech_mat/test/t-solve_tril_recursive.c deleted file mode 100644 index fd21e68619..0000000000 --- a/src/fq_zech_mat/test/t-solve_tril_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_zech.h" -#include "fq_zech_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_zech -#define CAP_T FQ_ZECH -#include "fq_mat_templates/test/t-solve_tril_recursive.c" -#undef CAP_T -#undef T diff --git a/src/fq_zech_mat/test/t-solve_triu_classical.c b/src/fq_zech_mat/test/t-solve_triu_classical.c deleted file mode 100644 index 42a1286ba4..0000000000 --- a/src/fq_zech_mat/test/t-solve_triu_classical.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_zech.h" -#include "fq_zech_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_zech -#define CAP_T FQ_ZECH -#include "fq_mat_templates/test/t-solve_triu_classical.c" -#undef CAP_T -#undef T diff --git a/src/fq_zech_mat/test/t-solve_triu_recursive.c b/src/fq_zech_mat/test/t-solve_triu_recursive.c deleted file mode 100644 index 8714eb930b..0000000000 --- a/src/fq_zech_mat/test/t-solve_triu_recursive.c +++ /dev/null @@ -1,23 +0,0 @@ -/* - Copyright (C) 2013 Mike Hansen - - This file is part of FLINT. - - FLINT is free software: you can redistribute it and/or modify it under - the terms of the GNU Lesser General Public License (LGPL) as published - by the Free Software Foundation; either version 3 of the License, or - (at your option) any later version. See . -*/ - -#include "fq_zech.h" -#include "fq_zech_mat.h" - -#ifdef T -#undef T -#endif - -#define T fq_zech -#define CAP_T FQ_ZECH -#include "fq_mat_templates/test/t-solve_triu_recursive.c" -#undef CAP_T -#undef T diff --git a/src/gr/fq.c b/src/gr/fq.c index ba9305c592..e8228a0db8 100644 --- a/src/gr/fq.c +++ b/src/gr/fq.c @@ -670,6 +670,24 @@ _gr_fq_mat_mul(fq_mat_t res, const fq_mat_t x, const fq_mat_t y, gr_ctx_t ctx) return GR_SUCCESS; } +int +_gr_fq_mat_nonsingular_solve_tril(fq_mat_t X, const fq_mat_t L, const fq_mat_t B, int unit, gr_ctx_t ctx) +{ + if (B->r < 64 || B->c < 64) + return gr_mat_nonsingular_solve_tril_classical((gr_mat_struct *) X, (const gr_mat_struct *) L, (const gr_mat_struct *) B, unit, ctx); + else + return gr_mat_nonsingular_solve_tril_recursive((gr_mat_struct *) X, (const gr_mat_struct *) L, (const gr_mat_struct *) B, unit, ctx); +} + +int +_gr_fq_mat_nonsingular_solve_triu(fq_mat_t X, const fq_mat_t U, const fq_mat_t B, int unit, gr_ctx_t ctx) +{ + if (B->r < 64 || B->c < 64) + return gr_mat_nonsingular_solve_triu_classical((gr_mat_struct *) X, (const gr_mat_struct *) U, (const gr_mat_struct *) B, unit, ctx); + else + return gr_mat_nonsingular_solve_triu_recursive((gr_mat_struct *) X, (const gr_mat_struct *) U, (const gr_mat_struct *) B, unit, ctx); +} + int _gr_fq_mat_charpoly(fq_struct * res, const fq_mat_t mat, gr_ctx_t ctx) { @@ -776,6 +794,8 @@ gr_method_tab_input _fq_methods_input[] = {GR_METHOD_POLY_MULLOW, (gr_funcptr) _gr_fq_poly_mullow}, {GR_METHOD_POLY_ROOTS, (gr_funcptr) _gr_fq_roots_gr_poly}, {GR_METHOD_MAT_MUL, (gr_funcptr) _gr_fq_mat_mul}, + {GR_METHOD_MAT_NONSINGULAR_SOLVE_TRIL, (gr_funcptr) _gr_fq_mat_nonsingular_solve_tril}, + {GR_METHOD_MAT_NONSINGULAR_SOLVE_TRIU, (gr_funcptr) _gr_fq_mat_nonsingular_solve_triu}, {GR_METHOD_MAT_CHARPOLY, (gr_funcptr) _gr_fq_mat_charpoly}, {GR_METHOD_MAT_REDUCE_ROW, (gr_funcptr) _gr_fq_mat_reduce_row}, {0, (gr_funcptr) NULL}, diff --git a/src/gr/fq_nmod.c b/src/gr/fq_nmod.c index 3c9f488e06..7f3164a32e 100644 --- a/src/gr/fq_nmod.c +++ b/src/gr/fq_nmod.c @@ -633,6 +633,24 @@ _gr_fq_nmod_mat_mul(fq_nmod_mat_t res, const fq_nmod_mat_t x, const fq_nmod_mat_ return GR_SUCCESS; } +int +_gr_fq_nmod_mat_nonsingular_solve_tril(fq_nmod_mat_t X, const fq_nmod_mat_t L, const fq_nmod_mat_t B, int unit, gr_ctx_t ctx) +{ + if (B->r < 64 || B->c < 64) + return gr_mat_nonsingular_solve_tril_classical((gr_mat_struct *) X, (const gr_mat_struct *) L, (const gr_mat_struct *) B, unit, ctx); + else + return gr_mat_nonsingular_solve_tril_recursive((gr_mat_struct *) X, (const gr_mat_struct *) L, (const gr_mat_struct *) B, unit, ctx); +} + +int +_gr_fq_nmod_mat_nonsingular_solve_triu(fq_nmod_mat_t X, const fq_nmod_mat_t U, const fq_nmod_mat_t B, int unit, gr_ctx_t ctx) +{ + if (B->r < 64 || B->c < 64) + return gr_mat_nonsingular_solve_triu_classical((gr_mat_struct *) X, (const gr_mat_struct *) U, (const gr_mat_struct *) B, unit, ctx); + else + return gr_mat_nonsingular_solve_triu_recursive((gr_mat_struct *) X, (const gr_mat_struct *) U, (const gr_mat_struct *) B, unit, ctx); +} + int _gr_fq_nmod_mat_charpoly(fq_nmod_struct * res, const fq_nmod_mat_t mat, gr_ctx_t ctx) { @@ -729,6 +747,8 @@ gr_method_tab_input _fq_nmod_methods_input[] = {GR_METHOD_POLY_ROOTS, (gr_funcptr) _gr_fq_nmod_roots_gr_poly}, {GR_METHOD_MAT_MUL, (gr_funcptr) _gr_fq_nmod_mat_mul}, + {GR_METHOD_MAT_NONSINGULAR_SOLVE_TRIL, (gr_funcptr) _gr_fq_nmod_mat_nonsingular_solve_tril}, + {GR_METHOD_MAT_NONSINGULAR_SOLVE_TRIU, (gr_funcptr) _gr_fq_nmod_mat_nonsingular_solve_triu}, {GR_METHOD_MAT_CHARPOLY, (gr_funcptr) _gr_fq_nmod_mat_charpoly}, {GR_METHOD_MAT_REDUCE_ROW, (gr_funcptr) _gr_fq_nmod_mat_reduce_row}, {0, (gr_funcptr) NULL}, diff --git a/src/gr/fq_zech.c b/src/gr/fq_zech.c index a6d693cfd8..528f2f7bec 100644 --- a/src/gr/fq_zech.c +++ b/src/gr/fq_zech.c @@ -521,6 +521,24 @@ _gr_fq_zech_mat_mul(fq_zech_mat_t res, const fq_zech_mat_t x, const fq_zech_mat_ return GR_SUCCESS; } +int +_gr_fq_zech_mat_nonsingular_solve_tril(fq_zech_mat_t X, const fq_zech_mat_t L, const fq_zech_mat_t B, int unit, gr_ctx_t ctx) +{ + if (B->r < 64 || B->c < 64) + return gr_mat_nonsingular_solve_tril_classical((gr_mat_struct *) X, (const gr_mat_struct *) L, (const gr_mat_struct *) B, unit, ctx); + else + return gr_mat_nonsingular_solve_tril_recursive((gr_mat_struct *) X, (const gr_mat_struct *) L, (const gr_mat_struct *) B, unit, ctx); +} + +int +_gr_fq_zech_mat_nonsingular_solve_triu(fq_zech_mat_t X, const fq_zech_mat_t U, const fq_zech_mat_t B, int unit, gr_ctx_t ctx) +{ + if (B->r < 64 || B->c < 64) + return gr_mat_nonsingular_solve_triu_classical((gr_mat_struct *) X, (const gr_mat_struct *) U, (const gr_mat_struct *) B, unit, ctx); + else + return gr_mat_nonsingular_solve_triu_recursive((gr_mat_struct *) X, (const gr_mat_struct *) U, (const gr_mat_struct *) B, unit, ctx); +} + int _gr_fq_zech_mat_charpoly(fq_zech_struct * res, const fq_zech_mat_t mat, gr_ctx_t ctx) { @@ -616,6 +634,8 @@ gr_method_tab_input _fq_zech_methods_input[] = {GR_METHOD_POLY_ROOTS, (gr_funcptr) _gr_fq_zech_roots_gr_poly}, {GR_METHOD_MAT_MUL, (gr_funcptr) _gr_fq_zech_mat_mul}, + {GR_METHOD_MAT_NONSINGULAR_SOLVE_TRIL, (gr_funcptr) _gr_fq_zech_mat_nonsingular_solve_tril}, + {GR_METHOD_MAT_NONSINGULAR_SOLVE_TRIU, (gr_funcptr) _gr_fq_zech_mat_nonsingular_solve_triu}, {GR_METHOD_MAT_CHARPOLY, (gr_funcptr) _gr_fq_zech_mat_charpoly}, {0, (gr_funcptr) NULL}, };