[petsc-users] About DMDA (and extracting its ordering)

Thu Feb 28 05:54:55 CST 2019

Hi Matthew,

I’ve just attached on that page https://bitbucket.org/petsc/petsc/issues/262 the same source code written in C that runs just fine from what I can see.

Which means there’s either a bug in ISLocalToGlobalMappingGetIndicesF90 or something missing in the documentation.

Best,

Thibaut

On 25 Feb 2019, at 12:03, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>> wrote:

On Mon, Feb 25, 2019 at 3:19 AM Appel, Thibaut <t.appel17 at imperial.ac.uk<mailto:t.appel17 at imperial.ac.uk>> wrote:
Hi Matthew,

Yes I need F90 and the syntax in the file / yours
PetscInt, pointer :: id_ltog(:)

Is the exact same as
PetscInt, dimension(:), pointer :: id_ltog

They’re both modern fortran ”correct”

Anyways I tried both and I still get the same error message.

Okay, that file I sent the link to is tested every night, so I suspect a compiler bug in your version.
Please send your entire example so I can run it here. Also, let us know what Fortran compiler
and version you are using again.

  Thanks,

     Matt

Thibaut

On 24 Feb 2019, at 23:38, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>> wrote:

On Sun, Feb 24, 2019 at 6:33 PM Appel, Thibaut <t.appel17 at imperial.ac.uk<mailto:t.appel17 at imperial.ac.uk>> wrote:
Is that for the id_ltog argument? I tried declaring it as IS, and you can’t declare an deferred-shape array as target. Same message. XXX(:) is syntactic sugar for dimension(:) :: XXX

I’m really confused because this
https://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex14f.F90.html
doesn’t use the xxxF90 variants of the routines we’re talking about?

You can use the F77 stuff if you want. I thought you needed F90.  For F90, I think you need to declare it

  PetscInt, pointer :: id_ltog(:)

like it is in the link I sent you.

  Thanks,

     Matt

Thibaut

On 24 Feb 2019, at 22:50, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>> wrote:

On Sun, Feb 24, 2019 at 4:28 PM Appel, Thibaut <t.appel17 at imperial.ac.uk<mailto:t.appel17 at imperial.ac.uk>> wrote:
Hi Matthew,

With the following data declaration and routines calls:

I am not a Fortran90 expert, but your declaration looks different than the ones I have seen work:

  https://bitbucket.org/petsc/petsc/src/master/src/dm/impls/plex/examples/tutorials/ex1f90.F90

  Thanks,

     Matt

  DM :: da
  ISLocalToGlobalMapping :: ltog
  PetscInt, DIMENSION(:), POINTER :: id_ltog

  CALL DMDACreate2D(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,(nx+1),(ny+1), &
                    PETSC_DECIDE,PETSC_DECIDE,neq,0,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
  CHKERRA(ierr)

  CALL DMSetUp(da,ierr)
  CHKERRA(ierr)

  CALL DMDAGetCorners(da,icx,icy,PETSC_NULL_INTEGER,lx,ly,PETSC_NULL_INTEGER,ierr)
  CHKERRA(ierr)

  CALL DMGetLocalToGlobalMapping(da,ltog,ierr)
  CHKERRA(ierr)

  CALL ISLocalToGlobalMappingGetIndicesF90(ltog,id_ltog,ierr)
  CHKERRA(ierr)

  CALL ISLocalToGlobalMappingRestoreIndicesF90(ltog,id_ltog,ierr)
  CHKERRA(ierr)

  CALL DMDestroy(da,ierr)
  CHKERRA(ierr)

I get, with the most recent Intel Fortran compiler and a fresh “git pull” PETsc:

forrtl: severe (408): fort: (7): Attempt to use pointer ID_LTOG when it is not associated with a target

Thibaut

On 22 Feb 2019, at 15:13, Matthew Knepley <knepley at gmail.com<mailto:knepley at gmail.com>> wrote:

On Fri, Feb 22, 2019 at 9:10 AM Thibaut Appel via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> wrote:
I reckon that what corresponds best to my situation is to create my own MPIAIJ matrix, recover the DMDA local to global mapping, preallocate myself and fill my values with MatSetValues.

However, what is the correct way to call ISLocalToGlobalMappingGetIndices in Fortran? It seems there's something missing in the documentation as my code just won't work.

Some examples include a PetscOffset argument which is not present in the documentation and there is no manual page for ISLocalToGlobalMappingGetIndicesF90: the only information I could find is on /src/vec/is/utils/f90-custom/zisltogf90.c. Are you supposed to allocate the index array yourself? I tried declaring it as a pointer too. I keep getting an error

1) There was a bug in that file, but I don't think it should have affected you. I will put in a PR

