On 6 April 2016 at 11:08, FRANCAVILLA MATTEO ALESSANDRO <d019...@polito.it> wrote:
> 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. Take a look at MatNest - it does exactly what you want. http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateNest.html http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MATNEST.html#MATNEST Thanks, Dave > 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 MatMultAdd(K_discretized,x2,y1,y1,IERR) > CALL MatMult(L_MR,x1,y2,IERR) > CALL MatMultAdd(K_MR,x2,y2,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 MatMultAdd(K_discretized,x2,y1,y1,IERR) > CALL MatMult(L_MR,x1,y2,IERR) > CALL MatMultAdd(K_MR,x2,y2,y2,IERR) > > but for doing that I need to be able to extract x1 and x2, and to > reassemble the partial results into y. > > Thanks in advance, > Alessandro >