Hi Matthew,
With the following data declaration and routines calls:
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
<[email protected]<mailto:[email protected]>> wrote:
On Fri, Feb 22, 2019 at 9:10 AM Thibaut Appel via petsc-users
<[email protected]<mailto:[email protected]>> 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
<[email protected]> 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/