2)  You are supposed to pass in a pointer, unallocated. We manage all the allocation, which is why you call Restore. What happens when you do that?

  Thanks,

     Matt
[0]PETSC ERROR: #1 User provided function() line 0 in User file
application called MPI_Abort(MPI_COMM_SELF, -42924) - process 0

Minimal example:

PROGRAM test_dmda

#include <petsc/finclude/petscdm.h>
#include <petsc/finclude/petscdmda.h>

  USE PetscDMDA

  IMPLICIT NONE

  PetscErrorCode :: ierr
  PetscInt       :: irank, icx, icy, lx, ly, n_ltog
  DM             :: da
  ISLocalToGlobalMapping :: ltog

  PetscInt, DIMENSION(:), ALLOCATABLE :: id_ltog
  !PetscInt, DIMENSION(:), POINTER :: id_ltog

  PetscInt, PARAMETER :: nx=24, ny=9, neq=4
  CHARACTER(LEN=100) :: err_msg

  CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)

  CALL MPI_COMM_RANK(PETSC_COMM_WORLD,irank,ierr)

  CALL DMDACreate2D(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,(nx+1),(ny+1), &
                    PETSC_DECIDE,PETSC_DECIDE,neq,0,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
  CHKERRA(ierr)

  CALL DMSetUp(da,ierr)
  CHKERRA(ierr)

  CALL DMDAGetCorners(da,icx,icy,PETSC_NULL_INTEGER,lx,ly,PETSC_NULL_INTEGER,ierr)
  CHKERRA(ierr)

  WRITE(*,'(*(1X,A,I3))') 'I am process #', irank, ' - my lower left corner is icx=', icx, &
                          'and icy=', icy, ' - the width is lx=', lx, ' and ly=', ly

  CALL DMGetLocalToGlobalMapping(da,ltog,ierr)
  CHKERRA(ierr)

  CALL ISLocalToGlobalMappingGetSize(ltog,n_ltog,ierr)
  CHKERRA(ierr)

  ALLOCATE(id_ltog(1:n_ltog),STAT=ierr,ERRMSG=err_msg)

  CALL ISLocalToGlobalMappingGetIndices(ltog,id_ltog,ierr)
  CHKERRA(ierr)

  CALL ISLocalToGlobalMappingRestoreIndices(ltog,id_ltog,ierr)
  CHKERRA(ierr)

  DEALLOCATE(id_ltog,STAT=ierr,ERRMSG=err_msg)

  CALL DMDestroy(da,ierr)
  CHKERRA(ierr)

  CALL PetscFinalize(ierr)

END PROGRAM test_dmda

On 21/02/2019 17:49, Matthew Knepley wrote:
On Thu, Feb 21, 2019 at 11:16 AM Thibaut Appel via petsc-users <petsc-users at mcs.anl.gov<mailto:petsc-users at mcs.anl.gov>> wrote:
Dear PETSc developers/users,

I’m solving linear PDEs on a regular grid with high-order finite differences, assembling an MPIAIJ matrix to solve linear systems or eigenvalue problems. I’ve been using vertex major, natural ordering for the parallelism with PetscSplitOwnership (yielding rectangular slices of the physical domain) and wanted to move to DMDA to have a more square-ish domain decomposition and minimize communication between processes.

However, my application is memory critical, and I have finely-tuned matrix preallocation routines for allocating memory “optimally”. It seems the memory of a DMDA matrix is allocated along the value of the stencil width of DMDACreate and the manual says about it

“These DMDA stencils have nothing directly to do with any finite difference stencils one might chose to use for a discretization”

And despite reading the manual pages there must be something I do not understand in the DM topology, what is that "stencil width" for then? I will not use ghost values for my FD-method, right?

What this is saying is, "You might be using some stencil that is not STAR or BOX, but we are preallocating according to one of those".
If you really care about how much memory is preallocated, which it seems you do, then you might be able to use

  https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMDA/DMDASetBlockFills.html

to tell use exactly how to preallocate.

I was then wondering if I could just create a MPIAIJ matrix, and with a PETSc routine get the global indices of the domain for each process: in other words, an equivalent of PetscSplitOwnership that gives me the DMDA unknown ordering. So I can feed and loop on that in my preallocation and assembly routines.

You can make an MPIAIJ matrix yourself of course. It should have the same division of rows as the DMDA division of dofs. Also, MatSetValuesStencil() will not work for a custom matrix.

  Thanks,

     Matt

Thanks very much,

Thibaut

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

--
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/<http://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/<http://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/<http://www.cse.buffalo.edu/~knepley/>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20190228/95257278/attachment-0001.html>