Thanks Jed, this is great. I was trying to implement what was proposed in http://lists.mcs.anl.gov/pipermail/petsc-users/2011-September/009991.html, by you as I just realised. For this the matrices would be of varying size, so for three local systems, I would need:
S \otimes Id \otimes Id ; Id \otimes S \otimes Id ; Id \otimes Id \otimes S. Realistic systems are of the order of at least 20 local systems, which would mean that I need the matrices of the above structure for all 20 possible positions. Again, I really appreciate the help. Thanks, Mathis On Tue, Apr 1, 2014 at 4:39 PM, Jed Brown <[email protected]> wrote: > Mathis Friesdorf <[email protected]> writes: > > > Thanks Jed! The libaries you are pointing to look very interesting > indeed. > > For the particular implementation I have in mind, I was hoping to get > away > > with something easier, as my time to work on this is a bit limited. All I > > really need is a way to construct a block-diagonal matrix where each > block > > is given by the same matrix and all blocks are on the diagonal. Am I > wrong > > to assume that this should be possible to implement with TAIJ? > > That is exactly TAIJ. How large are the diagonal blocks? The > implementation is intended to be fast only for small sizes, though I > could make it good for large sizes as well. In the special case of > > (I \otimes S) X > > where I is the sparse identity and S is the dense diagonal block, this > involves interpreting X as a matrix rather than a vector and using dgemm. >
