[petsc-users] MatMPIAIJGetSeqAIJ implemented for fortran?

[Smith, Barry F.]

    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 <mbuerkle at web.de> 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 <mbuerkle at web.de> 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 <mbuerkle at web.de> 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 <mbuerkle at web.de> 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 <mbuerkle at web.de> 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/
>