<div dir="ltr"><div>Yes, VecRestoreArray[Read,Write] should be called after you finish using the array.   This casting to (const double **) is bad. Otherwise, compiler could catch the error when you use a const pointer in a non-const way. </div><div><br></div><blockquote style="margin:0 0 0 40px;border:none;padding:0px"><div>PetscCall(VecGetArrayRead(ref, (const double **)(&pr))) ;</div></blockquote><div><br></div><div><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div>After I used the right VecGet/RestoreArrayWrite(), the code worked. </div><div dir="ltr"><br></div><div dir="ltr">--Junchao Zhang<br></div></div></div><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Jan 26, 2024 at 8:40 AM Pierre Jolivet <<a href="mailto:pierre@joliv.et">pierre@joliv.et</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><br><div><blockquote type="cite"><div>On 26 Jan 2024, at 3:11 PM, Pierre Jolivet <<a href="mailto:pierre@joliv.et" target="_blank">pierre@joliv.et</a>> wrote:</div><br><div><blockquote type="cite" style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><br>On 26 Jan 2024, at 3:03 PM, <a href="mailto:michael@paraffinalia.co.uk" target="_blank">michael@paraffinalia.co.uk</a> wrote:<br><br>On 2024-01-23 18:09, Junchao Zhang wrote:<br><blockquote type="cite">Do you have an example to reproduce it?<br>--Junchao Zhang<br></blockquote><br>I have put a minimum example on github:<br><br><a href="https://github.com/mjcarley/petsc-test" target="_blank">https://github.com/mjcarley/petsc-test</a><br><br>It does seem that the problem occurs if I do not use the PETSc interface to do a matrix multiplication.<br><br>In the original code, the PETSc matrix is a wrapper for a Fast Multipole Method evaluation; in the minimum example I have simulated this by using an array as a matrix. The sample code generates a randomised matrix A and reference solution vector ref, and generates a right hand side<br><br>b = A*ref<br><br>which is then supplied as the right hand side for the GMRES solver. If I use the PETSc matrix multiplication, the solver behaves as expected; if I generate b directly from the underlying array for the matrix, I get the result<br></blockquote><br style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><span style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;float:none;display:inline">You should not use VecGetArrayRead() if you change the Vec, but instead, VecGetArrayWrite().</span><br style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><span style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;float:none;display:inline">Does that solve the issue?</span><br style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"></div></blockquote><div><br></div><div>Sorry, I sent the message too fast, you are also missing a couple of calls to VecRestoreArray[Read,Write]().</div><div>These are the ones which will let PETSc know that the Vec has had its state being increased, and that the cached norm are not valid anymore, see <a href="https://petsc.org/release/src/vec/vec/interface/rvector.c.html#line2177" target="_blank">https://petsc.org/release/src/vec/vec/interface/rvector.c.html#line2177</a></div><div><br></div><div>Thanks,</div><div>Pierre</div><br><blockquote type="cite"><div><span style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;float:none;display:inline">Thanks,</span><br style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><span style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none;float:none;display:inline">Pierre</span><br style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><br style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none"><blockquote type="cite" style="font-family:Helvetica;font-size:12px;font-style:normal;font-variant-caps:normal;font-weight:400;letter-spacing:normal;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px;text-decoration:none">0 KSP Residual norm < 1.e-11<br>Linear solve converged due to CONVERGED_ATOL iterations 0</blockquote></div></blockquote></div><br></div></blockquote></div>