Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor fq_*_mat solving code using generics #2181

Merged
merged 1 commit into from
Jan 24, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 0 additions & 71 deletions doc/source/fq_mat.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
--------------------------------------------------------------------------------
Expand Down Expand Up @@ -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
Expand All @@ -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
--------------------------------------------------------------------------------
Expand Down
71 changes: 0 additions & 71 deletions doc/source/fq_nmod_mat.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
--------------------------------------------------------------------------------
Expand Down Expand Up @@ -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
Expand All @@ -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
--------------------------------------------------------------------------------
Expand Down
70 changes: 0 additions & 70 deletions doc/source/fq_zech_mat.rst
Original file line number Diff line number Diff line change
Expand Up @@ -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
--------------------------------------------------------------------------------
Expand Down Expand Up @@ -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
Expand All @@ -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
--------------------------------------------------------------------------------
Expand Down
6 changes: 0 additions & 6 deletions src/fq/mat_templates.c
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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"
Expand Down
14 changes: 2 additions & 12 deletions src/fq_mat/test/main.c
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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"
Expand All @@ -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),
Expand All @@ -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),
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
23 changes: 0 additions & 23 deletions src/fq_mat/test/t-lu_recursive.c

This file was deleted.

23 changes: 0 additions & 23 deletions src/fq_mat/test/t-solve_tril_classical.c

This file was deleted.

Loading
Loading