Merge pull request #2181 from fredrik-johansson/matrix5

Refactor fq_*_mat solving code using generics

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
 --------------------------------------------------------------------------------

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
 --------------------------------------------------------------------------------

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
 --------------------------------------------------------------------------------

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"

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),

src/fq_mat/test/t-lu_classical.c → 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