[petsc-users] Dropping single entries in matrix insertion for block matrices

[Smith, Barry F.]

   Lawrence,

      I understand what you want and it is a reasonable request. The problem is that currently ISLocalToGlobalMappingCreate() when used with block vectors and matrices is always based on blocks, that is, from the manual page, "There is one integer value in indices per block and it represents the actual indices bs*idx + j, where j=0,..,bs-1". If you look at ISLocalToGlobalMappingApply() it uses the code:

 out[i] = idx[in[i]/bs]*bs + (in[i] % bs);

to do the mapping; that is the idx[] mapping is stored only by block, not by point.

  But I think you may be able to get the effect you want by "managing the local to global mapping yourself before calling MatSetValues()". That is, you create a ISLocalToGlobal object yourself, not based on blocks, then when setting your Dirichlet conditions you call ISLocalToGlobalMapping() apply yourself and then call MatSetValues() using the resulting indices.


  Barry



> On Dec 9, 2018, at 9:55 AM, Lawrence Mitchell via petsc-users <petsc-users at mcs.anl.gov> wrote:
> 
> Dear petsc-users,
> 
> I have a matrix with a block size > 1. I would sometimes like to insert into it, applying Dirichlet conditions to one component of the block. Now, for normal block (or when the block size is 1) dirichlet conditions, I do this by swapping out the local to global map for one with negative entries. That is, I do:
> 
> MatCreate(..., &mat);
> 
> MatSetSizes(mat, ...);
> 
> MatSetUp(mat);
> 
> MatSetLocalToGlobalMapping(mat, rmap, cmap);
> // normal insertion
> MatSetValuesLocal(mat, ...);
> // or MatSetValuesBlockedLocal(mat, ...);
> 
> When I want to drop some entries I do:
> 
> ISLocalToGlobalMappingCreate(..., &rmap_with_bcs);
> 
> MatSetLocalToGlobalMapping(mat, rmap_with_bcs, cmap_with_bcs);
> 
> MatSetValuesLocal(mat, ...);
> 
> And those entries for which rmap_with_bcs maps to negative rows get dropped.
> 
> The problem arises when I want to do this for matrices with blocksize > 1.
> 
> Now the ISLocalToGlobalMapping must have the same block size as the matrix. But I can't then use it to drop single components.
> 
> In petsc4py speak I would like to be able to do:
> 
> from petsc4py import PETSc
> mat = PETSc.Mat().create()
> mat.setSizes((3, 3))
> mat.setBlockSize(3)
> mat.setUp()
> lgmap = PETSc.LGMap().create([0, 1, -1], bsize=1)
> mat.setLGMap(lgmap, lgmap)
> 
> This results in an error:
> 
> [0] MatSetLocalToGlobalMapping() line 2009 in matrix.c
> [0] PetscLayoutSetISLocalToGlobalMapping() line 247 in pmap.c
> [0] Petsc has generated inconsistent data
> [0] Blocksize of layout 3 must match that of mapping 1
> 
> I suppose this is to catch me from shooting myself in the foot when subsequently doing MatSetValuesBlockedLocal? But I can arrange that this doesn't occur.
> 
> Thoughts?
> 
> If I could do this, it would make my code a lot simpler, but I am aware it might make error checking more difficult.
> 
> Lawrence