[petsc-users] Confusion with MatGetLocalSubMatrix

[Lawrence Mitchell]

Hi folks,

On 19 Dec 2014, at 18:29, Lawrence Mitchell <lawrence.mitchell at imperial.ac.uk> wrote:

> 
> On 19 Dec 2014, at 18:19, Lawrence Mitchell <lawrence.mitchell at imperial.ac.uk> wrote:
> 
>> ...
>> 
>>>> Calling MatGetLocalSubMatrix results in an error:
>>>> 
>>>> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
>>>> [0]PETSC ERROR: Petsc has generated inconsistent data
>>>> [0]PETSC ERROR: Blocksize of localtoglobalmapping 1 must match that of layout 2
>>> 
>>> Hmm, I'm concerned that this might not work right after some recent
>>> changes to ISLocalToGlobalMapping.  Can you reproduce this with some
>>> code I can run?
> ...
> 
> Here's one that doesn't try to destroy sections that don't exist, and endeavours to print block sizes of all the sub matrices.  I expect to see (2, 2), (2, 1), (1, 2) and (1, 1) I think.

Any suggestions as to where to start looking to fix the problem this code demonstrates?  I'd like to be able to use MatGetLocalSubMatrix with ISes with block sizes > 1.

Cheers,

Lawrence

> Cheers,
> 
> Lawrence
> 
> #include <petsc.h>
> 
> int main(int argc, char **argv)
> {
>    PetscErrorCode ierr;
>    DM v;
>    DM p;
>    DM pack;
>    PetscInt rbs, cbs;
>    PetscInt srbs, scbs;
>    PetscInt i, j;
>    Mat mat;
>    Mat submat;
>    IS *ises;
> 
>    MPI_Comm c;
>    PetscInitialize(&argc, &argv, NULL, NULL);
> 
>    c = PETSC_COMM_WORLD;
>    ierr = DMDACreate1d(c, DM_BOUNDARY_NONE, 10, 2, 1, NULL, &v); CHKERRQ(ierr);
>    ierr = DMDACreate1d(c, DM_BOUNDARY_NONE, 10, 1, 1, NULL, &p); CHKERRQ(ierr);
>    ierr = DMSetFromOptions(v); CHKERRQ(ierr);
>    ierr = DMSetFromOptions(p); CHKERRQ(ierr);
> 
> 
>    ierr = DMCreateMatrix(v, &mat); CHKERRQ(ierr);
> 
>    ierr = MatGetBlockSizes(mat, &rbs, &cbs); CHKERRQ(ierr);
> 
>    ierr = PetscPrintf(c, "Global Mat block size (%d, %d)\n", rbs, cbs); CHKERRQ(ierr);
>    ierr = MatDestroy(&mat); CHKERRQ(ierr);
> 
>    ierr = DMCompositeCreate(c, &pack); CHKERRQ(ierr);
> 
>    ierr = DMCompositeAddDM(pack, v); CHKERRQ(ierr);
>    ierr = DMCompositeAddDM(pack, p); CHKERRQ(ierr);
>    ierr = DMSetFromOptions(pack); CHKERRQ(ierr);
> 
>    ierr = DMCompositeGetLocalISs(pack, &ises); CHKERRQ(ierr);
> 
>    ierr = DMCreateMatrix(pack, &mat); CHKERRQ(ierr);
> 
>    ierr = MatGetBlockSizes(mat, &rbs, &cbs); CHKERRQ(ierr);
> 
>    ierr = PetscPrintf(c, "Global Mat block size (%d, %d)\n", rbs, cbs); CHKERRQ(ierr);
> 
>    for (i=0; i < 2; i++ ) {
>        for (j=0; j < 2; j++ ) {
>            ierr = MatGetLocalSubMatrix(mat, ises[i], ises[j], &submat); CHKERRQ(ierr);
>            ierr = MatGetBlockSizes(submat, &srbs, &scbs); CHKERRQ(ierr);
>            ierr = PetscPrintf(c, "Local (%d, %d) block has block size (%d, %d)\n", i, j, srbs, scbs); CHKERRQ(ierr);
>            ierr = MatDestroy(&submat); CHKERRQ(ierr);
>        }
>    }
>    ierr = MatDestroy(&mat); CHKERRQ(ierr);
> 
>    ierr = DMDestroy(&pack); CHKERRQ(ierr);
>    ierr = DMDestroy(&v); CHKERRQ(ierr);
>    ierr = DMDestroy(&p); CHKERRQ(ierr);
> 
>    ierr = ISDestroy(&(ises[0])); CHKERRQ(ierr);
>    ierr = ISDestroy(&(ises[1])); CHKERRQ(ierr);
>    ierr = PetscFree(ises); CHKERRQ(ierr);
> 
>    PetscFinalize();
> 
>    return 0;
> }

-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 455 bytes
Desc: Message signed with OpenPGP using GPGMail
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150112/e3a25e9e/attachment.pgp>