On Fri, May 5, 2023 at 2:45 AM LEONARDO MUTTI < leonardo.mutt...@universitadipavia.it> 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 <knep...@gmail.com> ha scritto: > >> On Thu, May 4, 2023 at 1:43 PM LEONARDO MUTTI < >> leonardo.mutt...@universitadipavia.it> 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 < >>> knep...@gmail.com> ha scritto: >>> >>>> On Thu, May 4, 2023 at 11:24 AM LEONARDO MUTTI < >>>> leonardo.mutt...@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. >>>>> >>>> >>>> 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 <bsm...@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.mutt...@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/> >>>> >>> >> >> -- >> 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/>