[petsc-users] Retrieving a SubMatrix from a Matrix in parallel

[Barry Smith]

On Jun 17, 2014, at 5:03 PM, Arthur Kurlej <akurlej at gmail.com> wrote:

> Hi all,
> 
> I wish to extract a parallel AIJ submatrix (B) from another parallel AIJ matrix (A), in particular, if I have an NxN matrix, I wish to extract the first N-1 rows and N-1 columns of A.
> 
> I can do this perfectly fine when running my code sequentially using ISCreateStride and MatGetSubMatrix like so:
> ISCreateStride(PETSC_COMM_WORLD,N,1,1,&is);

   ISCreateStride() does not work this way. Each process should provide only the part of the stride space of rows that it wants. What you indicate above has each process taking the stride space from 1 to N.   

   So for your parallel matrix below you could have process 0 use in ISCreateStride(  , 2,0,1,&isrow) while process 1 has ( ,1,2,1,&isrow) so you are saying the first process wants rows 0 and 1 while the second process wants rows 2. Meanwhile both processes want columns 0,1,2 so they would each create an IS with ISCreateStride(  , 3,0,1,&iscol)

  Barry

> MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&B);
> 
> And it behaves as I would expect, but when running in parallel, the B matrix increases past it's allocated size of N-1xN-1. I've pasted the example of a simple example output below:
> 
> 
> ORIGINAL SEQUENTIAL MATRIX
> row 0: (0, 1)  (1, 2)  (2, 3)  (3, 4)
> row 1: (0, 1)  (1, 2)  (2, 3)  (3, 4)
> row 2: (0, 1)  (1, 2)  (2, 3)  (3, 4)
> row 3: (0, 1)  (1, 2)  (2, 3)  (3, 4)
> SEQUENTIAL SUBMATRIX
> row 0: (0, 1)  (1, 2)  (2, 3)
> row 1: (0, 1)  (1, 2)  (2, 3)
> row 2: (0, 1)  (1, 2)  (2, 3)
> 
> 
> ORIGINAL PARALLEL MATRIX
> row 0: (0, 1)  (1, 2)  (2, 3)  (3, 4)
> row 1: (0, 1)  (1, 2)  (2, 3)  (3, 4)
> row 2: (0, 1)  (1, 2)  (2, 3)  (3, 4)
> row 3: (0, 1)  (1, 2)  (2, 3)  (3, 4)
> PARALLEL "SUBMATRIX"
> row 0: (3, 1)  (4, 2)  (5, 3)
> row 1: (3, 1)  (4, 2)  (5, 3)
> row 2: (3, 1)  (4, 2)  (5, 3)
> row 3: (3, 1)  (4, 2)  (5, 3)
> row 4: (3, 1)  (4, 2)  (5, 3)
> row 5: (3, 1)  (4, 2)  (5, 3)
> 
> 
> I'm not really sure what is going on with the index set, and I would really appreciate some help.
> 
> 
> Thanks,
> Arthur
>