sparsemat / sprs

sparse linear algebra library for rust
Apache License 2.0
386 stars 45 forks source link

fixed a bug in mul_acc_mat_vec_csr and mul_acc_mat_vec_csc #218

Closed cxzheng closed 3 years ago

cxzheng commented 3 years ago

There is a bug in mul_acc_mat_vec_csr and mul_acc_mat_vec_csc. Basically, the output res_vec needs to be set to zero before accumulating values in the matrix vector multiplication.

I also btw did some small tweaks following the suggestions from cargo clippy

vbarrielle commented 3 years ago

Hello,

That's not a bug, that's the expected behaviour. The function accumulates in the output vector, so should preserve the contents it has on input. It probably should be stated more explicitly in the docs. Maybe some line like This function implements the operation y += A * x.

This accumulation behavior is desirable because it enables adding matrix/vector products with only one storage vector, ie the equivalent of y = A * x + B * x is something like

let mut y = [0; n];
mul_acc_mat_vec_csr(b.view(), x, &mut y);
mul_acc_mat_vec_csr(a.view(), x, &mut y);

One could argue this is not convenient, but this function is quite low-level, the high level interface is the *operator. The drawback of this high-level operator is that it allocates, and does not work on slices. Maybe we could add a version of mul_acc_mat_vec_cs{r,c} that has this zeroing behavior for that purpose, eg mul_mat_vec_cs{r,c}.

The clippy suggested fixes look nice.

cxzheng commented 3 years ago

Hello,

That's not a bug, that's the expected behaviour. The function accumulates in the output vector, so should preserve the contents it has on input. It probably should be stated more explicitly in the docs. Maybe some line like This function implements the operation y += A * x.

This accumulation behavior is desirable because it enables adding matrix/vector products with only one storage vector, ie the equivalent of y = A * x + B * x is something like

let mut y = [0; n];
mul_acc_mat_vec_csr(b.view(), x, &mut y);
mul_acc_mat_vec_csr(a.view(), x, &mut y);

One could argue this is not convenient, but this function is quite low-level, the high level interface is the *operator. The drawback of this high-level operator is that it allocates, and does not work on slices. Maybe we could add a version of mul_acc_mat_vec_cs{r,c} that has this zeroing behavior for that purpose, eg mul_mat_vec_cs{r,c}.

The clippy suggested fixes look nice.

Sorry for my misunderstanding of the two functions' behaviors. What you explained makes sense. I changed them back.

cxzheng commented 3 years ago

Thanks for your quick iteration @cxzheng. Do you want me to add the non-accumulation versions mul_acc_mat_vec_cs{r,c}, or are you using the * operator?

It's okay for me. I'm just gonna set the vector to zero before calling mul_acc_mat_vec_cs{r,c}. I think it makes sense to stay with the same style like saxpy in BLAS in this case.