Yes, once you have preallocated the real matrix you can destroy the preallocation matrix whose only job is to gather the preallocation information
> On Aug 18, 2022, at 12:52 PM, Tu, Jiannan <jiannan...@uml.edu> wrote: > > Thanks, Barry. > > So I need actually two matrixes, one for preallocator and one for actual > matrix that can be passed to KSPSetOperators(). -mat_type preallocator option > is used to speed up doing insertion into preallocator, then use > MatPreallocatorPreallocate() to preallocate actual matrix and do actual > insertion of values into it., right? > > The code runs in parallel. Each process owns number of rows that equals to > number of unknowns (that is, xm in 1D DM) it owns. > > Jiannan > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Thursday, August 18, 2022 11:37 AM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Cc: 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. > > > The preallocator MatType matrix cannot be passed to the KSPSetOperators(), > you need to create an actual matrix, for example MATAIJ and use the > preallocator to set its preallocation and then fill up the MATAIJ matrix with > the usual MatSetValues and pass that matrix to the KSPSetOperators. > > > On Aug 18, 2022, at 11:30 AM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > Hi Barry, > > The MATPREALLOCATOR solved problem of slow matrix assembly. But the solver > failed because of “Matrix type preallocator does not have a multiply defined”. > > I guess this is because the matrix-free is used. I am wondering how a > preconditioner is applied to the matrix-vector multiply in Petsc. > > Thank you, > Jiannan > > [0]PETSC ERROR: No support for this operation for this object type > [0]PETSC ERROR: Matrix type preallocator does not have a multiply defined > [0]PETSC ERROR: See https://petsc.org/release/faq/ > <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Frelease%2Ffaq%2F&data=05%7C01%7CJiannan_Tu%40uml.edu%7C59dae21b58e64ed9a89108da812f8660%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964338374493241%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=n2BHnbiXekFhF%2FfhzCBMrOhsMLMiHdBzY4eknlnqJi0%3D&reserved=0> > for trouble shooting. > [0]PETSC ERROR: Petsc Release Version 3.16.6, Mar 30, 2022 > [0]PETSC ERROR: ./iesolver on a named REXNET-WS4 by jiannantu Thu Aug 18 > 11:21:55 2022 > [0]PETSC ERROR: Configure options --download-f2cblaslapack=yes > --with-mpi-dir=/usr/local --with-fc=0 --prefix=/home/jiannantu/petsc > --with-scalar-type=complex --with-64-bit-indices=1 > [0]PETSC ERROR: #1 MatMult() at > /home/jiannantu/petsc-3.16.6/src/mat/interface/matrix.c:2437 > [0]PETSC ERROR: #2 PCApply_Mat() at > /home/jiannantu/petsc-3.16.6/src/ksp/pc/impls/mat/pcmat.c:9 > [0]PETSC ERROR: #3 PCApply() at > /home/jiannantu/petsc-3.16.6/src/ksp/pc/interface/precon.c:445 > [0]PETSC ERROR: #4 KSP_PCApply() at > /home/jiannantu/petsc-3.16.6/include/petsc/private/kspimpl.h:382 > [0]PETSC ERROR: #5 KSPInitialResidual() at > /home/jiannantu/petsc-3.16.6/src/ksp/ksp/interface/itres.c:65 > [0]PETSC ERROR: #6 KSPSolve_BCGS() at > /home/jiannantu/petsc-3.16.6/src/ksp/ksp/impls/bcgs/bcgs.c:43 > [0]PETSC ERROR: #7 KSPSolve_Private() at > /home/jiannantu/petsc-3.16.6/src/ksp/ksp/interface/itfunc.c:914 > [0]PETSC ERROR: #8 KSPSolve() at > /home/jiannantu/petsc-3.16.6/src/ksp/ksp/interface/itfunc.c:1086 > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Thursday, August 18, 2022 8:35 AM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Cc: 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. > > > Slow assembly is usually due to under preallocation. You can change to > using > > MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); > > to detect if you are under preallocating. See > https://petsc.org/main/docs/manual/mat/#sec-matsparse > <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Fmain%2Fdocs%2Fmanual%2Fmat%2F%23sec-matsparse&data=05%7C01%7CJiannan_Tu%40uml.edu%7C59dae21b58e64ed9a89108da812f8660%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964338374493241%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=ZYbqedJGTeexabiVb6cLro4wb4tl9GMLrOsDsqtXNNk%3D&reserved=0>. > That section clearly needs some sprucing up. If it is difficult to determine > good preallocation values you can use > https://petsc.org/main/docs/manualpages/Mat/MATPREALLOCATOR/ > <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Fmain%2Fdocs%2Fmanualpages%2FMat%2FMATPREALLOCATOR%2F&data=05%7C01%7CJiannan_Tu%40uml.edu%7C59dae21b58e64ed9a89108da812f8660%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964338374493241%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=3wlhpm2mV1T4m%2B7EOqKhJAMUxE%2F0mkAx%2B8F%2F47oaNVw%3D&reserved=0> > to compute the needed preallocation efficiently. > > Are you running in parallel? If so how are you determining which rows of > entries to compute on each MPI rank? You will want most of the rows to be > computed on the rank where the values are stored, you can determine this with > https://petsc.org/main/docs/manualpages/Mat/MatGetOwnershipRange.html > <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Fmain%2Fdocs%2Fmanualpages%2FMat%2FMatGetOwnershipRange.html&data=05%7C01%7CJiannan_Tu%40uml.edu%7C59dae21b58e64ed9a89108da812f8660%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964338374493241%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=bviDcY78U9cAT2dvLADU49yCivArb4K%2F4Og2bAGFCRg%3D&reserved=0> > > Barry > > > > > > > > > > On Aug 18, 2022, at 12:50 AM, Tu, Jiannan <jiannan...@uml.edu > <mailto:jiannan...@uml.edu>> wrote: > > I implemented the preconditioner. The solver converges quickly when the > problem size is small. But when the size increases, say to 100k unknowns, the > code hangs at assembly of the preconditioner matrix. The function call > sequence is like this > > MatCreate(MPI_COMM_WORLD, &B); > MatSetType(B, MATMPIAIJ); > MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, N, N); > MatSetFromOptions(B); > > // The number of diagonal and off diagonal non zeros only can be > determined run-time. > MatMPIAIJSetPreallocation(B, 0, d_nnz, 0, o_nnz); > MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); > > //here insert values row by row. > MatSetValues(B, 1, &m, nn, nCols, matCols, INSERT_VALUES); > > MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY); > MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY); > > Could you please tell me what I have done wrong? > > Thank you, > Jiannan > > From: Barry Smith <bsm...@petsc.dev <mailto:bsm...@petsc.dev>> > Sent: Tuesday, August 16, 2022 2:00 PM > To: Tu, Jiannan <jiannan...@uml.edu <mailto:jiannan...@uml.edu>> > Cc: 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. > > > 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 > <mailto: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%7C59dae21b58e64ed9a89108da812f8660%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964338374493241%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=2aLlxTwp4PnIfLRef2P8X60XEYgLo1FvhHp13QmXT7E%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%7C59dae21b58e64ed9a89108da812f8660%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964338374493241%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=qPnmCge3WYNb0OIQZMDODm8MpMF9AQliMPlWCjyDjDE%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
