Understanding index sets for PCGASM
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
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 <[email protected]> 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
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 <[email protected]> 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 < [email protected]> 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
Looks like we don't have Fortran interfaces for these. We need to add them. Barry
On May 4, 2023, at 11:23 AM, LEONARDO MUTTI <[email protected]> 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 <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> 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
On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI < [email protected]> 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.
I just want to make sure you really want GASM. It sounded like you might able to do what you want just with ASM. Can you tell me again what you want to do overall? Thanks, Matt
Thank you again. Best, Leonardo
Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <[email protected]> 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 < [email protected]> 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/>
Of course, I'll try to explain. I am solving a parabolic equation with space-time FEM and I want an efficient solver/preconditioner for the resulting system. The corresponding matrix, call it X, has an e.g. block bi-diagonal structure, if the cG(1)-dG(0) method is used (i.e. implicit Euler solved in batch). Every block-row of X corresponds to a time instant. I want to introduce parallelism in time by subdividing X into overlapping submatrices of e.g 2x2 or 3x3 blocks, along the block diagonal. For instance, call X_i the individual blocks. The submatrices would be, for various i, (X_{i-1,i-1},X_{i-1,i};X_{i,i-1},X_{i,i}). I'd like each submatrix to be solved in parallel, to combine the various results together in an ASM like fashion. Every submatrix has thus a predecessor and a successor, and it overlaps with both, so that as far as I could understand, GASM has to be used in place of ASM. Hope this helps. Best, Leonardo Il giorno gio 4 mag 2023 alle ore 18:05 Matthew Knepley <[email protected]> ha scritto:
On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI < [email protected]> 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.
I just want to make sure you really want GASM. It sounded like you might able to do what you want just with ASM. Can you tell me again what you want to do overall?
Thanks,
Matt
Thank you again. Best, Leonardo
Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <[email protected]> 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 < [email protected]> 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/>
On Thu, May 4, 2023 at 1:43 PM LEONARDO MUTTI < [email protected]> wrote:
Of course, I'll try to explain.
I am solving a parabolic equation with space-time FEM and I want an efficient solver/preconditioner for the resulting system. The corresponding matrix, call it X, has an e.g. block bi-diagonal structure, if the cG(1)-dG(0) method is used (i.e. implicit Euler solved in batch). Every block-row of X corresponds to a time instant.
I want to introduce parallelism in time by subdividing X into overlapping submatrices of e.g 2x2 or 3x3 blocks, along the block diagonal. For instance, call X_i the individual blocks. The submatrices would be, for various i, (X_{i-1,i-1},X_{i-1,i};X_{i,i-1},X_{i,i}). I'd like each submatrix to be solved in parallel, to combine the various results together in an ASM like fashion. Every submatrix has thus a predecessor and a successor, and it overlaps with both, so that as far as I could understand, GASM has to be used in place of ASM.
Yes, ordered that way you need GASM. I wonder if inverting the ordering would be useful, namely putting the time index on the inside. Then the blocks would be over all time, but limited space, which is more the spirit of ASM I think. Have you considered waveform relaxation for this problem? Thanks, Matt
Hope this helps. Best, Leonardo
Il giorno gio 4 mag 2023 alle ore 18:05 Matthew Knepley <[email protected]> ha scritto:
On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI < [email protected]> 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.
I just want to make sure you really want GASM. It sounded like you might able to do what you want just with ASM. Can you tell me again what you want to do overall?
Thanks,
Matt
Thank you again. Best, Leonardo
Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <[email protected]> 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 < [email protected]> 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/>
-- 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/>
Interesting, a priori I'm not sure this will work better, mainly because I'd lose the compact band structure. As for waveform relaxation: I excluded it at first since it appears to be requiring too many CPUs than I have to beat sequential solvers, plus it is more complicated and I have very limited time at this project. For both suggestions, because of the way the space-time matrix is generated, it is much more convenient for me to mess with the time dimension than with space. Overall GASM seems a simpler way to go before trying other things. Please let me know if you decide to add the GASM interfaces. Thanks again. Best, Leonardo Il ven 5 mag 2023, 00:01 Matthew Knepley <[email protected]> ha scritto:
On Thu, May 4, 2023 at 1:43 PM LEONARDO MUTTI < [email protected]> wrote:
Of course, I'll try to explain.
I am solving a parabolic equation with space-time FEM and I want an efficient solver/preconditioner for the resulting system. The corresponding matrix, call it X, has an e.g. block bi-diagonal structure, if the cG(1)-dG(0) method is used (i.e. implicit Euler solved in batch). Every block-row of X corresponds to a time instant.
I want to introduce parallelism in time by subdividing X into overlapping submatrices of e.g 2x2 or 3x3 blocks, along the block diagonal. For instance, call X_i the individual blocks. The submatrices would be, for various i, (X_{i-1,i-1},X_{i-1,i};X_{i,i-1},X_{i,i}). I'd like each submatrix to be solved in parallel, to combine the various results together in an ASM like fashion. Every submatrix has thus a predecessor and a successor, and it overlaps with both, so that as far as I could understand, GASM has to be used in place of ASM.
Yes, ordered that way you need GASM. I wonder if inverting the ordering would be useful, namely putting the time index on the inside. Then the blocks would be over all time, but limited space, which is more the spirit of ASM I think.
Have you considered waveform relaxation for this problem?
Thanks,
Matt
Hope this helps. Best, Leonardo
Il giorno gio 4 mag 2023 alle ore 18:05 Matthew Knepley < [email protected]> ha scritto:
On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI < [email protected]> 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.
I just want to make sure you really want GASM. It sounded like you might able to do what you want just with ASM. Can you tell me again what you want to do overall?
Thanks,
Matt
Thank you again. Best, Leonardo
Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <[email protected]> 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 < [email protected]> 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/>
-- 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/>
On Fri, May 5, 2023 at 2:45 AM LEONARDO MUTTI < [email protected]> wrote:
Interesting, a priori I'm not sure this will work better, mainly because I'd lose the compact band structure.
As for waveform relaxation: I excluded it at first since it appears to be requiring too many CPUs than I have to beat sequential solvers, plus it is more complicated and I have very limited time at this project.
For both suggestions, because of the way the space-time matrix is generated, it is much more convenient for me to mess with the time dimension than with space.
Overall GASM seems a simpler way to go before trying other things. Please let me know if you decide to add the GASM interfaces.
We will add them. It might take until the end of classes this term. Thanks, Matt
Thanks again. Best, Leonardo
Il ven 5 mag 2023, 00:01 Matthew Knepley <[email protected]> ha scritto:
On Thu, May 4, 2023 at 1:43 PM LEONARDO MUTTI < [email protected]> wrote:
Of course, I'll try to explain.
I am solving a parabolic equation with space-time FEM and I want an efficient solver/preconditioner for the resulting system. The corresponding matrix, call it X, has an e.g. block bi-diagonal structure, if the cG(1)-dG(0) method is used (i.e. implicit Euler solved in batch). Every block-row of X corresponds to a time instant.
I want to introduce parallelism in time by subdividing X into overlapping submatrices of e.g 2x2 or 3x3 blocks, along the block diagonal. For instance, call X_i the individual blocks. The submatrices would be, for various i, (X_{i-1,i-1},X_{i-1,i};X_{i,i-1},X_{i,i}). I'd like each submatrix to be solved in parallel, to combine the various results together in an ASM like fashion. Every submatrix has thus a predecessor and a successor, and it overlaps with both, so that as far as I could understand, GASM has to be used in place of ASM.
Yes, ordered that way you need GASM. I wonder if inverting the ordering would be useful, namely putting the time index on the inside. Then the blocks would be over all time, but limited space, which is more the spirit of ASM I think.
Have you considered waveform relaxation for this problem?
Thanks,
Matt
Hope this helps. Best, Leonardo
Il giorno gio 4 mag 2023 alle ore 18:05 Matthew Knepley < [email protected]> ha scritto:
On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI < [email protected]> 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.
I just want to make sure you really want GASM. It sounded like you might able to do what you want just with ASM. Can you tell me again what you want to do overall?
Thanks,
Matt
Thank you again. Best, Leonardo
Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <[email protected]> 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 < [email protected]> 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/>
-- 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/>
I will add the two interfaces you requested today. (Likely you may need more also). Barry
On May 4, 2023, at 6:01 PM, Matthew Knepley <[email protected]> wrote:
On Thu, May 4, 2023 at 1:43 PM LEONARDO MUTTI <[email protected] <mailto:[email protected]>> wrote:
Of course, I'll try to explain.
I am solving a parabolic equation with space-time FEM and I want an efficient solver/preconditioner for the resulting system. The corresponding matrix, call it X, has an e.g. block bi-diagonal structure, if the cG(1)-dG(0) method is used (i.e. implicit Euler solved in batch). Every block-row of X corresponds to a time instant.
I want to introduce parallelism in time by subdividing X into overlapping submatrices of e.g 2x2 or 3x3 blocks, along the block diagonal. For instance, call X_i the individual blocks. The submatrices would be, for various i, (X_{i-1,i-1},X_{i-1,i};X_{i,i-1},X_{i,i}). I'd like each submatrix to be solved in parallel, to combine the various results together in an ASM like fashion. Every submatrix has thus a predecessor and a successor, and it overlaps with both, so that as far as I could understand, GASM has to be used in place of ASM.
Yes, ordered that way you need GASM. I wonder if inverting the ordering would be useful, namely putting the time index on the inside. Then the blocks would be over all time, but limited space, which is more the spirit of ASM I think.
Have you considered waveform relaxation for this problem?
Thanks,
Matt
Hope this helps. Best, Leonardo
Il giorno gio 4 mag 2023 alle ore 18:05 Matthew Knepley <[email protected] <mailto:[email protected]>> ha scritto:
On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI <[email protected] <mailto:[email protected]>> 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.
I just want to make sure you really want GASM. It sounded like you might able to do what you want just with ASM. Can you tell me again what you want to do overall?
Thanks,
Matt
Thank you again. Best, Leonardo
Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> 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/>
-- 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/>
Thanks a lot. If this can help, we should need (not much more than) the functionalities from https://petsc.org/release/src/ksp/ksp/tutorials/ex62.c.html: - PCGASMSetSubdomains - PCGASMDestroySubdomains - PCGASMGetSubKSP Best, Leonardo Il giorno ven 5 mag 2023 alle ore 17:00 Barry Smith <[email protected]> ha scritto:
I will add the two interfaces you requested today. (Likely you may need more also).
Barry
On May 4, 2023, at 6:01 PM, Matthew Knepley <[email protected]> wrote:
On Thu, May 4, 2023 at 1:43 PM LEONARDO MUTTI < [email protected]> wrote:
Of course, I'll try to explain.
I am solving a parabolic equation with space-time FEM and I want an efficient solver/preconditioner for the resulting system. The corresponding matrix, call it X, has an e.g. block bi-diagonal structure, if the cG(1)-dG(0) method is used (i.e. implicit Euler solved in batch). Every block-row of X corresponds to a time instant.
I want to introduce parallelism in time by subdividing X into overlapping submatrices of e.g 2x2 or 3x3 blocks, along the block diagonal. For instance, call X_i the individual blocks. The submatrices would be, for various i, (X_{i-1,i-1},X_{i-1,i};X_{i,i-1},X_{i,i}). I'd like each submatrix to be solved in parallel, to combine the various results together in an ASM like fashion. Every submatrix has thus a predecessor and a successor, and it overlaps with both, so that as far as I could understand, GASM has to be used in place of ASM.
Yes, ordered that way you need GASM. I wonder if inverting the ordering would be useful, namely putting the time index on the inside. Then the blocks would be over all time, but limited space, which is more the spirit of ASM I think.
Have you considered waveform relaxation for this problem?
Thanks,
Matt
Hope this helps. Best, Leonardo
Il giorno gio 4 mag 2023 alle ore 18:05 Matthew Knepley < [email protected]> ha scritto:
On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI < [email protected]> 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.
I just want to make sure you really want GASM. It sounded like you might able to do what you want just with ASM. Can you tell me again what you want to do overall?
Thanks,
Matt
Thank you again. Best, Leonardo
Il giorno sab 29 apr 2023 alle ore 20:30 Barry Smith <[email protected]> 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 < [email protected]> 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/>
-- 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/>
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 <[email protected]> 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 <[email protected] <mailto:[email protected]>> 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 <[email protected] <mailto:[email protected]>> 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
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....: in particular each element of the matrix A corresponds to a number from 0 to 63. 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 <[email protected]> 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 < [email protected]> 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 <[email protected]> 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 < [email protected]> 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
On Tue, May 9, 2023 at 10:05 AM LEONARDO MUTTI < [email protected]> 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....: 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 <[email protected]> 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 < [email protected]> 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 <[email protected]> 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 < [email protected]> 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/>
participants (3)
-
Barry Smith -
LEONARDO MUTTI -
Matthew Knepley