The communicator used to construct the IS must be the same communicator as
used to construct the matrix, in your case PETSC_COMM_WORLD.
Barry
I will add an error check in the maint branch to catch these bugs in the future.
It just so happens that you can use the PETSC_COMM_SELF communicator for AIJ
matrices.
> On Jul 11, 2018, at 9:54 PM, Marius Buerkle <[email protected]> wrote:
>
> here some example code which works for MPIAIJ but fails for MPIDENSE. Maybe
> you can tell me what I am doing wrong
>
>
> program submatrix
> #include <petsc/finclude/petscmat.h>
> use petscmat
> implicit none
>
> Mat :: A,B
> PetscInt ::n,m,nl1,nl2,nlc1,nlc2,nproc,iproc,i,ncols,nrow,ncol
> PetscReal :: info(MAT_INFO_SIZE)
> PetscErrorCode :: ierr
> IS :: isrow,iscol
> integer,allocatable :: irow(:),icol(:),cols(:)
> character(256) :: stype,sdummy
>
> call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
>
> call MPI_COMM_SIZE(PETSC_COMM_WORLD, nproc,ierr)
>
> call MPI_Comm_rank(PETSC_COMM_WORLD, iproc, ierr )
>
> if (nproc.ne.3) then
> if (iproc.eq.0) write(0,*) "nproc.ne.3 ",nproc
> call PetscFinalize(ierr)
> stop
> end if
>
> n=3
> m=n
>
> call
> MatCreateDense(PETSC_COMM_WORLD,1,PETSC_DECIDE,n,m,PETSC_NULL_SCALAR,A,ierr)
> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
>
> allocate(irow(1),icol(1))
>
> !~ call MatConvert(A,MATAIJ, MAT_INPLACE_MATRIX,A,ierr)
>
> nrow=0
> ncol=1
> icol=iproc
> if (iproc.eq.0) then
> irow=iproc
> nrow=1
> end if
> call ISCreateGeneral(PETSC_COMM_SELF,nrow,irow,&
> & PETSC_COPY_VALUES,isrow,ierr)
>
> call ISCreateGeneral(PETSC_COMM_SELF,ncol,icol,&
> & PETSC_COPY_VALUES,iscol,ierr)
>
> call MatCreateSubMatrix(A,isrow,iscol,MAT_INITIAL_MATRIX,B,ierr)
>
> call MatGetType(B,stype,ierr)
> call MatGetInfo(B, MAT_LOCAL, info,ierr)
> call MatGetOwnershipRange(B,nl1,nl2,ierr)
> call MatGetOwnershipRangeColumn(B,nlc1,nlc2,ierr)
> call MatGetSize(B,m,n,ierr)
> write(sdummy,fmt='(A,25i16)')
> trim(stype),iproc,int(info(mat_info_nz_allocated)),int(info(mat_info_nz_used)),nl1,nl2,nlc1,nlc2,m,n
> call
> PetscSynchronizedPrintf(PETSC_COMM_WORLD,trim(sdummy)//NEW_LINE('A'),ierr)
> call PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT,ierr)
>
>
> call PetscFinalize(ierr)
>
> end program submatrix
>
>
>
>
>
> i see. but i doesn't work for me. maybe I made some mistake. If not I will
> try to make some example code where it fails.
>
>
>
> > On Jul 10, 2018, at 11:50 PM, Marius Buerkle <[email protected]> wrote:
> >
> > Sorry, I have yet another question. While it works well for MPIAIJ, I end
> > up in a deadlock if I use MatCreateSubmatrix on a MPIDENSE matrix. Should
> > MatCreateSubmatrix work on dense matricies ?
>
> yes
>
> > On Mon, Jul 9, 2018 at 10:02 PM Marius Buerkle <[email protected]> wrote:
> > MatGetSubmatrix is in the current PETSc released called MatCreateSubmatrix.
> > Is this correct ?
> >
> >
> > Yep. I am old and cannot follow API changes anymore.
> >
> > Matt
> >
> > On Mon, Jul 9, 2018 at 8:38 PM Marius Buerkle <[email protected]> wrote:
> > I see. What I want to do is to calculate the matrix product C=A*B' between
> > two sparse matrices A=(A11 0 , A21 0) and B=(B11 0 , B21 0) where C will be
> > dense in the end but I just want to calculate some selected entries C_ij of
> > C. At the moment I extract submatricies for the corresponding rows and
> > columns,
> >
> > I think that is the right way, but you should only need MatGetSubmatrix for
> > that.
> >
> > Thanks,
> >
> > Matt
> >
> > so I was wondering if there is a simpler or performancer-wise faster way. I
> > assume there is not such thing as a restricted MatMatMul which just
> > calculated the lets say predefined nonzero entries of C.
> >
> >
> > These are "internal" routines that we rarely expect end users to use since
> > they are specific for particular matrix implementations. As such they are
> > also kind of strange to use from Fortran since none of the AIJ data
> > structures can be made visible to Fortran.
> >
> > Could you explain why you want them from Fortran and maybe we'll have
> > alternative suggestions on how you can achieve the same effect.
> >
> > Barry
> >
> >
> > > On Jul 5, 2018, at 3:05 AM, Marius Buerkle <[email protected]> wrote:
> > >
> > > or MatMPIAIJGetLocalMatCondensed for that matter.
> > >
> > >
> > >
> > > Hi !
> > >
> > > Is MatMPIAIJGetSeqAIJ implemented for fortran?
> > >
> > > best,
> > > Marius
> >
> >
> >
> > --
> > What most experimenters take for granted before they begin their
> > experiments is infinitely more interesting than any results to which their
> > experiments lead.
> > -- Norbert Wiener
> >
> > https://www.cse.buffalo.edu/~knepley/
> >
> >
> > --
> > What most experimenters take for granted before they begin their
> > experiments is infinitely more interesting than any results to which their
> > experiments lead.
> > -- Norbert Wiener
> >
> > https://www.cse.buffalo.edu/~knepley/
>
