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.
PetscCall(VecGetArrayRead(ref, (const double **)(&pr))) ; After I used the right VecGet/RestoreArrayWrite(), the code worked. --Junchao Zhang On Fri, Jan 26, 2024 at 8:40 AM Pierre Jolivet <pie...@joliv.et> wrote: > > On 26 Jan 2024, at 3:11 PM, Pierre Jolivet <pie...@joliv.et> wrote: > > > On 26 Jan 2024, at 3:03 PM, mich...@paraffinalia.co.uk wrote: > > On 2024-01-23 18:09, Junchao Zhang wrote: > > Do you have an example to reproduce it? > --Junchao Zhang > > > I have put a minimum example on github: > > https://github.com/mjcarley/petsc-test > > It does seem that the problem occurs if I do not use the PETSc interface > to do a matrix multiplication. > > 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 > > b = A*ref > > 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 > > > You should not use VecGetArrayRead() if you change the Vec, but instead, > VecGetArrayWrite(). > Does that solve the issue? > > > Sorry, I sent the message too fast, you are also missing a couple of calls > to VecRestoreArray[Read,Write](). > 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 > https://petsc.org/release/src/vec/vec/interface/rvector.c.html#line2177 > > Thanks, > Pierre > > Thanks, > Pierre > > 0 KSP Residual norm < 1.e-11 > Linear solve converged due to CONVERGED_ATOL iterations 0 > > >