On 22 Feb 2019, at 15:13, Matthew Knepley
<knep...@gmail.com <mailto:knep...@gmail.com>> wrote:
On Fri, Feb 22, 2019 at 9:10 AM Thibaut Appel via
petsc-users <petsc-users@mcs.anl.gov
<mailto:petsc-users@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@mcs.anl.gov
<mailto:petsc-users@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/
