Sparse3

[kartikv]

trafficstars

Just the start, but converting the sparse2 stuff to the FLINT 3 generic ring system.

[fredrik-johansson]

One subtle thing is to handle "non-canonical" rings, i.e. rings with inexact zeros (like arb or gr_series) or undecidable or non-canonically represented zeros (like ca).

For example, when checking equality of two matrices or vectors, we cannot generally first compare the indices, then compare the entries. If A[i] exists but B[i] doesn't, we have to check if B[i] is possibly zero.

The good news is that we can check whether this is an issue with gr_ctx_is_canonical.

So we can do something like:

if (gr_ctx_is_canonical(ctx) == T_TRUE)
{
    ... fast algorithm comparing indices, then entries
}
else
{
    ... slow algorithm
}

Some functions (solvers?) are probably going to be a bit messy to implement for non-canonical rings, so in those functions you could just put in

if (gr_ctx_is_canonical(ctx) != T_TRUE)
    return GR_UNABLE;

and do the canonical case, leaving a non-canonical version for later.

[fredrik-johansson]

BTW, there is a use case for sparse vectors with bignum indices, represented by fmpz rather than ulong. In particular, some mpoly code could use such vectors. However, I think this is fringe enough that it can be done separately when the need arises; no need to complicate things with bignums here.

At a first glance, the data structures here look good.

[kartikv]

I think I did it the proper generic way now: explicitly check for all zeroes, ones, and equality, and return

• T_FALSE if any check returns T_FALSE; otherwise
• T_UNKNOWN if any check returns T_UNKNOWN; otherwise
• T_TRUE I'll start on testing tomorrow, I'll start with a ring I understand to deal with the standard bugs, but would be good to have a suite of examples of a noncanonical ring to see what is broken.

[fredrik-johansson]

This looks great so far. Obviously, the blocker for merging is that there are tests that crash and that the documentation doesn't build cleanly.

A few nits to pick:

