Hi Mat, Thanks for your reply. It seems I have to use "VecSetValues" to assign the values of my ghost vector "petsc_dcsv". and then call VecAssemblyBegin/End. If I do it this way, the ghost cells are exchanged correctly.
Besides, I notice that, when I run my code sequentially or with multiple processors, the produced eigenvalues are similar, but the number of iterations are different to reach the specified "-eps_tol" and the relative residuals are also slightly different. Is this normal? I am using the default Krylov-Schur solver and double precision. Thanks, Feng ________________________________ From: Matthew Knepley <[email protected]> Sent: 09 October 2022 12:11 To: feng wang <[email protected]> Cc: Jose E. Roman <[email protected]>; [email protected] <[email protected]> Subject: Re: [petsc-users] Slepc, shell matrix, parallel, halo exchange On Fri, Oct 7, 2022 at 5:48 PM feng wang <[email protected]<mailto:[email protected]>> wrote: Hi Mat, I've tried the suggested approach. The halo cells are not exchanged somehow. Below is how I do it, have I missed anything? I create a ghost vector petsc_dcsv and it is a data member of the class cFdDomain, which is a context of the shell matrix. PetscCall(VecCreateGhostBlock(*A_COMM_WORLD, blocksize, blocksize*nlocal, PETSC_DECIDE ,nghost, ighost, &petsc_dcsv)); blocksize and nv have the same value. nlocal is number of local cells and nghost is number of halo cells. ighost contains the ghost cell index. Below is how I compute a matrix-vector product with a shell matrix PetscErrorCode cFdDomain::mymult_slepc(Mat m ,Vec x, Vec y) { void *ctx; cFdDomain *myctx; PetscErrorCode ierr; MatShellGetContext(m, &ctx); myctx = (cFdDomain*)ctx; //matrix-vector product ierr = myctx->myfunc(x, y); CHKERRQ(ierr); ierr = 0; return ierr; } PetscErrorCode cFdDomain::myfunc(Vec in, Vec out) { //some declaration ierr = VecGetArray(petsc_dcsv,&array_g); CHKERRQ(ierr); ierr = VecGetArrayRead(in, &array); CHKERRQ(ierr); //assign in to petsc_dcsv, only local cells for(iv=0; iv<nv; iv++) { for(iq=0; iq<nlocal; iq++) { array_g[iv+nv*iq] = array[iv + nv*iq]; } } ierr = VecRestoreArray(petsc_dcsv,&array_g); CHKERRQ(ierr); ierr = VecRestoreArrayRead(in, &array); CHKERRQ(ierr); //update halo cells? PetscCall(VecGhostUpdateBegin(petsc_dcsv, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecGhostUpdateEnd(petsc_dcsv, INSERT_VALUES, SCATTER_FORWARD)); PetscCall(VecGhostGetLocalForm(petsc_dcsv,&veclocal)); //read in v ierr = VecGetArray(veclocal,&array_ghost); CHKERRQ(ierr); for(iv=0; iv<nv; iv++) { for(iq=0; iq<nlocal; iq++) { jq = ilocal[iq]; dq[iv][jq] = array_ghost[iv + nv*iq]; } for(iq=nlocal; iq<nlocal+nghost; iq++) { jq = ighost_local[iq-nlocal]; dq[iv][jq] = array_ghost[iv + nv*iq]; } } ierr = VecRestoreArray(veclocal,&array_ghost); CHKERRQ(ierr); //some computations PetscCall(VecGhostRestoreLocalForm(petsc_dcsv,&veclocal)); } so I fill the local part of the ghost vector petsc_dcsv for each rank and then call ghost update, and think this will update the halo cells. it seems not doing that. I can only think you are misinterpreting the result. There are many examples, such src/vec/tutorials/ex9.c (and ex9f.F) I would start there and try to change that into the communication you want, since it definitely works. I cannot see a problem with the code snippet above. Thanks, Matt Thanks, Feng ________________________________ From: Matthew Knepley <[email protected]<mailto:[email protected]>> Sent: 21 September 2022 14:36 To: feng wang <[email protected]<mailto:[email protected]>> Cc: Jose E. Roman <[email protected]<mailto:[email protected]>>; [email protected]<mailto:[email protected]> <[email protected]<mailto:[email protected]>> Subject: Re: [petsc-users] Slepc, shell matrix, parallel, halo exchange On Wed, Sep 21, 2022 at 10:35 AM feng wang <[email protected]<mailto:[email protected]>> wrote: Hi Jose, For your 2nd suggestion on halo exchange, I get the idea and roughly know how to do it, but there are some implementation details which I am not quite sure. If I understand it correctly, in MatMult(Mat m ,Vec x, Vec y), Vec x is a normal parallel vector and it does not contain halo values. Suppose I create an auxiliary ghost vector x_g, then I assign the values of x to x_g. The values of the halo for each partition will not be assigned at this stage. But If I call VecGhostUpdateBegin/End(x_g, INSERT_VALUES, SCATTER_FORWARD), this will fill the values of the halo cells of x_g for each partition. Then x_g has local and halo cells assigned correctly and I can use x_g to do my computation. Is this what you mean? Yes Matt Thanks, Feng ________________________________ From: Jose E. Roman <[email protected]<mailto:[email protected]>> Sent: 21 September 2022 13:07 To: feng wang <[email protected]<mailto:[email protected]>> Cc: Matthew Knepley <[email protected]<mailto:[email protected]>>; [email protected]<mailto:[email protected]> <[email protected]<mailto:[email protected]>> Subject: Re: [petsc-users] Slepc, shell matrix, parallel, halo exchange > El 21 sept 2022, a las 14:47, feng wang > <[email protected]<mailto:[email protected]>> escribió: > > Thanks Jose, I will try this and will come back to this thread if I have any > issue. > > Besides, for EPSGetEigenpair, I guess each rank gets its portion of the > eigenvector, and I need to put them together afterwards? Eigenvectors are stored in parallel vectors, which are used in subsequent parallel computation in most applications. If for some reason you need to gather them in a single MPI process you can use e.g. VecScatterCreateToZero() > > Thanks, > Feng > > From: Jose E. Roman <[email protected]<mailto:[email protected]>> > Sent: 21 September 2022 12:34 > To: feng wang <[email protected]<mailto:[email protected]>> > Cc: Matthew Knepley <[email protected]<mailto:[email protected]>>; > [email protected]<mailto:[email protected]> > <[email protected]<mailto:[email protected]>> > Subject: Re: [petsc-users] Slepc, shell matrix, parallel, halo exchange > > If you define the MATOP_CREATE_VECS operation in your shell matrix so that it > creates a ghost vector, then all vectors within EPS will be ghost vectors, > including those that are received as arguments of MatMult(). Not sure if this > will work. > > A simpler solution is that you store a ghost vector in the context of your > shell matrix, and then in MatMult() you receive a regular parallel vector x, > then update the ghost points using the auxiliary ghost vector, do the > computation and store the result in the regular parallel vector y. > > Jose > > > > El 21 sept 2022, a las 14:09, feng wang > > <[email protected]<mailto:[email protected]>> escribió: > > > > Thanks for your reply. > > > > For GMRES, I create a ghost vector and give it to KSPSolve. For Slepc, it > > only takes the shell matrix for EPSSetOperators. Suppose the shell matrix > > of the eigensolver defines MatMult(Mat m ,Vec x, Vec y), how does it know > > Vec x is a ghost vector and how many ghost cells there are? > > > > Thanks, > > Feng > > From: Matthew Knepley <[email protected]<mailto:[email protected]>> > > Sent: 21 September 2022 11:58 > > To: feng wang <[email protected]<mailto:[email protected]>> > > Cc: [email protected]<mailto:[email protected]> > > <[email protected]<mailto:[email protected]>> > > Subject: Re: [petsc-users] Slepc, shell matrix, parallel, halo exchange > > > > On Wed, Sep 21, 2022 at 7:41 AM feng wang > > <[email protected]<mailto:[email protected]>> wrote: > > Hello, > > > > I am using Slepc with a shell matrix. The sequential version seems working > > and now I am trying to make it run in parallel. > > > > The partition of the domain is done, I am not sure how to do the halo > > exchange in the shell matrix in Slepc. I have a parallel version of > > matrix-free GMRES in my code with Petsc. I was using VecCreateGhostBlock to > > create vector with ghost cells, and then used VecGhostUpdateBegin/End for > > the halo exchange in the shell matrix, would this be the same for Slepc? > > > > That will be enough for the MatMult(). You would also have to use a SLEPc > > EPS that only needed MatMult(). > > > > Thanks, > > > > Matt > > > > Thanks, > > Feng > > > > > > > > > > -- > > What most experimenters take for granted before they begin their > > experiments is infinitely more interesting than any results to which their > > experiments lead. > > -- Norbert Wiener > > > > https://www.cse.buffalo.edu/~knepley/ -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/> -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/<http://www.cse.buffalo.edu/~knepley/>
