On Thu, Jul 21, 2022 at 6:28 AM Emile Soutter <emile.sout...@corintis.com> wrote:
> Dear all, > > I am struggling with the simple following problem : Having a first matrix > B1 of size n1xm1, a second matrix B2 of size n2 x m2, build a matrix M of > size (n1+n2)x(m1+m2) where the blocks B1 and B2 are the "diagonal" of M > (M[0:n1,0:m1]=B1, M[n1:(n1+n2),m1:(m1+m2)]=B2). In my case, the blocks B1 > and B2 are obtained from another routine, directly in the petsc matrix form > (or pyop2.Sparsity form). However the blocks are not squared (n1,n2,m1,m2 > are all different integers). The operation is easy to do with the SetValues > option. However, it takes a large amount of time (too much) when the system > becomes large. I struggle to do it efficiently and in parallel. What method > do you recommend to use to do this as fast as possible? > > Thanks you for any tips, > I think it depends on what you want to do with the final matrix. If you only want MatMult, then I think you can just use MatNest https://petsc.org/main/docs/manualpages/Mat/MatCreateNest/ which will wrap up the submatrices. However, if you want to manipulate the values (factorization, relaxation, etc) then you need to assemble a monolithic matrix. For this you could create the global matrix, and then use https://petsc.org/main/docs/manualpages/Mat/MatCreateLocalRef/ to get a submatrix to assemble directly into, which you pass to your assembly routine. Clearly this is more complicated, but sometimes necessary. Thanks, Matt > Emile > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>