<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
</head>
<body text="#000000" bgcolor="#FFFFFF">
<p>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.</p>
<p>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.</p>
<p>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</p>
<p>[0]PETSC ERROR: #1 User provided function() line 0 in User file<br>
application called MPI_Abort(MPI_COMM_SELF, -42924) - process 0</p>
<p>Minimal example:</p>
<p><tt>PROGRAM test_dmda<br>
<br>
#include <petsc/finclude/petscdm.h><br>
#include <petsc/finclude/petscdmda.h><br>
<br>
USE PetscDMDA<br>
<br>
IMPLICIT NONE<br>
<br>
PetscErrorCode :: ierr<br>
PetscInt :: irank, icx, icy, lx, ly, n_ltog<br>
DM :: da<br>
ISLocalToGlobalMapping :: ltog<br>
<br>
PetscInt, DIMENSION(:), ALLOCATABLE :: id_ltog<br>
!PetscInt, DIMENSION(:), POINTER :: id_ltog<br>
<br>
PetscInt, PARAMETER :: nx=24, ny=9, neq=4<br>
CHARACTER(LEN=100) :: err_msg<br>
<br>
CALL PetscInitialize(PETSC_NULL_CHARACTER,ierr)<br>
<br>
CALL MPI_COMM_RANK(PETSC_COMM_WORLD,irank,ierr)<br>
<br>
CALL
DMDACreate2D(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,(nx+1),(ny+1),
&<br>
PETSC_DECIDE,PETSC_DECIDE,neq,0,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)<br>
CHKERRA(ierr)<br>
<br>
CALL DMSetUp(da,ierr)<br>
CHKERRA(ierr)<br>
<br>
CALL
DMDAGetCorners(da,icx,icy,PETSC_NULL_INTEGER,lx,ly,PETSC_NULL_INTEGER,ierr)<br>
CHKERRA(ierr)<br>
<br>
WRITE(*,'(*(1X,A,I3))') 'I am process #', irank, ' - my lower
left corner is icx=', icx, &<br>
'and icy=', icy, ' - the width is
lx=', lx, ' and ly=', ly<br>
<br>
CALL DMGetLocalToGlobalMapping(da,ltog,ierr)<br>
CHKERRA(ierr)<br>
<br>
CALL ISLocalToGlobalMappingGetSize(ltog,n_ltog,ierr)<br>
CHKERRA(ierr)<br>
<br>
ALLOCATE(id_ltog(1:n_ltog),STAT=ierr,ERRMSG=err_msg)<br>
<br>
CALL ISLocalToGlobalMappingGetIndices(ltog,id_ltog,ierr)<br>
CHKERRA(ierr)<br>
<br>
CALL ISLocalToGlobalMappingRestoreIndices(ltog,id_ltog,ierr)<br>
CHKERRA(ierr)<br>
<br>
DEALLOCATE(id_ltog,STAT=ierr,ERRMSG=err_msg)<br>
<br>
CALL DMDestroy(da,ierr)<br>
CHKERRA(ierr)<br>
<br>
CALL PetscFinalize(ierr)<br>
<br>
END PROGRAM test_dmda</tt><br>
</p>
<p><br>
</p>
On 21/02/2019 17:49, Matthew Knepley wrote:<br>
<blockquote type="cite"
cite="mid:CAMYG4GnUKkjLm4s-WrexvmtOW0h7GHzH4=qnruWdad-6pnFrdg@mail.gmail.com">
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<div dir="ltr">
<div dir="ltr">
<div dir="ltr">On Thu, Feb 21, 2019 at 11:16 AM Thibaut Appel
via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov"
moz-do-not-send="true">petsc-users@mcs.anl.gov</a>>
wrote:<br>
</div>
<div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex">
<div bgcolor="#FFFFFF"> <font size="-1" face="Arial,
sans-serif"><span style="font-size:10pt">Dear PETSc
developers/users,</span></font><br>
<font size="-1" face="Arial, sans-serif"> </font><br>
<font size="-1" face="Arial, sans-serif"><span
style="font-size:10pt">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.</span></font><br>
<font size="-1" face="Arial, sans-serif"><span
style="font-size:10pt"></span></font><br>
<font size="-1" face="Arial, sans-serif"><span
style="font-size:10pt">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</span></font><br>
<font size="-1" face="Arial, sans-serif"> </font><br>
<font size="-1" face="Arial, sans-serif"><span
style="font-size:10pt">“</span>These DMDA stencils
have nothing directly to do with any finite difference
stencils one might chose to use for a discretization”</font><br>
<font size="-1" face="Arial, sans-serif"> </font><br>
<font size="-1" face="Arial, sans-serif">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?</font></div>
</blockquote>
<div><br>
</div>
<div>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".</div>
<div>If you really care about how much memory is
preallocated, which it seems you do, then you might be
able to use</div>
<div><br>
</div>
<div> <a
href="https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMDA/DMDASetBlockFills.html"
moz-do-not-send="true">https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/DMDA/DMDASetBlockFills.html</a></div>
<div><br>
</div>
<div>to tell use exactly how to preallocate.</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex">
<div bgcolor="#FFFFFF"> <font size="-1" face="Arial,
sans-serif">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.</font></div>
</blockquote>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div> Thanks,</div>
<div><br>
</div>
<div> Matt</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left:1px solid
rgb(204,204,204);padding-left:1ex">
<div bgcolor="#FFFFFF"> <font size="-1" face="Arial,
sans-serif">T</font><font size="-1" face="Arial,
sans-serif">hanks very much,<span
style="font-size:10pt"></span></font><br>
<font size="-1" face="Arial, sans-serif"><span
style="font-size:10pt"> </span></font><br>
<font size="-1" face="Arial, sans-serif"> <span
style="font-size:10pt">Thibaut</span></font> </div>
</blockquote>
</div>
<br clear="all">
<div><br>
</div>
-- <br>
<div dir="ltr" class="gmail_signature">
<div dir="ltr">
<div>
<div dir="ltr">
<div>
<div dir="ltr">
<div>What most experimenters take for granted
before they begin their experiments is
infinitely more interesting than any results to
which their experiments lead.<br>
-- Norbert Wiener</div>
<div><br>
</div>
<div><a
href="http://www.cse.buffalo.edu/~knepley/"
target="_blank" moz-do-not-send="true">https://www.cse.buffalo.edu/~knepley/</a><br>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</body>
</html>