# [petsc-users] How to extract array portions

FRANCAVILLA MATTEO ALESSANDRO d019598 at polito.it
Wed Apr 6 04:08:24 CDT 2016

```Hi,

I'm trying to set a linear problem with a 2x2 block matrix (say A=[A11 A12;
A21 A22] in matlab notation). Obviously I could generate a single matrix,
but I would like to keep it as a 2x2 block matrix with the 4 blocks
generated and used independently (e.g., because at some point I need to
change some of those blocks to matrix-free). My plan was then to generate a
MatShell and set a MATOP_MULT shell operation to explicitly compute y=A*x by
splitting x=[x1;x2] and y=[y1;y2] and carrying out the single matrix-vector
multiplications, and reassembling the final result. My problem is that my
efforts in extracting x1 and x2 from the original x vector have failed. The
following Fortran code is what seemed to me the most reasonable thing to do
to set the matrix-vector multiplication, but it does not work

!======================================================
subroutine MyMult(A,x,y,IERR)

implicit none

#include <petsc/finclude/petscsys.h>
#include <petsc/finclude/petscvec.h>
#include <petsc/finclude/petscvec.h90>
#include <petsc/finclude/petscmat.h>

Mat     A
Vec     x,y, x1, x2, y1, y2
IS      :: ix1, ix2, iy1, iy2
PetscErrorCode ierr
PetscInt    M, N, II, NumberOfUnknowns, NumberOfMeasures
PetscInt, ALLOCATABLE         :: INDICES(:)
PetscMPIInt                   :: RANK

CALL VecGetSize(x,N,IERR)
NumberOfUnknowns = N / 2
CALL VecGetSize(y,N,IERR)
NumberOfMeasures = N - NumberOfUnknowns

N = NumberOfUnknowns + MAX(NumberOfUnknowns,NumberOfMeasures)
ALLOCATE(INDICES(N))
INDICES = (/ (ii, ii=0, N-1) /)

CALL
ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(1:NumberOfUnknowns),PETSC_COPY_VALUES,ix1,IERR)
CALL VecGetSubVector(x,ix1,x1,IERR)
CALL
ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(NumberOfUnknowns+1:2*NumberOfUnknowns),PETSC_COPY_VALUES,ix2,IERR)
CALL VecGetSubVector(x,ix2,x2,IERR)
CALL
ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(1:NumberOfUnknowns),PETSC_COPY_VALUES,iy1,IERR)
CALL VecGetSubVector(y,iy1,y1,IERR)
CALL
ISCreateGeneral(MPI_COMM_SELF,NumberOfUnknowns,INDICES(NumberOfUnknowns+1:NumberOfUnknowns+NumberOfMeasures),PETSC_COPY_VALUES,iy2,IERR)
CALL VecGetSubVector(y,iy2,y2,IERR)

CALL MatMult(L_discretized,x1,y1,IERR)
CALL MatMult(L_MR,x1,y2,IERR)

CALL VecRestoreSubVector(y,iy1,y1,IERR)
CALL VecRestoreSubVector(y,iy2,y2,IERR)

return
end subroutine MyMult
!======================================================

Obviously the sequence of calls to ISCreateGeneral and VecGetSubVector does
not do what I expect, as the errors I'm getting are in the following MatMult
multiplications (the global sizes of x1 and x2 SHOULD BE both 771, while y1
and y2 should be 771 and 2286):

1) executed with 1 MPI process:

[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Mat mat,Vec y: global dim 2286 771
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015
[0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named
alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr  6 10:54:03 2016
[0]PETSC ERROR: Configure options
--with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl
--with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc
--with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90
--with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec
--with-scalar-type=complex
[0]PETSC ERROR: #1 MatMult() line 2215 in
/home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Mat mat,Vec v3: local dim 2286 771
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015
[0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named
alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr  6 10:54:03 2016
[0]PETSC ERROR: Configure options
--with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl
--with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc
--with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90
--with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec
--with-scalar-type=complex
[0]PETSC ERROR: #2 MatMultAdd() line 2396 in
/home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c

2) executed with 4 MPI processes (I'm copying the output of process 0 only,
they're all the same):

....
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Mat mat,Vec x: global dim 771 3084
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015
[0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named
alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr  6 10:49:13 2016
[0]PETSC ERROR: Configure options
--with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl
--with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc
--with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90
--with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec
--with-scalar-type=complex
[0]PETSC ERROR: #1 MatMult() line 2214 in
/home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Mat mat,Vec v1: global dim 771 3084
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015
[0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named
alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr  6 10:49:13 2016
[0]PETSC ERROR: Configure options
--with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl
--with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc
--with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90
--with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec
--with-scalar-type=complex
[0]PETSC ERROR: #2 MatMultAdd() line 2393 in
/home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Mat mat,Vec x: global dim 771 3084
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015
[0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named
alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr  6 10:49:13 2016
[0]PETSC ERROR: Configure options
--with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl
--with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc
--with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90
--with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec
--with-scalar-type=complex
[0]PETSC ERROR: #3 MatMult() line 2214 in
/home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c
[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Nonconforming object sizes
[0]PETSC ERROR: Mat mat,Vec v1: global dim 771 3084
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for
trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015
[0]PETSC ERROR: ./ISS on a arch-linux2-c-debug named
alessandro-HP-ProDesk-490-G2-MT by alessandro Wed Apr  6 10:49:13 2016
[0]PETSC ERROR: Configure options
--with-blas-lapack-dir=/opt/intel/composer_xe_2015.2.164/mkl
--with-mpi-cc=/home/alessandro/.openmpi-1.8.3/bin/mpicc
--with-mpi-f90=/home/alessandro/.openmpi-1.8.3/bin/mpif90
--with-mpiexec=/home/alessandro/.openmpi-1.8.3/bin/mpiexec
--with-scalar-type=complex
[0]PETSC ERROR: #4 MatMultAdd() line 2393 in
/home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c
....

Any idea? I would like to keep the flexibility of having a sequence of
mat-vec multiplications as in

CALL MatMult(L_discretized,x1,y1,IERR)
CALL MatMult(L_MR,x1,y2,IERR)

but for doing that I need to be able to extract x1 and x2, and to reassemble
the partial results into y.

Alessandro
```