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
#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
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
& 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
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/
> 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/
