Create another matrix B that contains "Selected parts of the matrix can be pre-calculated and stored to serve as the preconditioner." When you call KSPSetOperators() pass in B as the second matrix argument.
Now if literally application of B is the preconditioner you want to use, then simply use PCMAT as the preconditioner, and it will apply B each time it needs to apply the preconditioner. With this, you do not even need a PCShell. If your preconditioner uses B to construct your own specific data structure needed to apply your preconditioner, then you can use PCShellSetSetUp() and have your routine call PCGetOperators() to pull out the B matrix and use it to construct your preconditioner. If you want PETSc to construct a standard preconditioner from B, then you don't need the PCShell; you simply pass in the B, as above, and use -pc_type type to set the preconditioner you want to use. Barry > On Aug 16, 2022, at 12:31 PM, Tu, Jiannan <jiannan...@uml.edu> wrote: > > Barry, > > Thank you very much for your instructions. I am sorry I may not ask clearer > questions. > > I don't use SNES. What I am trying to solve is a very large linear equation > system with dense and ill-conditioned matrix. Selected parts of the matrix > can be pre-calculated and stored to serve as the preconditioner. My question > is how I can apply this matrix as the preconditioner. > > Here is the part of code from my linear equations solver. > MatrixFreePreconditioner is supposed to be the routine to provide the > preconditioner. But I don't understand how to use it. Should I use > SNESSetJacobian() even if this is a linear problem? > > MatCreateMFFD(MPI_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, N, &A); > MatMFFDSetFunction(A, formfunction, ¶ms); > MatSetFromOptions(A); > MatMFFDSetBase(A, X, NULL); > > KSPCreate(MPI_COMM_WORLD, &ksp); > KSPSetOperators(ksp, A, A); > KSPSetFromOptions(ksp); > > PC pc; > KSPGetPC(ksp,&pc); > PCSetType(pc,PCSHELL); > PCShellSetApply(pc,MatrixFreePreconditioner); > > Thank you, > Jiannan > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Tuesday, August 16, 2022 10:10 AM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> > <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> > Subject: Re: [petsc-users] Using matrix-free with KSP > > CAUTION: This email was sent from outside the UMass Lowell network. > > > I don't fully understand your question so I will answer questions that I do > know the answer to :-) > > The function you provide to PCShellSetApply() applies the preconditioner to > the input vector and puts the result in the output vector. This input vector > changes at every application of the preconditioner so should not be used in > constructing the preconditioner. > > If you are using the standard PETSc matrix-free matrix-vector product via > finite difference, -snes_mf_operator or via MatCreateSNESMF() or you are > using a MatShell, the matrix-vector product routines input changes for each > multiply and should not be used in constructing the preconditioner. > > The function you provide with SNESSetJacobian() is where you can compute > anything you need to help build your preconditioner. The input vector to that > function is the location at which the Jacobian is needed and is what you > should use to build your approximate/part of Jacobian. You can store your > approximate/part of Jacobian that you need to compute in the pmat and then in > a PCShellSetUp() function that you provide you will compute what you will > need latter to apply the preconditioner in your PCShellSetApply() function. > Note that the values you put in the pmat do not need to be the full Jacobian; > they can be anything you want since this pmat is not used to apply the matrix > vector product. > > If this is unclear, feel free to ask additional questions. > > Barry > > >> On Aug 16, 2022, at 9:40 AM, Tu, Jiannan <jiannan...@uml.edu >> <mailto:jiannan...@uml.edu>> wrote: >> >> I am trying to apply a user-defined preconditioner to speed up the >> convergency of matrix-free KSP solver. The user's manual provides an example >> of how to call the user-defined preconditioner routine. The routine has PC, >> input Vec and output Vec as arguments. I want to use part of the Jacobian as >> preconditioner matrix, the elements of which can be calculated using input >> Vec. My question is what the role of output Vec from the routine is in >> matrix-vector product? Or I just have the preconditioner matrix elements >> calculated in the routine and then Petsc handles applying the preconditioner >> to matrix-vector product? I used PCSetType(pc, PCSHELL) and PCShellSetApply. >> >> Thank you and appreciate your help! >> >> Jiannan >> From: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> >> Sent: Wednesday, July 6, 2022 2:34 PM >> To: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> >> Cc: Jed Brown <j...@jedbrown.org <mailto:j...@jedbrown.org>>; >> petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> >> <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> >> Subject: Re: [petsc-users] Using matrix-free with KSP >> >> Barry, >> >> VecScatterCreateToAll solves the problem! The code now gives the correct >> answer with multiple processes. >> >> Thank you all so much! >> >> Jiannan >> From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> >> Sent: Wednesday, July 6, 2022 2:08 PM >> To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> >> Cc: Jed Brown <j...@jedbrown.org <mailto:j...@jedbrown.org>>; >> petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> >> <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> >> Subject: Re: [petsc-users] Using matrix-free with KSP >> >> >> So your operator is a "dense" operator in that its matrix representation >> would be essentially dense. In that case, you can use specialized routines >> to do this efficiently. You can use VecScatterCreateToAll() and then >> VecScatterBegin/End to do the communication. >> >> Barry >> >> >> >>> On Jul 6, 2022, at 1:39 PM, Tu, Jiannan <jiannan...@uml.edu >>> <mailto:jiannan...@uml.edu>> wrote: >>> >>> Jed, thank you very much for the reply. >>> >>> I'm not sure GlobaltoLocal will work. Say the solution vector is of length >>> 10. Two processes then each process can see 5 elements of the vector. For A >>> x = b (10 equations), each process performs matrix-vector product for five >>> equations. And for each equation i , sum_j A_ij * x_j = b_i requires the >>> process has access to all 10 unknows x_i (i=0,1,...,9). >>> >>> Can VecScatter and PetscSF make the entire vector accessible to each >>> process? I am reading the manual and trying to understand how VecScatter >>> and PetscSF work. >>> >>> Jiannan >>> From: Jed Brown <j...@jedbrown.org <mailto:j...@jedbrown.org>> >>> Sent: Wednesday, July 6, 2022 10:09 AM >>> To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>>; Barry >>> Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> >>> Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> >>> <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> >>> Subject: Re: [petsc-users] Using matrix-free with KSP >>> >>> You'll usually have a GlobalToLocal operation for each rank to get the halo >>> data it needs, then either a LocalToGlobal to collect the result (e.g., >>> finite element methods) or the local compute will write directly into the >>> owned portion of the global vector. If you're doing the communication >>> "raw", then you may find VecScatter or PetscSF useful to perform the >>> necessary communication. >>> >>> "Tu, Jiannan" <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> writes: >>> >>> > Hi Barry, >>> > >>> > Following your instructions I implemented the matrix-free method to solve >>> > a large linear equation system resulted from a surface integration >>> > equation. The KSP solver works fine for single process, but it is >>> > problematic with multiple processes. The problem is that each process >>> > only can access its own part of solution vector so that each process only >>> > conducts part of matrix-vector multiplication. MPI can be used to >>> > assemble these partial matrix-vector multiplication together. Does Petsc >>> > provide any ways to implement multi-process matrix-free method? >>> > >>> > Thanks, >>> > Jiannan >>> > ________________________________ >>> > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> >>> > Sent: Tuesday, May 24, 2022 2:12 PM >>> > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> >>> > Cc: petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov> >>> > <petsc-users@mcs.anl.gov <mailto:petsc-users@mcs.anl.gov>> >>> > Subject: Re: [petsc-users] Using matrix-free with KSP >>> > >>> > This e-mail originated from outside the UMass Lowell network. >>> > ________________________________ >>> > >>> > You can use MatCreateMFFD >>> > https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Fmain%2Fdocs%2Fmanualpages%2FMat%2FMatCreateMFFD%2F&data=05%7C01%7CJiannan_Tu%40uml.edu%7Cab0859d8d154471590bc08da5f5918d6%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637927133525745809%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=A7wqegTfh94No5BpDiWLK3VxOuR44U2wlWHVm2k7l60%3D&reserved=0<https://urldefense.com/v3/__https://petsc.org/main/docs/manualpages/Mat/MatCreateMFFD/__;!!PVKG_VDCxu5g!stxq2NWGVjUihvHBQC9Rlk0aHVoowcy0PTQYMcQxQ4DqOEC05KSw6TmstSXKckiUgelHzy4ue-d10-zlXkw$> >>> > >>> > <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Fmain%2Fdocs%2Fmanualpages%2FMat%2FMatCreateMFFD%2F&data=05%7C01%7CJiannan_Tu%40uml.edu%7C9cf9b4c72f7c4f0cd92a08da7f9120f4%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637962558562603998%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=BdPTAQcdmpnMYPbWECCabNQSzJfqpiH5G7iq1d%2FlY4k%3D&reserved=0> >>> > MatMFFDSetFunction MatSetFromOptions MatMFFDSetBase and provide the >>> > matrix to KSP. Note you will need to use -pc_type none or provide your >>> > own custom preconditioner >>> > withhttps://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Fmain%2Fdocs%2Fmanualpages%2FPC%2FPCSHELL%2F&data=05%7C01%7CJiannan_Tu%40uml.edu%7Cab0859d8d154471590bc08da5f5918d6%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637927133525745809%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=Gb%2F9uFG3jp%2FbfWaSh2cSEQLItHS9CuLwarzN0KcNiJY%3D&reserved=0<https://urldefense.com/v3/__https://petsc.org/main/docs/manualpages/PC/PCSHELL/__;!!PVKG_VDCxu5g!stxq2NWGVjUihvHBQC9Rlk0aHVoowcy0PTQYMcQxQ4DqOEC05KSw6TmstSXKckiUgelHzy4ue-d1cA4fKGw$> >>> > >>> > <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Fmain%2Fdocs%2Fmanualpages%2FPC%2FPCSHELL%2F&data=05%7C01%7CJiannan_Tu%40uml.edu%7C9cf9b4c72f7c4f0cd92a08da7f9120f4%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637962558562603998%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=6XJ3yp6S8IbTevj%2BeUOHDjgCbeTq346j%2Fini%2BLGmkag%3D&reserved=0> >>> > >>> > >>> > >>> > On May 24, 2022, at 1:21 PM, Tu, Jiannan <jiannan...@uml.edu >>> > <mailto:jiannan...@uml.edu><mailto:jiannan...@uml.edu >>> > <mailto:jiannan...@uml.edu>>> wrote: >>> > >>> > I want to use a matrix-free matrix to solve a large linear equation >>> > system because the matrix is too large to be stored. Petsc user manual >>> > describes matrix-free method for SNES with examples. The matrix-free >>> > matrices section explains how to set up such a matrix, but I can't find >>> > any example of matrix-free method with KSP. I am wondering how to apply >>> > the method to KSP iterative solver in parallel. >>> > >>> > I greatly appreciate your help for me to understand this topic. >>> > >>> > Jiannan Tu
