Additional functionality: Trilinos BlockSparseMatrix

[Jean-Paul Pelteret]

To the deal.II developers,

I've written a couple of additions for consideration for the Trilinos BlockSparseMatrix class. Please find attached the patch, and within the zip files a test program.

I've added a print command, as well as add functions fo BlockSparseMatrix. So far, this functionality includes the addition of a general index value to the block system and a multiple of another BlockSparseMatrix. There are several "add" functions that need to be filled out, but I haven't yet investigated how they work (e.g. in SparseMatrix) so I've left them as TODO's. However, their prototypes need to be there for deal.II to compile (some other function wants to call the add function of this matrix). The one function that is complete has bounds checking, so to make sure that the global indices provided can possibly exist within the global system. It then searches for the correct block that may hold this index and tries to add to the correct local entry in that block. If that entry does not exist, an error is thrown. This appears to work properly in parallel as well. When adding a multiple of another matrix, we just check that the number of block rows and columns is identical and then the seperate blocks check their own sparsity patterns against one another.

I hope that this is not too late to be considered for the next stable release.

Many thanks,

Jean-Paul

[Wolfgang Bangerth]

> I've written a couple of additions for consideration for the Trilinos
> BlockSparseMatrix class. Please find attached the patch, and within the zip
> files a test program.

Nice, thanks!

Can I ask you to make a couple of changes?
1/ The row/col_index_within_bounds() functions are actually equivalent to the
condition
i < this->m()
or
j < this->n()
Rather than duplicating this functionality, please use the conditions above
and remove the new functions.

2/ You also don't really need the local_block_indices() function. To get the
block_index/index_within_block pair, for example for a row index 'i', you can call
this->row_block_indices.global_to_local(i).

3/ If write_ascii is not yet implemented, leave it out for now.

4/ Looking at the add functions you introduce, there are really only two that
you currently implement:
- The one that adds a single value to position (i,j). But this one already
exists in the base class, see block_matrix_base.h around line 2025-2050.
What does your function do that you don't?
- The one that adds a whole block matrix. Could you try to just put this
function into the BlockMatrixBase class where the argument is then of
type BlockMatrixBase<MatrixType> instead of BlockSparseMatrix. This
functionality would seem generic for all possible matrix types.
If it's true that you don't need the first function and that the second can be
put into the base class then you don't need any of the overloads of add() in
the TrilinosWrappers::BlockSparseMatrix class.


As I said, many thanks for posting the patch. I'm sure we can get this in
before the release!

Cheers
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Jean-Paul Pelteret]

Dear Wolfgang,

Thanks very much for the detailed response. Drats, that was a bit silly - I don't know why I didn't pay more attention to the functionality in the base class. Yes, now I see that BlockMatrixBase does contain all but one of the add functions. Thank you for indicating the index checking methodology as well. So I agree that one simply needs to add the single add() and the print function to the base class itself. I've made these additions and attached another patch.

There was already a (conflicting) prototype for an add (factor, matrix) function on line 778, but it had no definition. I simply commented out the prototype as it suggested another implementation of the the function and I wasn't sure whether there were plans to implement it in some more robust manner.

Best regards,

Jean-Paul

[Timo Heister]

Hi Jean-Paul,

I am applying and testing your patch. Thanks.

How did you generate that .patch file, btw.? I had to apply it by hand
because "svn patch" and "patch" would not eat it.

To generate a patch: svn diff >my_patch.patch
(you can also only pick certain files by listing them)
To apply a patch: patch -p0 <my_patch.patch
(or with a new svn): svn patch my_patch.patch

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>


--
Timo Heister
http://www.math.tamu.edu/~heister/

[Jean-Paul Pelteret]

Hi Timo,

Apologies for that - I just used diff itself. Would you like me to remake it? Sorry for the novice mistake, I don't use svn for any purpose other than deal.II.

Kinds regards,

Jean-Paul

[Jean-Paul Pelteret]

Hi Timo,

If you haven't already done it, here's the amended version. You learn something new every day...

Kind regards,

J-P