Creating block sparse matrices

Hi, I’m looking to initialize a block-sparse matrix (conceptually sparse - it is fine that it will end up storing the 0s in memory). However, I am not sure which is the most idiomatic way to do that. Of course I could start from primitives again and use from_fn, but that’s awkward. My attempt was as follows:

use nalgebra as na;

fn main() {
    let expected = na::DMatrix::from_row_slice(4, 4, &[
        1., 2., 0., 0., 
        3., 4., 0., 0., 
        0., 0., 1., 2.,
        0., 0., 3., 4.,
    ]);

    let block = na::DMatrix::from_row_slice(2, 2, &[
        1., 2.,
        3., 4.
    ]);

    let mut bsparse = na::DMatrix::zeros(4, 4);

    //bsparse.index_mut((..2, ..2)); // this is not an invalid expression
    //*bsparse.index_mut((..2, ..2)) = block; // <very_complicated_type> cannot be dereferenced
    bsparse.index_mut((..2, ..2)) = block;
    bsparse.index_mut((2.., 2..)) = block;

    println!("{}\n{}", bsparse, expected);

    assert_eq!(bsparse, expected);
}

However it produces two kinds of errors:

bsparse.index_mut((..2, ..2)) = block;
                                ^^^^^ expected struct `nalgebra::base::matrix_slice::SliceStorageMut`, found struct `nalgebra::base::vec_storage::VecStorage`

which suggests that in order to assign to a mutable slice I need to also construct a slice, not just an owned matrix. Normally I would try to dereference the mutable borrow when assigning, but I am being told that the type returned by index_mut does not support dereferencing.

error[E0614]: type `nalgebra::base::matrix::Matrix<_, nalgebra::base::dimension::Dynamic, nalgebra::base::dimension::Dynamic, nalgebra::base::matrix_slice::SliceStorageMut<'_, _, nalgebra::base::dimension::Dynamic, nalgebra::base::dimension::Dynamic, nalgebra::base::dimension::U1, nalgebra::base::dimension::Dynamic>>` cannot be dereferenced

I can resolve the type mismatch error by matching the exact type on both sides

bsparse.index_mut((..2, ..2)) = block.index_mut((.., ..))

which shuts the compiler up in this regard, as little sense as this incantation makes. However, at all times there is another error

bsparse.index_mut((..2, ..2)) = block;
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ left-hand of expression not valid

which is slightly surprising me, since removing the = block; part results in an expression which compiles (although it does nothing).

I would be grateful for both explanation how to best form block-sparse matrices (in numpy I could use stack/concatenate functions for instance, but I couldn’t find those in nalgebra) and how I would go about getting my current approach to work.

Hi!

What you are searching for is the .copy_from method that you can use insead of the assignation:

    bsparse.index_mut((..2, ..2)).copy_from(&block);
    bsparse.index_mut((2.., 2..)).copy_from(&block);

Here is why what you tried did not compile:

This happens because bsparse.index_mut((..2, ..2)) returns a matrix slice which is a value (in particular, it is a struct containing, among other things, a pointer). So what you are attempting to white here is similar to trying to write 10 = 20;. The left-hand side of an assignation can either be a pointer you dereference *variable = 20 or a mutable variable (or a field on a mutable struct). You can’t have a plain value like 10 in this example, or like bsparse.index_mut((..2, ..2)).

Something like the following would compile however:

let mut target = bsparse.index_mut((..2, ..2));
target = block.index_mut((.., ..));

But that won’t do what you want. It will simply make output become the block matrix viewed as a slice. It won’t copy the content of block into your target slice. This is because rust does not allow overloading assignation. That’s why what you need to use is:

let mut output = bsparse.index_mut((..2, ..2));
output.copy_from(&block);

// One-line version:
bsparse.index_mut((..2, ..2)).copy_from(&block);

Finaly, I think what you are trying to do here is the best way of initializing a block-diagonal matrix. Feel free to open an issue on github, suggesting other matrix construction functions (or even open a pull request if you desire to add them to nalgebra yourself)!

Small addition to my previous answer: if your block is small and has a size known at compile-time, you may want to consider using a fixed-size matrix instead of DMatrix:

    let block = na::Matrix2::new(1., 2.,
                                 3., 4.)

That way you will avoid an unnecessary allocation. The .copy_from method will still work with a Matrix2.

Thank you for your detailed answer, this solves my problem.

I think there was a bit of misunderstanding between us - I understand the rules of rust type checking but I was expecting the type resulting from index_mut to implement DerefMut which would enable syntax like *bsparse.index_mut((..2, ..2)) = block;. I think this is a more or less idiomatic way of implementing such “advanced getters”. I could indeed check if it does in the documentation (I admit I did not, due to some brainfart).

Please don’t take it as any offence or nitpicking, but I am having some discoverability issues (vide my previous topic) and I am hoping my detailed explanations can come in handy for someone working on the docs, somewhere in the future. So far I am finding every feature I need, yet the API is structured quite differently than I would expect and the ocean of generic parameters and bounds makes exploring the docs rather difficult. At times I wish I had a rust REPL, so that I could grab the type I have at hand and ask “ok, so what is the API of that particular type?”.

The use of DerefMut is considered idiomatic only for smart-pointers as stated in the std doc. I don’t think there is any way to make *bsparse.index_mut((..2, ..2)) = block; work anyway because there is no Deref::Target type that would be suitable given that the matrix view components are not necessarily contiguous in memory.

No offense taken. Discoverability issue is a known problem with nalgebra because of the generics and because of a bug on rustdoc. We try to make it easier by having a user-guide; and this forum is also a good way for others to find answers to already-asked questions as you said. But it is still a problem. It will likely get much better once rust supports specialization (because this will simplify some trait bounds).

Don’t hesitate to open issues describing what you would expect for the API if something feels odd to you. We may or may not change the API as a result, but in both case this will give us a better idea of what people expect or are used to.

Yeah, that would be great! I think this is kind of what fixing the rustdoc bug I mentioned above would do. That way, whenever you click on a type (be it an alias or not) on the documentation, you will get all the methods available for it.