<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=us-ascii">
<meta name="Generator" content="Microsoft Word 15 (filtered medium)">
<!--[if !mso]><style>v\:* {behavior:url(#default#VML);}
o\:* {behavior:url(#default#VML);}
w\:* {behavior:url(#default#VML);}
.shape {behavior:url(#default#VML);}
</style><![endif]--><style><!--
/* Font Definitions */
@font-face
{font-family:Helvetica;
panose-1:2 11 6 4 2 2 2 2 2 4;}
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:DengXian;
panose-1:2 1 6 0 3 1 1 1 1 1;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:"\@DengXian";
panose-1:2 1 6 0 3 1 1 1 1 1;}
@font-face
{font-family:"Segoe UI";
panose-1:2 11 5 2 4 2 4 2 2 3;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0in;
font-size:11.0pt;
font-family:"Calibri",sans-serif;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
span.apple-converted-space
{mso-style-name:apple-converted-space;}
span.xapple-converted-space
{mso-style-name:xapple-converted-space;}
span.xxxapple-converted-space
{mso-style-name:xxxapple-converted-space;}
span.EmailStyle21
{mso-style-type:personal-reply;
font-family:"Calibri",sans-serif;
color:windowtext;}
.MsoChpDefault
{mso-style-type:export-only;
font-size:10.0pt;}
@page WordSection1
{size:8.5in 11.0in;
margin:1.0in 1.0in 1.0in 1.0in;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]-->
</head>
<body lang="EN-US" link="blue" vlink="purple" style="word-wrap:break-word">
<div class="WordSection1">
<p class="MsoNormal"><span style="font-size:12.0pt">Barry, <o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">Thanks. The method works. Assembly of actual matrix now completes quickly. The convergency is much better for the first iteration. The residual normal is in the order of 10^-10 instead of 10^-2. However, there
is an error message after the first iteration: “Object is in wrong state. Not for unassembled matrix” even though both preallocator and actual matrix have been assembled (see the code snippet below). Could you please give me some hints what’s going wrong?
I use KSP type bcgs and PC type mat<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">Thank you,<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">Jiannan<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">-----------------------------------------------------<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">Solve linear equations ...<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> 0 KSP Residual norm 9.841179147519e-10<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: Object is in wrong state<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: Not for unassembled matrix<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: Petsc Release Version 3.16.6, Mar 30, 2022<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: ./iesolver on a named REXNET-WS4 by jiannantu Thu Aug 18 14:48:43 2022<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[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<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #1 MatMult() at /home/jiannantu/petsc-3.16.6/src/mat/interface/matrix.c:2425<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #2 PCApplyBAorAB() at /home/jiannantu/petsc-3.16.6/src/ksp/pc/interface/precon.c:730<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #3 KSP_PCApplyBAorAB() at /home/jiannantu/petsc-3.16.6/include/petsc/private/kspimpl.h:421<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: Linear equations not solved. Diverged reason: 0<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">#4 KSPSolve_BCGS() at /home/jiannantu/petsc-3.16.6/src/ksp/ksp/impls/bcgs/bcgs.c:87<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #5 KSPSolve_Private() at /home/jiannantu/petsc-3.16.6/src/ksp/ksp/interface/itfunc.c:914<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #6 KSPSolve() at /home/jiannantu/petsc-3.16.6/src/ksp/ksp/interface/itfunc.c:1086<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt">-----------------------------------------------<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatCreate(MPI_COMM_WORLD, &B);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetType(B, MATMPIAIJ);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, N, N);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> Mat preall;<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatCreate(MPI_COMM_WORLD, &preall);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetType(preall, MATPREALLOCATOR);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetSizes(preall, PETSC_DECIDE, PETSC_DECIDE, N, N);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatMPIAIJSetPreallocation(preall, 0, d_nnz, 0, o_nnz);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetUp(preall);<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatSetOption(preall, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">//set values row by row<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatSetValues(preall, 1, &m, 1, nCols, matCols, INSERT_VALUES);<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatAssemblyBegin(preall, MAT_FINAL_ASSEMBLY);<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatAssemblyEnd(preall, MAT_FINAL_ASSEMBLY);<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">//set values for actual matrix row by row<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatSetValues(B, 1, &m, 1, nCols, matCols, INSERT_VALUES);<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);<o:p></o:p></span></p>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> KSPCreate(MPI_COMM_WORLD, &ksp);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> KSPSetOperators(ksp, A, B);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> KSPSetFromOptions(ksp);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> PC pc;<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> KSPGetPC(ksp, &pc);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"> PCSetFromOptions(pc);<o:p></o:p></span></p>
<p class="MsoNormal"><span style="font-size:12.0pt"><o:p> </o:p></span></p>
<div>
<div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in">
<p class="MsoNormal"><b>From:</b> Barry Smith <bsmith@petsc.dev> <br>
<b>Sent:</b> Thursday, August 18, 2022 12:59 PM<br>
<b>To:</b> Tu, Jiannan <Jiannan_Tu@uml.edu><br>
<b>Cc:</b> petsc-users@mcs.anl.gov<br>
<b>Subject:</b> Re: [petsc-users] Using matrix-free with KSP<o:p></o:p></p>
</div>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<p class="MsoNormal" style="line-height:16.0pt;background:#FFCACA"><b><span style="font-size:12.0pt;color:black">CAUTION:</span></b><span style="font-size:12.0pt;color:black"> This email was sent from outside the UMass Lowell network.<o:p></o:p></span></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
<p class="MsoNormal"> Yes, once you have preallocated the real matrix you can destroy the preallocation matrix whose only job is to gather the preallocation information<o:p></o:p></p>
<div>
<p class="MsoNormal"><br>
<br>
<o:p></o:p></p>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<p class="MsoNormal">On Aug 18, 2022, at 12:52 PM, Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>> wrote:<o:p></o:p></p>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Thanks, Barry.</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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?</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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.</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Jiannan</span><o:p></o:p></p>
</div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
<div>
<div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in">
<div>
<p class="MsoNormal"><b>From:</b><span class="apple-converted-space"> </span>Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>><span class="apple-converted-space"> </span><br>
<b>Sent:</b><span class="apple-converted-space"> </span>Thursday, August 18, 2022 11:37 AM<br>
<b>To:</b><span class="apple-converted-space"> </span>Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>><br>
<b>Cc:</b><span class="apple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a><br>
<b>Subject:</b><span class="apple-converted-space"> </span>Re: [petsc-users] Using matrix-free with KSP<o:p></o:p></p>
</div>
</div>
</div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
<div>
<div>
<p class="MsoNormal" style="line-height:16.0pt;background:#FFCACA"><b><span style="font-size:12.0pt;color:black">CAUTION:</span></b><span class="apple-converted-space"><span style="font-size:12.0pt;color:black"> </span></span><span style="font-size:12.0pt;color:black">This
email was sent from outside the UMass Lowell network.</span><o:p></o:p></p>
</div>
</div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
<div>
<p class="MsoNormal"> 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.<o:p></o:p></p>
</div>
<div>
<div>
<p class="MsoNormal"><br>
<br>
<br>
<o:p></o:p></p>
</div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<div>
<p class="MsoNormal">On Aug 18, 2022, at 11:30 AM, Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>> wrote:<o:p></o:p></p>
</div>
</div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Hi Barry,</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">The MATPREALLOCATOR solved problem of slow matrix assembly. But the solver failed because of “Matrix type preallocator does not have a multiply defined”.</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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.</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Thank you,</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Jiannan</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: No support for this operation for this object type</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: Matrix type preallocator does not have a multiply defined</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: See<span class="apple-converted-space"> </span></span><a href="https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Frelease%2Ffaq%2F&data=05%7C01%7CJiannan_Tu%40uml.edu%7C92781c0483d44c98a24508da813b02f8%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964387727402803%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=vfijXLof%2FmAMQ%2F6qysJ9Z%2Fx7bG7hwN7z7da6NHBncFE%3D&reserved=0"><span style="font-size:12.0pt">https://petsc.org/release/faq/</span></a><span class="apple-converted-space"><span style="font-size:12.0pt"> </span></span><span style="font-size:12.0pt">for
trouble shooting.</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: Petsc Release Version 3.16.6, Mar 30, 2022</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: ./iesolver on a named REXNET-WS4 by jiannantu Thu Aug 18 11:21:55 2022</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[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</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #1 MatMult() at /home/jiannantu/petsc-3.16.6/src/mat/interface/matrix.c:2437</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #2 PCApply_Mat() at /home/jiannantu/petsc-3.16.6/src/ksp/pc/impls/mat/pcmat.c:9</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #3 PCApply() at /home/jiannantu/petsc-3.16.6/src/ksp/pc/interface/precon.c:445</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #4 KSP_PCApply() at /home/jiannantu/petsc-3.16.6/include/petsc/private/kspimpl.h:382</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #5 KSPInitialResidual() at /home/jiannantu/petsc-3.16.6/src/ksp/ksp/interface/itres.c:65</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #6 KSPSolve_BCGS() at /home/jiannantu/petsc-3.16.6/src/ksp/ksp/impls/bcgs/bcgs.c:43</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #7 KSPSolve_Private() at /home/jiannantu/petsc-3.16.6/src/ksp/ksp/interface/itfunc.c:914</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">[0]PETSC ERROR: #8 KSPSolve() at /home/jiannantu/petsc-3.16.6/src/ksp/ksp/interface/itfunc.c:1086</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
<div>
<div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in">
<div>
<div>
<p class="MsoNormal"><b>From:</b><span class="apple-converted-space"> </span>Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>><span class="apple-converted-space"> </span><br>
<b>Sent:</b><span class="apple-converted-space"> </span>Thursday, August 18, 2022 8:35 AM<br>
<b>To:</b><span class="apple-converted-space"> </span>Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>><br>
<b>Cc:</b><span class="apple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a><br>
<b>Subject:</b><span class="apple-converted-space"> </span>Re: [petsc-users] Using matrix-free with KSP<o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="line-height:16.0pt;background:#FFCACA"><b><span style="font-size:12.0pt;color:black">CAUTION:</span></b><span class="apple-converted-space"><span style="font-size:12.0pt;color:black"> </span></span><span style="font-size:12.0pt;color:black">This
email was sent from outside the UMass Lowell network.</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<p class="MsoNormal"> Slow assembly is usually due to under preallocation. You can change to using <o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);</span><o:p></o:p></p>
</div>
</div>
</div>
</blockquote>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> to detect if you are under preallocating. See<span class="apple-converted-space"> </span><a href="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%7C92781c0483d44c98a24508da813b02f8%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964387727402803%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=YV7wu69YyMaCNJS8KxY4Klhg1zR928MQDOdIvj6cSz4%3D&reserved=0">https://petsc.org/main/docs/manual/mat/#sec-matsparse</a>.
That section clearly needs some sprucing up. If it is difficult to determine good preallocation values you can use <a href="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%7C92781c0483d44c98a24508da813b02f8%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964387727402803%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=xAB2%2BlqbKiXrY0wKdTYJ72wI%2Bj1TrYpNCscaUTt166s%3D&reserved=0">https://petsc.org/main/docs/manualpages/Mat/MATPREALLOCATOR/</a> to
compute the needed preallocation efficiently.<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> 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 <a href="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%7C92781c0483d44c98a24508da813b02f8%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964387727402803%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=v8ueYRddZ0zw8oVm9OpseHsDcQnnGZOndwVF294v4Mc%3D&reserved=0">https://petsc.org/main/docs/manualpages/Mat/MatGetOwnershipRange.html</a><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> Barry<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-family:"Segoe UI",sans-serif;color:#333333"><br>
<br>
<br>
<br>
</span><o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><br>
<br>
<br>
<br>
<o:p></o:p></p>
</div>
</div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<div>
<div>
<p class="MsoNormal">On Aug 18, 2022, at 12:50 AM, Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>> wrote:<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatCreate(MPI_COMM_WORLD, &B);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetType(B, MATMPIAIJ);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetSizes(B, PETSC_DECIDE, PETSC_DECIDE, N, N);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetFromOptions(B);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> // The number of diagonal and off diagonal non zeros only can be determined run-time.</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatMPIAIJSetPreallocation(B, 0, d_nnz, 0, o_nnz);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetOption(B, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> //here insert values row by row.</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatSetValues(B, 1, &m, nn, nCols, matCols, INSERT_VALUES);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">Could you please tell me what I have done wrong?</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">Thank you,</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal" style="text-indent:10.5pt"><span style="font-size:12.0pt">Jiannan</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> <span class="apple-converted-space"> </span></span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div style="border:none;border-top:solid #E1E1E1 1.0pt;padding:3.0pt 0in 0in 0in">
<div>
<div>
<div>
<p class="MsoNormal"><b>From:</b><span class="apple-converted-space"> </span>Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>><span class="apple-converted-space"> </span><br>
<b>Sent:</b><span class="apple-converted-space"> </span>Tuesday, August 16, 2022 2:00 PM<br>
<b>To:</b><span class="apple-converted-space"> </span>Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>><br>
<b>Cc:</b><span class="apple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a><br>
<b>Subject:</b><span class="apple-converted-space"> </span>Re: [petsc-users] Using matrix-free with KSP<o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="line-height:16.0pt;background:#FFCACA"><b><span style="font-size:12.0pt;color:black">CAUTION:</span></b><span class="apple-converted-space"><span style="font-size:12.0pt;color:black"> </span></span><span style="font-size:12.0pt;color:black">This
email was sent from outside the UMass Lowell network.</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> Create another matrix B that contains "<span style="font-size:12.0pt">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.</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> 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. </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> 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. </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> 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. </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> Barry</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><br>
<br>
<br>
<br>
<br>
<o:p></o:p></p>
</div>
</div>
</div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<div>
<div>
<div>
<p class="MsoNormal">On Aug 16, 2022, at 12:31 PM, Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>> wrote:<o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"> <o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Barry,</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Thank you very much for your instructions. I am sorry I may not ask clearer questions.</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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.<span class="apple-converted-space"> </span></span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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?<span class="apple-converted-space"> </span></span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatCreateMFFD(MPI_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, N, N, &A);</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatMFFDSetFunction(A, formfunction, ¶ms);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatSetFromOptions(A);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> MatMFFDSetBase(A, X, NULL);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> KSPCreate(MPI_COMM_WORLD, &ksp);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> KSPSetOperators(ksp, A, A);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> KSPSetFromOptions(ksp);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> PC pc;</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> KSPGetPC(ksp,&pc);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> PCSetType(pc,PCSHELL);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> PCShellSetApply(pc,MatrixFreePreconditioner);</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Thank you,</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Jiannan</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div class="MsoNormal" align="center" style="text-align:center">
<hr size="2" width="1025" style="width:768.8pt" align="center">
</div>
<div id="divRplyFwdMsg">
<div>
<div>
<div>
<p class="MsoNormal"><b>From:</b><span class="apple-converted-space"> </span>Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>><br>
<b>Sent:</b><span class="apple-converted-space"> </span>Tuesday, August 16, 2022 10:10 AM<br>
<b>To:</b><span class="apple-converted-space"> </span>Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>><br>
<b>Cc:</b><span class="apple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a><span class="apple-converted-space"> </span><<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b><span class="apple-converted-space"> </span>Re: [petsc-users] Using matrix-free with KSP<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal" style="line-height:16.0pt;background:#FFCACA"><b><span style="font-size:12.0pt;color:black">CAUTION:</span></b><span class="apple-converted-space"><span style="font-size:12.0pt;color:black"> </span></span><span style="font-size:12.0pt;color:black">This
email was sent from outside the UMass Lowell network.</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> I don't fully understand your question so I will answer questions that I do know the answer to :-)</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> 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.</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> 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.</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> 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.</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> If this is unclear, feel free to ask additional questions.</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> Barry</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"><br>
<br>
<br>
<br>
<br>
</span><o:p></o:p></p>
</div>
</div>
</div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">On Aug 16, 2022, at 9:40 AM, Tu, Jiannan <</span><a href="mailto:Jiannan_Tu@uml.edu"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">Jiannan_Tu@uml.edu</span></a><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">>
wrote:</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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.</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Thank you and appreciate your help!</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Jiannan<span class="xapple-converted-space"> </span></span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div class="MsoNormal" align="center" style="text-align:center"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">
<hr size="2" width="1025" style="width:768.8pt" align="center">
</span></div>
<div id="x_divRplyFwdMsg">
<div>
<div>
<div>
<p class="MsoNormal"><b>From:</b><span class="xapple-converted-space"> </span>Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>><br>
<b>Sent:</b><span class="xapple-converted-space"> </span>Wednesday, July 6, 2022 2:34 PM<br>
<b>To:</b><span class="xapple-converted-space"> </span>Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>><br>
<b>Cc:</b><span class="xapple-converted-space"> </span>Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>>;<span class="xapple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a><span class="xapple-converted-space"> </span><<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b><span class="xapple-converted-space"> </span>Re: [petsc-users] Using matrix-free with KSP<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Barry,</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">VecScatterCreateToAll solves the problem! The code now gives the correct answer with multiple processes.<span class="xapple-converted-space"> </span></span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Thank you all so much!</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Jiannan</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div class="MsoNormal" align="center" style="text-align:center"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">
<hr size="2" width="1025" style="width:768.8pt" align="center">
</span></div>
<div id="x_x_divRplyFwdMsg">
<div>
<div>
<div>
<p class="MsoNormal"><b>From:</b><span class="xapple-converted-space"> </span>Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>><br>
<b>Sent:</b><span class="xapple-converted-space"> </span>Wednesday, July 6, 2022 2:08 PM<br>
<b>To:</b><span class="xapple-converted-space"> </span>Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>><br>
<b>Cc:</b><span class="xapple-converted-space"> </span>Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>>;<span class="xapple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a><span class="xapple-converted-space"> </span><<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b><span class="xapple-converted-space"> </span>Re: [petsc-users] Using matrix-free with KSP<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> 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.</span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> Barry</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"><br>
<br>
<br>
<br>
<br>
</span><o:p></o:p></p>
</div>
</div>
</div>
<blockquote style="margin-top:5.0pt;margin-bottom:5.0pt">
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">On Jul 6, 2022, at 1:39 PM, Tu, Jiannan <</span><a href="mailto:Jiannan_Tu@uml.edu"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">Jiannan_Tu@uml.edu</span></a><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">>
wrote:</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Jed, thank you very much for the reply.<span class="xxxapple-converted-space"> </span></span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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).<span class="xxxapple-converted-space"> </span></span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">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.</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:12.0pt">Jiannan</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
<div class="MsoNormal" align="center" style="text-align:center"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif">
<hr size="2" width="1262" style="width:946.7pt" align="center">
</span></div>
<div id="x_x_x_divRplyFwdMsg">
<div>
<div>
<div>
<p class="MsoNormal"><b>From:</b><span class="xxxapple-converted-space"> </span>Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>><br>
<b>Sent:</b><span class="xxxapple-converted-space"> </span>Wednesday, July 6, 2022 10:09 AM<br>
<b>To:</b><span class="xxxapple-converted-space"> </span>Tu, Jiannan <<a href="mailto:Jiannan_Tu@uml.edu">Jiannan_Tu@uml.edu</a>>; Barry Smith <<a href="mailto:bsmith@petsc.dev">bsmith@petsc.dev</a>><br>
<b>Cc:</b><span class="xxxapple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a><span class="xxxapple-converted-space"> </span><<a href="mailto:petsc-users@mcs.anl.gov">petsc-users@mcs.anl.gov</a>><br>
<b>Subject:</b><span class="xxxapple-converted-space"> </span>Re: [petsc-users] Using matrix-free with KSP<o:p></o:p></p>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-size:13.5pt;font-family:"Helvetica",sans-serif"> </span><o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
<div>
<div>
<div>
<div>
<div>
<p class="MsoNormal"><span style="font-family:"Helvetica",sans-serif">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.<br>
<br>
"Tu, Jiannan" <</span><a href="mailto:Jiannan_Tu@uml.edu"><span style="font-family:"Helvetica",sans-serif">Jiannan_Tu@uml.edu</span></a><span style="font-family:"Helvetica",sans-serif">> writes:<br>
<br>
> Hi Barry,<br>
><br>
> 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?<br>
><br>
> Thanks,<br>
> Jiannan<br>
> ________________________________<br>
> From: Barry Smith <</span><a href="mailto:bsmith@petsc.dev"><span style="font-family:"Helvetica",sans-serif">bsmith@petsc.dev</span></a><span style="font-family:"Helvetica",sans-serif">><br>
> Sent: Tuesday, May 24, 2022 2:12 PM<br>
> To: Tu, Jiannan <</span><a href="mailto:Jiannan_Tu@uml.edu"><span style="font-family:"Helvetica",sans-serif">Jiannan_Tu@uml.edu</span></a><span style="font-family:"Helvetica",sans-serif">><br>
> Cc:<span class="xxxapple-converted-space"> </span></span><a href="mailto:petsc-users@mcs.anl.gov"><span style="font-family:"Helvetica",sans-serif">petsc-users@mcs.anl.gov</span></a><span class="xxxapple-converted-space"><span style="font-family:"Helvetica",sans-serif"> </span></span><span style="font-family:"Helvetica",sans-serif"><</span><a href="mailto:petsc-users@mcs.anl.gov"><span style="font-family:"Helvetica",sans-serif">petsc-users@mcs.anl.gov</span></a><span style="font-family:"Helvetica",sans-serif">><br>
> Subject: Re: [petsc-users] Using matrix-free with KSP<br>
><br>
> This e-mail originated from outside the UMass Lowell network.<br>
> ________________________________<br>
><br>
> You can use MatCreateMFFD<span class="xxxapple-converted-space"> </span></span><a href="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%7C92781c0483d44c98a24508da813b02f8%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964387727402803%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=JTSUN1mODwX3H7TTnT2QlNnz01QNs1GXMtj6BuzH1b8%3D&reserved=0"><span style="font-family:"Helvetica",sans-serif">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$></span></a><span class="xxxapple-converted-space"><span style="font-family:"Helvetica",sans-serif"> </span></span><span style="font-family:"Helvetica",sans-serif">MatMFFDSetFunction
MatSetFromOptions MatMFFDSetBase and provide the matrix to KSP. Note you will need to use -pc_type none or provide your own custom preconditioner with</span><a href="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%7C92781c0483d44c98a24508da813b02f8%7C4c25b8a617f746f983f054734ab81fb1%7C0%7C0%7C637964387727559059%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=hVqa56fqRpHKVHfb9nzD7deTWPOTrTOSSYHT59%2FsyVc%3D&reserved=0"><span style="font-family:"Helvetica",sans-serif">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%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$></span></a><span style="font-family:"Helvetica",sans-serif"><br>
><br>
><br>
><br>
> On May 24, 2022, at 1:21 PM, Tu, Jiannan <</span><a href="mailto:Jiannan_Tu@uml.edu"><span style="font-family:"Helvetica",sans-serif">Jiannan_Tu@uml.edu</span></a><span style="font-family:"Helvetica",sans-serif"><</span><a href="mailto:Jiannan_Tu@uml.edu"><span style="font-family:"Helvetica",sans-serif">mailto:Jiannan_Tu@uml.edu</span></a><span style="font-family:"Helvetica",sans-serif">>>
wrote:<br>
><br>
> 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.<br>
><br>
> I greatly appreciate your help for me to understand this topic.<br>
><br>
> Jiannan Tu</span><o:p></o:p></p>
</div>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</div>
</div>
</blockquote>
</div>
</div>
</div>
</blockquote>
</div>
<p class="MsoNormal"><o:p> </o:p></p>
</div>
</div>
</body>
</html>