• gr_*_mat_mul_mat should just be gr_*_mat_mul and gr_*_mat_mul should be named gr_lil_mat_mul_entrywise
• "Addition and subtraction are not provided because they would create dense output." Well no, matrix + scalar should just fill in the main diagonal.
• gr_*_mat_mul_vec should maybe be named _gr_*_mat_mul_vec with a corresponding gr_*_mat_mul_vec taking a gr_vec_t as input rather than a raw pointer to elements. The documentation doesn't match the code here. There could be such dual versions for other functions as well, e.g. solvers.
• gr_*_mat_nz_product: "set res to the product of the nonzero entries in mat". This should return GR_UNABLE if encountering any entry that is maybe zero. I suggest adding helper functions to check if sparse vectors/matrices are in a strong canonical form in the sense that there are no maybe-zero entries.
• In the solvers: if (M->r != M->c) return 0; /* TODO */ should return GR_UNABLE (if you plan to support nonsquare systems some time in the future) or GR_DOMAIN (if this isn't allowed at all). You should also clarify in the documentation whether the matrices need to be full-rank and what happens otherwise.
• gr_coo_mat_is_zero should not mutate its argument; you could have a separate _gr_coo_mat_is_zero that does.
• There are some places in the code where int is used as an index (should be slong).

BTW, a great way to test the sparse vectors and matrices is to wrap them for generics like gr_vec_t and gr_mat_t (see gr/vector.c and gr/matrix.c). You can then test a wide range of operations with gr_test_ring.

[kartikv]

This looks great so far. Obviously, the blocker for merging is that there are tests that crash and that the documentation doesn't build cleanly.

I may need some help with this, since the tests do not crash for me (I would not push otherwise)...are there settings I need to try when I configure/build flint to see these issues? I can try to get some other platforms to test on, I'm on an M3 Mac.

A few nits to pick:

• gr_*_mat_mul_mat should just be gr_*_mat_mul and gr_*_mat_mul should be named gr_lil_mat_mul_entrywise

Happy with the latter, but the former was intended to convey that I'm multiplying by a gr_mat, not another gr_lil_mat: my assumption is that that with no other descriptor, one would assume that the two arguments would be the same type.

• "Addition and subtraction are not provided because they would create dense output." Well no, matrix + scalar should just fill in the main diagonal.

Sorry, that was copied over from a previous version, yes that will change.

• gr_*_mat_mul_vec should maybe be named _gr_*_mat_mul_vec with a corresponding gr_*_mat_mul_vec taking a gr_vec_t as input rather than a raw pointer to elements. The documentation doesn't match the code here. There could be such dual versions for other functions as well, e.g. solvers.

Sounds good, will do.

• gr_*_mat_nz_product: "set res to the product of the nonzero entries in mat". This should return GR_UNABLE if encountering any entry that is maybe zero. I suggest adding helper functions to check if sparse vectors/matrices are in a strong canonical form in the sense that there are no maybe-zero entries.

That's fine.

• In the solvers: if (M->r != M->c) return 0; /* TODO */ should return GR_UNABLE (if you plan to support nonsquare systems some time in the future) or GR_DOMAIN (if this isn't allowed at all). You should also clarify in the documentation whether the matrices need to be full-rank and what happens otherwise.

Sorry, those should have all been replaced with is_square, will take care of that. And nothing needs to be full rank, I'll clarify in the documentation but Wiedemann will just return a solution and Lanzcos a pseudo-solution (if one exists).

• gr_coo_mat_is_zero should not mutate its argument; you could have a separate _gr_coo_mat_is_zero that does.

Ok, is it alright if it just returns GR_UNABLE if the matrix is not in canonical form?

• There are some places in the code where int is used as an index (should be slong).

Thanks, I'll dig out all of those.

BTW, a great way to test the sparse vectors and matrices is to wrap them for generics like gr_vec_t and gr_mat_t (see gr/vector.c and gr/matrix.c). You can then test a wide range of operations with gr_test_ring.

Will do, once things are in a good shape for the rings they definitely should work for. :-)

[kartikv]

"without having probably no solution"? "Without probably having no solution"? Happy for other suggestions.

On Thu, Apr 4, 2024, 7:42 AM Dima Pasechnik @.***> wrote:

@.**** commented on this pull request.

In doc/source/gr_sparse_mat.rst https://github.com/flintlib/flint/pull/1845#discussion_r1551820695:

+.. function:: int gr_csr_mat_mul_mat(gr_mat_t C, const gr_csr_mat_t A, const gr_mat_t B, gr_ctx_t ctx)

•          int gr_lil_mat_mul_mat(gr_mat_t C, const gr_lil_mat_t A, const gr_mat_t B, gr_ctx_t ctx)
  

• Set C equal to A \cdot B, i.e., perform right multiplication by B.

+Solving, nullvector, and nullspace computation +-------------------------------------------------------------------------------- + +.. function:: int gr_lil_mat_solve_lanczos(gr_ptr x, const gr_lil_mat_t M, gr_srcptr b, flint_rand_t state, gr_ctx_t ctx)

•          int gr_lil_mat_solve_block_lanczos(gr_ptr x, const gr_lil_mat_t M, gr_srcptr b, slong block_size, flint_rand_t state, gr_ctx_t ctx)
  

•          int gr_lil_mat_solve_wiedemann(gr_ptr x, const gr_lil_mat_t M, gr_srcptr b, gr_ctx_t ctx)
  

•          int gr_lil_mat_solve_block_wiedemann(gr_ptr x, const gr_lil_mat_t M, gr_srcptr b, slong block_size, flint_rand_t state, gr_ctx_t ctx)
  

• Solve Mx = b for a sparse matrix M and M->c long column vector b, using either Lanczos' or Wiedemann's algorithm.
• Both are randomized algorithms which use a given random state machine, and thus may fail without provably no solution

"fail without ... no solution" is very hard to read, sort of triple negation

— Reply to this email directly, view it on GitHub https://github.com/flintlib/flint/pull/1845#pullrequestreview-1980138488, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAKDUIIBMP53TNWV6T4E3ZLY3VRDXAVCNFSM6AAAAABE3JHKI6VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMYTSOBQGEZTQNBYHA . You are receiving this because you authored the thread.Message ID: @.***>

[dimpase]

Let's be more precise here. Randomized algorithms are of two types - Las Vegas and Monte Carlo. https://en.wikipedia.org/wiki/Randomized_algorithm

What type do we have here?

[kartikv]

Neither, at least in general:

For a floating point represented field (which I don't think we have?), I think it's essentially a Monte Carlo algorithm, in that it will produce an epsilon approximation of the correct answer.

For exact (eg finite) fields, it will give the correct answer or just fail, mostly independent of the random input it uses and dependent instead on whether the matrix is well conditioned...the randomness is mostly so that, for things like kernel finding, it will eventually get the whole null space. But there are no guarantees whatsoever.

On Thu, Apr 4, 2024, 11:11 AM Dima Pasechnik @.***> wrote:

Let's be more precise here. Randomized algorithms are of two types - Las Vegas and Monte Carlo. https://en.wikipedia.org/wiki/Randomized_algorithm

What type do we have here?

— Reply to this email directly, view it on GitHub https://github.com/flintlib/flint/pull/1845#issuecomment-2037874364, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAKDUILFSDE4CL5F6D3Y2JLY3WJURAVCNFSM6AAAAABE3JHKI6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAMZXHA3TIMZWGQ . You are receiving this because you authored the thread.Message ID: @.***>

[dimpase]

an $\epsilon$-approximation is a correct answer in the case of ball arithmetic (with balls of this radius).

[munuxi]

In the macro GR_SPV_RFL_TEMPLATE of gr_sparse_vec/arith.c, if a_nnz or b_nnz is zero, the inital a_nz_idx or a_nz_idx is -1 and it's impossible to visit (A_VEC)->inds[a_nz_idx] (or B) in this case:

    ...
    a_nz_idx = a_nnz-1;                                                                                    \
    b_nz_idx = b_nnz-1;                                                                                    \
    dest_nz_idx = new_nnz-1;                                                                               \
    while (a_nz_idx >= -1 && b_nz_idx >= -1 && status == GR_SUCCESS)                                       \
    {                                                                                                      \
        if (a_nz_idx == -1 && b_nz_idx == -1) break;                                                       \
        a_ind = (A_VEC)->inds[a_nz_idx];                                                                   \
        b_ind = (B_VEC)->inds[b_nz_idx];                                                                   \
        if (b_nz_idx == -1 || (a_nz_idx >= 0 && a_ind > b_ind))                                                                              \
        ...

This could appear when a row of a sparse matrix is empty, which leads to the crash of some tests of sparse matrix.

[albinahlback]

Is there any I/O-functionality?

[fredrik-johansson]

In the macro GR_SPV_RFL_TEMPLATE of gr_sparse_vec/arith.c, if a_nnz or b_nnz is zero, the inital a_nz_idx or a_nz_idx is -1 and it's impossible to visit (A_VEC)->inds[a_nz_idx] (or B) in this case:

    ...
    a_nz_idx = a_nnz-1;                                                                                    \
    b_nz_idx = b_nnz-1;                                                                                    \
    dest_nz_idx = new_nnz-1;                                                                               \
    while (a_nz_idx >= -1 && b_nz_idx >= -1 && status == GR_SUCCESS)                                       \
    {                                                                                                      \
        if (a_nz_idx == -1 && b_nz_idx == -1) break;                                                       \
        a_ind = (A_VEC)->inds[a_nz_idx];                                                                   \
        b_ind = (B_VEC)->inds[b_nz_idx];                                                                   \
        if (b_nz_idx == -1 || (a_nz_idx >= 0 && a_ind > b_ind))                                                                              \
        ...

This could appear when a row of a sparse matrix is empty, which leads to the crash of some tests of sparse matrix.

Good catch. I pushed a fix and the sparse_mat tests now pass as follows:

gr_sparse_mat_init...
gr_sparse_mat_init                                0.00   (PASS)
gr_sparse_mat_conversion...
gr_sparse_mat_conversion                          0.22   (PASS)
gr_sparse_mat_randtest...
gr_sparse_mat_randtest                            0.11   (PASS)
gr_sparse_mat_arith...
gr_sparse_mat_arith                               1.79   (PASS)
gr_sparse_mat_mul...
gr_sparse_mat_mul                                 0.06   (PASS)
gr_sparse_mat_solve...
solving Ax = b........................PASS
Solved 0 with rref
Solved 0 with lu
Solved 100 with Lanczos
Solved 0 with block Lanczos
	Found no solution for 100/100 examples
Solved 100 with Wiedemann
Solved 100 with block Wiedemann
gr_sparse_mat_solve                               2.73   (PASS)

All tests passed for gr_sparse_mat.

Those "Solved 0" look suspect?

[kartikv]

Yes, I'm going to finish and check in the rref and lu stuff this weekend, and will figure out what is going on with block lanzos. Thanks for catching the zero length bug!

On Sat, Apr 27, 2024, 5:33 AM Fredrik Johansson @.***> wrote:

In the macro GR_SPV_RFL_TEMPLATE of gr_sparse_vec/arith.c, if a_nnz or b_nnz is zero, the inital a_nz_idx or a_nz_idx is -1 and it's impossible to visit (A_VEC)->inds[a_nz_idx] (or B) in this case:

...
a_nz_idx = a_nnz-1;                                                                                    \
b_nz_idx = b_nnz-1;                                                                                    \
dest_nz_idx = new_nnz-1;                                                                               \
while (a_nz_idx >= -1 && b_nz_idx >= -1 && status == GR_SUCCESS)                                       \
{                                                                                                      \
    if (a_nz_idx == -1 && b_nz_idx == -1) break;                                                       \
    a_ind = (A_VEC)->inds[a_nz_idx];                                                                   \
    b_ind = (B_VEC)->inds[b_nz_idx];                                                                   \
    if (b_nz_idx == -1 || (a_nz_idx >= 0 && a_ind > b_ind))                                                                              \
    ...

This could appear when a row of a sparse matrix is empty, which leads to the crash of some tests of sparse matrix.

Good catch. I pushed a fix and the sparse_mat tests now pass as follows:

gr_sparse_mat_init... gr_sparse_mat_init 0.00 (PASS) gr_sparse_mat_conversion... gr_sparse_mat_conversion 0.22 (PASS) gr_sparse_mat_randtest... gr_sparse_mat_randtest 0.11 (PASS) gr_sparse_mat_arith... gr_sparse_mat_arith 1.79 (PASS) gr_sparse_mat_mul... gr_sparse_mat_mul 0.06 (PASS) gr_sparse_mat_solve... solving Ax = b........................PASS Solved 0 with rref Solved 0 with lu Solved 100 with Lanczos Solved 0 with block Lanczos Found no solution for 100/100 examples Solved 100 with Wiedemann Solved 100 with block Wiedemann gr_sparse_mat_solve 2.73 (PASS)

All tests passed for gr_sparse_mat.

Those "Solved 0" look suspect?

— Reply to this email directly, view it on GitHub https://github.com/flintlib/flint/pull/1845#issuecomment-2080549019, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAKDUIPI6ALPI6FYSNDLQZ3Y7OLKDAVCNFSM6AAAAABE3JHKI6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAOBQGU2DSMBRHE . You are receiving this because you authored the thread.Message ID: @.***>

[kartikv]

Only write and print functions, nothing binary so far: happy to take proposals, I don't know what gr does.

On Sat, Apr 27, 2024, 5:19 AM Albin Ahlbäck @.***> wrote:

Is there any I/O-functionality?

— Reply to this email directly, view it on GitHub https://github.com/flintlib/flint/pull/1845#issuecomment-2080520135, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAKDUINVCKIOXOEH24DEZUDY7OJUVAVCNFSM6AAAAABE3JHKI6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAOBQGUZDAMJTGU . You are receiving this because you authored the thread.Message ID: @.***>

[kartikv]

Only while I'm still debugging/adding requested functionality: they will be removed before the merge.

On Sat, Apr 27, 2024, 5:17 AM Albin Ahlbäck @.***> wrote:

@.**** commented on this pull request.

In src/gr_sparse_mat/test/t-conversion.c https://github.com/flintlib/flint/pull/1845#discussion_r1581821149:

•    // flint_printf("\n\ncoo_mat = "); status |= gr_coo_mat_print_nz(coo_mat, ctx); flint_printf("\nnnz = %d\n", coo_mat->nnz);
  

•    status |= gr_lil_mat_set_coo_mat(lil_mat, coo_mat, ctx);
  

•    // flint_printf("\n\nlil_mat = "); status |= gr_lil_mat_print_nz(lil_mat, ctx); flint_printf("\nnnz = %d\n", lil_mat->nnz);
  

•    status |= gr_lil_mat_set(lil_mat2, lil_mat, ctx);
  

•    // flint_printf("\n\nlil_mat = "); status |= gr_lil_mat_print_nz(lil_mat2, ctx); flint_printf("\nnnz = %d\n", lil_mat2->nnz);
  

•    status |= gr_csr_mat_set_lil_mat(csr_mat, lil_mat2, ctx);
  

•    // flint_printf("\n\ncsr_mat = "); status |= gr_csr_mat_print_nz(csr_mat, ctx); flint_printf("\nnnz = %d\n", csr_mat->nnz);
  

•    status |= gr_mat_set_csr_mat(dmat, csr_mat, ctx);
  

•    // flint_printf("\n\nmat = "); status |= gr_mat_print(dmat, ctx);
  

•    status |= gr_lil_mat_set_mat(lil_mat2, dmat, ctx);
  

•    // flint_printf("\n\nlil_mat = "); status |= gr_lil_mat_print_nz(lil_mat2, ctx); flint_printf("\nnnz = %d\n", lil_mat2->nnz);
  

•    status |= gr_coo_mat_set_lil_mat(coo_mat2, lil_mat2, ctx);
  

•    // flint_printf("\n\ncoo_mat = "); status |= gr_coo_mat_print_nz(coo_mat2, ctx); flint_printf("\nnnz = %d\n", coo_mat2->nnz);
  

Is it beneficial to have theses print statement left in the code?

— Reply to this email directly, view it on GitHub https://github.com/flintlib/flint/pull/1845#pullrequestreview-2026640580, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAKDUIIUB3TPV4643BDVA5LY7OJOTAVCNFSM6AAAAABE3JHKI6VHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMZDAMRWGY2DANJYGA . You are receiving this because you authored the thread.Message ID: @.***>

[albinahlback]

Only write and print functions, nothing binary so far: happy to take proposals, I don't know what gr does.

Ah, I see now. What's the current print format?

[kartikv]

Primarily comma-separated list of column: value for each row, since sparse row form is the core implementation, though I added a coo presentation too.

On Sat, Apr 27, 2024, 7:32 AM Albin Ahlbäck @.***> wrote:

Only write and print functions, nothing binary so far: happy to take proposals, I don't know what gr does.

Ah, I see now. What's the current print format?

— Reply to this email directly, view it on GitHub https://github.com/flintlib/flint/pull/1845#issuecomment-2080792186, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAKDUILYJO66PEKQB3Z4OBTY7OZGTAVCNFSM6AAAAABE3JHKI6VHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDAOBQG44TEMJYGY . You are receiving this because you authored the thread.Message ID: @.***>