[petsc-users] How to extract array portions
Dave May
dave.mayhem23 at gmail.com
Wed Apr 6 09:27:55 CDT 2016
On 6 April 2016 at 16:20, FRANCAVILLA MATTEO ALESSANDRO <d019598 at polito.it>
wrote:
> Fantastic, it definitely solved my problem! Thank you!
>
> Should I always create the vectors with MatCreateVecs when I use them for
> matrix-vector multiplications? I used to create vectors with VecCreate, and
> it has always worked fine with "standard" matrices (matrices not created
> with MATNEST structure).
>
Yes
>
> In the example at
> http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex2.c.html
> I read the following:
>
> "PETSc automatically generates appropriately partitioned matrices and
> vectors when MatCreate() and VecCreate() are used with the same
> communicator."
>
This comment also assumes that you use PETSC_DECIDE for the local sizes. So
if you create a matrix with global size M x M and a vector of global size M
and you (i) use PETSC_DECIDE for the local size for the row/col with
MatCreate() and PETSC_DECIDE for the local row for VecCreate(); and (ii)
you use the same communicator - then the row partitioning adopted by the
matrix will be consistent with the row partitioning for the vector.
> It is not clear to me when I can simply use VecCreate (as I've always done
> so far), and when to use MatCreateVecs...
>
You should always use MatCreateVecs(). The implementation of Mat will give
you back a consistent vector.
Thanks,
Dave
>
>
>
> On Wed, 6 Apr 2016 14:38:12 +0200
>
> Dave May <dave.mayhem23 at gmail.com> wrote:
>
>> On 6 April 2016 at 14:12, FRANCAVILLA MATTEO ALESSANDRO <
>> d019598 at polito.it>
>> wrote:
>>
>> Thanks Dave, that's exactly what I was looking for, and it (partially)
>>> solved my problem.
>>>
>>> The MatMult operation now works fine, unfortunately not for every number
>>> of MPI processes (while the MatMult involving each submatrix works fine
>>> for
>>> any number of processes):
>>>
>>> ....
>>> matArray(1) = L_discretized
>>> matArray(2) = K_discretized
>>> matArray(3) = L_MR
>>> matArray(4) = K_MR
>>> CALL
>>>
>>> MatCreateNest(PETSC_COMM_WORLD,2,PETSC_NULL_OBJECT,2,PETSC_NULL_OBJECT,matArray,TheMatrix,IERR)
>>> CALL MatSetFromOptions(TheMatrix,IERR)
>>> !
>>> CALL VecCreate(PETSC_COMM_WORLD,x,IERR)
>>> CALL VecSetSizes(x,PETSC_DECIDE,2*NumberOfUnknowns,IERR)
>>> CALL VecSetFromOptions(x,IERR)
>>> CALL VecSet(x,0.+0.*PETSC_i,IERR)
>>> CALL VecSetValues(x,1,0,1.+0.*PETSC_i,INSERT_VALUES,IERR)
>>> CALL VecAssemblyBegin(x,IERR)
>>> CALL VecAssemblyEnd(x,IERR)
>>> CALL VecCreate(PETSC_COMM_WORLD,y,IERR)
>>> CALL VecSetSizes(y,PETSC_DECIDE,NumberOfUnknowns+NumberOfMeasures,IERR)
>>> CALL VecSetFromOptions(y,IERR)
>>> CALL MatMult(TheMatrix,x,y,IERR)
>>> ....
>>>
>>>
>>> I guess there is something wrong with the distribution of vectors between
>>> the MPI processes. With the specific sizes of my problem everything works
>>> fine with 1-2-3-6 processes, otherwise I get "Nonconforming object sizes"
>>> error, e.g. with 8 processes:
>>>
>>> [2]PETSC ERROR: --------------------- Error Message
>>> --------------------------------------------------------------
>>> [2]PETSC ERROR: Nonconforming object sizes
>>> [2]PETSC ERROR: Mat mat,Vec y: local dim 383 382
>>> [2]PETSC ERROR: Petsc Release Version 3.6.3, Dec, 03, 2015
>>> [2]PETSC ERROR: Nonconforming object sizes
>>> ./ISS on a arch-linux2-c-debug named alessandro-HP-ProDesk-490-G2-MT by
>>> alessandro Wed Apr 6 13:31:45 2016
>>> [2]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
>>> [2]PETSC ERROR: #1 MatMult() line 2216 in
>>> /home/alessandro/CODE/PETSc/petsc-3.6.3/src/mat/interface/matrix.c
>>> ...
>>> [6]PETSC ERROR: Mat mat,Vec y: local dim 381 382
>>> ...
>>> [7]PETSC ERROR: Mat mat,Vec y: local dim 381 382
>>> ...
>>> [1]PETSC ERROR: Mat mat,Vec y: local dim 383 382
>>> ...
>>>
>>>
>>> I noticed that the code works (runs and provides correct results) when
>>> the
>>> number of processes divides the length of y2 (number of rows of A21 and
>>> A22); however when I apply MatMult separately on A21 and A22 everything
>>> works fine with any number of processes
>>>
>>> Any idea where I'm messing up here?
>>>
>>>
>> Sounds like the large vector you are using to multiple with the matnest
>> object wasn't created via MatCreateVecs(). If you use this method you will
>> obtain a vector with the correct local and global sizes.
>>
>>
>>
>>
>>
>>>
>>> Thanks,
>>> Alessandro
>>>
>>>
>>>
>>>
>>> On Wed, 6 Apr 2016 11:16:43 +0200
>>> Dave May <dave.mayhem23 at gmail.com> wrote:
>>>
>>> On 6 April 2016 at 11:08, FRANCAVILLA MATTEO ALESSANDRO <
>>>> d019598 at 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
>>>>>
>>>>>
>>>>> ___________________________________
>>> Matteo Alessandro Francavilla, PhD
>>> Antenna and EMC Lab (LACE)
>>> Istituto Superiore Mario Boella (ISMB)
>>> Torino, Italy
>>> phone +39 011 2276 711
>>>
>>>
> ___________________________________
> Matteo Alessandro Francavilla, PhD
> Antenna and EMC Lab (LACE)
> Istituto Superiore Mario Boella (ISMB)
> Torino, Italy
> phone +39 011 2276 711
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160406/2a0e420a/attachment-0001.html>
More information about the petsc-users
mailing list