[petsc-users] Fwd: Understanding index sets for PCGASM
LEONARDO MUTTI
leonardo.mutti01 at universitadipavia.it
Tue May 9 11:31:29 CDT 2023
---------- Forwarded message ---------
Da: LEONARDO MUTTI <leonardo.mutti01 at universitadipavia.it>
Date: mar 9 mag 2023 alle ore 18:29
Subject: Re: [petsc-users] Understanding index sets for PCGASM
To: Matthew Knepley <knepley at gmail.com>
Thank you for your answer, but I am still confused, sorry.
Consider
https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90 on
one processor.
Let M=12 for the sake of simplicity, i.e. we deal with a 12x12 2D grid,
hence, a 144x144 matrix.
Let NSubx = 3, so that on the grid we do 3 vertical and 3 horizontal
subdivisions.
We should obtain 9 subdomains that are grids of 4x4 nodes each, thus
corresponding to 9 submatrices of size 16x16.
In my run I obtain NSub = 9 (great) and subdomain_IS(i), i=1,...,9, reads:
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 0*
*1 1*
*2 2*
*3 3*
*4 12*
*5 13*
*6 14*
*7 15*
*8 24*
*9 25*
*10 26*
*11 27*
*12 36*
*13 37*
*14 38*
*15 39*
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 4*
*1 5*
*2 6*
*3 7*
*4 16*
*5 17*
*6 18*
*7 19*
*8 28*
*9 29*
*10 30*
*11 31*
*12 40*
*13 41*
*14 42*
*15 43*
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 8*
*1 9*
*2 10*
*3 11*
*4 20*
*5 21*
*6 22*
*7 23*
*8 32*
*9 33*
*10 34*
*11 35*
*12 44*
*13 45*
*14 46*
*15 47*
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 48*
*1 49*
*2 50*
*3 51*
*4 60*
*5 61*
*6 62*
*7 63*
*8 72*
*9 73*
*10 74*
*11 75*
*12 84*
*13 85*
*14 86*
*15 87*
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 52*
*1 53*
*2 54*
*3 55*
*4 64*
*5 65*
*6 66*
*7 67*
*8 76*
*9 77*
*10 78*
*11 79*
*12 88*
*13 89*
*14 90*
*15 91*
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 56*
*1 57*
*2 58*
*3 59*
*4 68*
*5 69*
*6 70*
*7 71*
*8 80*
*9 81*
*10 82*
*11 83*
*12 92*
*13 93*
*14 94*
*15 95*
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 96*
*1 97*
*2 98*
*3 99*
*4 108*
*5 109*
*6 110*
*7 111*
*8 120*
*9 121*
*10 122*
*11 123*
*12 132*
*13 133*
*14 134*
*15 135*
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 100*
*1 101*
*2 102*
*3 103*
*4 112*
*5 113*
*6 114*
*7 115*
*8 124*
*9 125*
*10 126*
*11 127*
*12 136*
*13 137*
*14 138*
*15 139*
*IS Object: 1 MPI process*
* type: general*
*Number of indices in set 16*
*0 104*
*1 105*
*2 106*
*3 107*
*4 116*
*5 117*
*6 118*
*7 119*
*8 128*
*9 129*
*10 130*
*11 131*
*12 140*
*13 141*
*14 142*
*15 143*
As you said, no number here reaches 144.
But the number stored in subdomain_IS are 9x16= #subdomains x 16, whereas I
would expect, also given your latest reply, 9x16x16x2=#subdomains x
submatrix height x submatrix width x length of a (row,column) pair.
It would really help me if you could briefly explain how the output above
encodes the subdivision into subdomains.
Many thanks again,
Leonardo
Il giorno mar 9 mag 2023 alle ore 16:24 Matthew Knepley <knepley at gmail.com>
ha scritto:
> On Tue, May 9, 2023 at 10:05 AM LEONARDO MUTTI <
> leonardo.mutti01 at universitadipavia.it> wrote:
>
>> Great thanks! I can now successfully run
>> https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/ex71f.F90.
>>
>> Going forward with my experiments, let me post a new code snippet (very
>> similar to ex71f.F90) that I cannot get to work, probably I must be
>> setting up the IS objects incorrectly.
>>
>> I have an 8x8 matrix A=diag(1,1,2,2,...,2) and a vector b=(0.5,...,0.5).
>> We have only one processor, and I want to solve Ax=b using GASM. In
>> particular, KSP is set to preonly, GASM is the preconditioner and it uses
>> on each submatrix an lu direct solver (sub_ksp = preonly, sub_pc = lu).
>>
>> For the GASM algorithm, I divide A into diag(1,1) and diag(2,2,...,2).
>> For simplicity I set 0 overlap. Now I want to use GASM to solve Ax=b. The
>> code follows.
>>
>> #include <petsc/finclude/petscmat.h>
>> #include <petsc/finclude/petscksp.h>
>> #include <petsc/finclude/petscpc.h>
>> USE petscmat
>> USE petscksp
>> USE petscpc
>> USE MPI
>>
>> Mat :: A
>> Vec :: b, x
>> PetscInt :: M, I, J, ISLen, NSub
>> PetscMPIInt :: size
>> PetscErrorCode :: ierr
>> PetscScalar :: v
>> KSP :: ksp
>> PC :: pc
>> IS :: subdomains_IS(2), inflated_IS(2)
>> PetscInt,DIMENSION(4) :: indices_first_domain
>> PetscInt,DIMENSION(36) :: indices_second_domain
>>
>> call PetscInitialize(PETSC_NULL_CHARACTER, ierr)
>> call MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)
>>
>>
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> ! INTRO: create matrix and right hand side
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>>
>> WRITE(*,*) "Assembling A,b"
>>
>> M = 8
>> call MatCreateAIJ(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,
>> & M, M, PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,
>> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,A, ierr)
>> DO I=1,M
>> DO J=1,M
>> IF ((I .EQ. J) .AND. (I .LE. 2 )) THEN
>> v = 1
>> ELSE IF ((I .EQ. J) .AND. (I .GT. 2 )) THEN
>> v = 2
>> ELSE
>> v = 0
>> ENDIF
>> call MatSetValue(A, I-1, J-1, v, INSERT_VALUES, ierr)
>> END DO
>> END DO
>>
>> call MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)
>> call MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)
>>
>> call VecCreate(PETSC_COMM_WORLD,b,ierr)
>> call VecSetSizes(b, PETSC_DECIDE, M,ierr)
>> call VecSetFromOptions(b,ierr)
>>
>> do I=1,M
>> v = 0.5
>> call VecSetValue(b,I-1,v, INSERT_VALUES,ierr)
>> end do
>>
>> call VecAssemblyBegin(b,ierr)
>> call VecAssemblyEnd(b,ierr)
>>
>>
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> ! FIRST KSP/PC SETUP
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>>
>> WRITE(*,*) "KSP/PC first setup"
>>
>> call KSPCreate(PETSC_COMM_WORLD, ksp, ierr)
>> call KSPSetOperators(ksp, A, A, ierr)
>> call KSPSetType(ksp, 'preonly', ierr)
>> call KSPGetPC(ksp, pc, ierr)
>> call KSPSetUp(ksp, ierr)
>> call PCSetType(pc, PCGASM, ierr)
>>
>>
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> ! GASM, SETTING SUBDOMAINS
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>>
>> WRITE(*,*) "Setting GASM subdomains"
>>
>> ! Let's create the subdomain IS and inflated_IS
>> ! They are equal if no overlap is present
>> ! They are 1: 0,1,8,9
>> ! 2: 10,...,15,18,...,23,...,58,...,63
>>
>> indices_first_domain = [0,1,8,9] ! corresponds to diag(1,1)
>> do I=0,5
>> do J=0,5
>> indices_second_domain(I*6+1+J) = 18 + J + 8*I ! corresponds
>> to diag(2,2,...,2)
>> !WRITE(*,*) I*6+1+J, 18 + J + 8*I
>> end do
>> end do
>>
>> ! Convert into IS
>> ISLen = 4
>> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain,
>> & PETSC_COPY_VALUES, subdomains_IS(1), ierr)
>> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_first_domain,
>> & PETSC_COPY_VALUES, inflated_IS(1), ierr)
>> ISLen = 36
>> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain,
>> & PETSC_COPY_VALUES, subdomains_IS(2), ierr)
>> call ISCreateGeneral(PETSC_COMM_WORLD,ISLen,indices_second_domain,
>> & PETSC_COPY_VALUES, inflated_IS(2), ierr)
>>
>> NSub = 2
>> call PCGASMSetSubdomains(pc,NSub,
>> & subdomains_IS,inflated_IS,ierr)
>> call PCGASMDestroySubdomains(NSub,
>> & subdomains_IS,inflated_IS,ierr)
>>
>>
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> ! GASM: SET SUBSOLVERS
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>>
>> WRITE(*,*) "Setting subsolvers for GASM"
>>
>> call PCSetUp(pc, ierr) ! should I add this?
>>
>> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,
>> & "-sub_pc_type", "lu", ierr)
>> call PetscOptionsSetValue(PETSC_NULL_OPTIONS,
>> & "-sub_ksp_type", "preonly", ierr)
>>
>> call KSPSetFromOptions(ksp, ierr)
>> call PCSetFromOptions(pc, ierr)
>>
>>
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>> ! DUMMY SOLUTION: DID IT WORK?
>> !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
>>
>> WRITE(*,*) "Solve"
>>
>> call VecDuplicate(b,x,ierr)
>> call KSPSolve(ksp,b,x,ierr)
>>
>> call MatDestroy(A, ierr)
>> call KSPDestroy(ksp, ierr)
>> call PetscFinalize(ierr)
>>
>> This code is failing in multiple points. At call PCSetUp(pc, ierr) it
>> produces:
>>
>> *[0]PETSC ERROR: Argument out of range*
>> *[0]PETSC ERROR: Scatter indices in ix are out of range*
>> *...*
>> *[0]PETSC ERROR: #1 VecScatterCreate() at
>> ***\src\vec\is\sf\INTERF~1\vscat.c:736*
>> *[0]PETSC ERROR: #2 PCSetUp_GASM() at
>> ***\src\ksp\pc\impls\gasm\gasm.c:433*
>> *[0]PETSC ERROR: #3 PCSetUp() at ***\src\ksp\pc\INTERF~1\precon.c:994*
>>
>> And at call KSPSolve(ksp,b,x,ierr) it produces:
>>
>> *forrtl: severe (157): Program Exception - access violation*
>>
>>
>> The index sets are setup coherently with the outputs of e.g.
>> https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/ksp/tests/output/ex71f_1.out:
>> in particular each element of the matrix A corresponds to a number from 0
>> to 63.
>>
>
> This is not correct, I believe. The indices are row/col indices, not
> indices into dense blocks, so for
> your example, they are all in [0, 8].
>
> Thanks,
>
> Matt
>
>
>> Note that each submatrix does not represent some physical subdomain, the
>> subdivision is just at the algebraic level.
>> I thus have the following questions:
>>
>> - is this the correct way of creating the IS objects, given my
>> objective at the beginning of the email? Is the ordering correct?
>> - what am I doing wrong that is generating the above errors?
>>
>> Thanks for the patience and the time.
>> Best,
>> Leonardo
>>
>> Il giorno ven 5 mag 2023 alle ore 18:43 Barry Smith <bsmith at petsc.dev>
>> ha scritto:
>>
>>>
>>> Added in *barry/2023-05-04/add-pcgasm-set-subdomains *see also
>>> https://gitlab.com/petsc/petsc/-/merge_requests/6419
>>>
>>> Barry
>>>
>>>
>>> On May 4, 2023, at 11:23 AM, LEONARDO MUTTI <
>>> leonardo.mutti01 at universitadipavia.it> wrote:
>>>
>>> Thank you for the help.
>>> Adding to my example:
>>>
>>>
>>> * call PCGASMSetSubdomains(pc,NSub, subdomains_IS,
>>> inflated_IS,ierr) call
>>> PCGASMDestroySubdomains(NSub,subdomains_IS,inflated_IS,ierr)*
>>> results in:
>>>
>>> * Error LNK2019 unresolved external symbol PCGASMDESTROYSUBDOMAINS
>>> referenced in function ... *
>>>
>>> * Error LNK2019 unresolved external symbol PCGASMSETSUBDOMAINS
>>> referenced in function ... *
>>> I'm not sure if the interfaces are missing or if I have a compilation
>>> problem.
>>> Thank you again.
>>> Best,
>>> Leonardo
>>>
>>> Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <bsmith at petsc.dev>
>>> ha scritto:
>>>
>>>>
>>>> Thank you for the test code. I have a fix in the branch
>>>> barry/2023-04-29/fix-pcasmcreatesubdomains2d
>>>> <https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d> with
>>>> merge request https://gitlab.com/petsc/petsc/-/merge_requests/6394
>>>>
>>>> The functions did not have proper Fortran stubs and interfaces so I
>>>> had to provide them manually in the new branch.
>>>>
>>>> Use
>>>>
>>>> git fetch
>>>> git checkout barry/2023-04-29/fix-pcasmcreatesubdomains2d
>>>> <https://gitlab.com/petsc/petsc/-/commits/barry/2023-04-29/fix-pcasmcreatesubdomains2d>
>>>> ./configure etc
>>>>
>>>> Your now working test code is in src/ksp/ksp/tests/ex71f.F90 I had
>>>> to change things slightly and I updated the error handling for the latest
>>>> version.
>>>>
>>>> Please let us know if you have any later questions.
>>>>
>>>> Barry
>>>>
>>>>
>>>>
>>>>
>>>> On Apr 28, 2023, at 12:07 PM, LEONARDO MUTTI <
>>>> leonardo.mutti01 at universitadipavia.it> wrote:
>>>>
>>>> Hello. I am having a hard time understanding the index sets to feed
>>>> PCGASMSetSubdomains, and I am working in Fortran (as a PETSc novice). To
>>>> get more intuition on how the IS objects behave I tried the following
>>>> minimal (non) working example, which should tile a 16x16 matrix into 16
>>>> square, non-overlapping submatrices:
>>>>
>>>> #include <petsc/finclude/petscmat.h>
>>>> #include <petsc/finclude/petscksp.h>
>>>> #include <petsc/finclude/petscpc.h>
>>>> USE petscmat
>>>> USE petscksp
>>>> USE petscpc
>>>>
>>>> Mat :: A
>>>> PetscInt :: M, NSubx, dof, overlap, NSub
>>>> INTEGER :: I,J
>>>> PetscErrorCode :: ierr
>>>> PetscScalar :: v
>>>> KSP :: ksp
>>>> PC :: pc
>>>> IS :: subdomains_IS, inflated_IS
>>>>
>>>> call PetscInitialize(PETSC_NULL_CHARACTER , ierr)
>>>>
>>>> !-----Create a dummy matrix
>>>> M = 16
>>>> call MatCreateAIJ(MPI_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE,
>>>> & M, M,
>>>> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,
>>>> & PETSC_DEFAULT_INTEGER, PETSC_NULL_INTEGER,
>>>> & A, ierr)
>>>>
>>>> DO I=1,M
>>>> DO J=1,M
>>>> v = I*J
>>>> CALL MatSetValue (A,I-1,J-1,v,
>>>> & INSERT_VALUES , ierr)
>>>> END DO
>>>> END DO
>>>>
>>>> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY , ierr)
>>>> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY , ierr)
>>>>
>>>> !-----Create KSP and PC
>>>> call KSPCreate(PETSC_COMM_WORLD,ksp, ierr)
>>>> call KSPSetOperators(ksp,A,A, ierr)
>>>> call KSPSetType(ksp,"bcgs",ierr)
>>>> call KSPGetPC(ksp,pc,ierr)
>>>> call KSPSetUp(ksp, ierr)
>>>> call PCSetType(pc,PCGASM, ierr)
>>>> call PCSetUp(pc , ierr)
>>>>
>>>> !-----GASM setup
>>>> NSubx = 4
>>>> dof = 1
>>>> overlap = 0
>>>>
>>>> call PCGASMCreateSubdomains2D(pc,
>>>> & M, M,
>>>> & NSubx, NSubx,
>>>> & dof, overlap,
>>>> & NSub, subdomains_IS, inflated_IS, ierr)
>>>>
>>>> call ISView(subdomains_IS, PETSC_VIEWER_STDOUT_WORLD, ierr)
>>>>
>>>> call KSPDestroy(ksp, ierr)
>>>> call PetscFinalize(ierr)
>>>>
>>>> Running this on one processor, I get NSub = 4.
>>>> If PCASM and PCASMCreateSubdomains2D are used instead, I get NSub = 16
>>>> as expected.
>>>> Moreover, I get in the end "forrtl: severe (157): Program Exception -
>>>> access violation". So:
>>>> 1) why do I get two different results with ASM, and GASM?
>>>> 2) why do I get access violation and how can I solve this?
>>>> In fact, in C, subdomains_IS, inflated_IS should pointers to IS
>>>> objects. As I see on the Fortran interface, the arguments to
>>>> PCGASMCreateSubdomains2D are IS objects:
>>>>
>>>> subroutine PCGASMCreateSubdomains2D(a,b,c,d,e,f,g,h,i,j,z)
>>>> import tPC,tIS
>>>> PC a ! PC
>>>> PetscInt b ! PetscInt
>>>> PetscInt c ! PetscInt
>>>> PetscInt d ! PetscInt
>>>> PetscInt e ! PetscInt
>>>> PetscInt f ! PetscInt
>>>> PetscInt g ! PetscInt
>>>> PetscInt h ! PetscInt
>>>> IS i ! IS
>>>> IS j ! IS
>>>> PetscErrorCode z
>>>> end subroutine PCGASMCreateSubdomains2D
>>>> Thus:
>>>> 3) what should be inside e.g., subdomains_IS? I expect it to contain,
>>>> for every created subdomain, the list of rows and columns defining the
>>>> subblock in the matrix, am I right?
>>>>
>>>> Context: I have a block-tridiagonal system arising from space-time
>>>> finite elements, and I want to solve it with GMRES+PCGASM preconditioner,
>>>> where each overlapping submatrix is on the diagonal and of size 3x3 blocks
>>>> (and spanning multiple processes). This is PETSc 3.17.1 on Windows.
>>>>
>>>> Thanks in advance,
>>>> Leonardo
>>>>
>>>>
>>>>
>>>
>
> --
> 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/20230509/98d0ce22/attachment-0001.html>
More information about the petsc-users
mailing list