[petsc-users] Petsc using VecDuplicate in solution process
Barry Smith
bsmith at petsc.dev
Wed Jun 7 08:10:01 CDT 2023
Can you please present the all output that callgrind is outputing to you that provides this indication.
> On Jun 7, 2023, at 4:11 AM, Pichler, Franz <franz.pichler at v2c2.at> wrote:
>
> Hello thanks for the reply,
>
> I created a working minimal example (as minimal as I can think of?) that I include here, even though I am not sure which is the best format to do this, I just add some plain text:
> //##########################################################################################
> #include <petscksp.h>
> #include <petsc/petscpc.h>
> #include <petscksp.h>
> #include <petscpc.h>
> #include <petscsys.h>
> #include <petscsystypes.h>
> #include <petsc/private/kspimpl.h> /*I "petscksp.h" I*/
> #include <cmath>
> #include <vector>
> #include <iomanip>
> #include <iostream>
>
> class petsc_solver{
> Vec petsc_x, petsc_b;
> Mat petsc_A;
> KSP petsc_ksp;
> PC petsc_pc;
> int linear_iter;
> KSPConvergedReason reason;
> bool first_time;
> int n_rows;
> int number_of_pc_rebuilds=0;
> public:
> petsc_solver() {
> KSPCreate(PETSC_COMM_WORLD, &petsc_ksp);
> KSPSetFromOptions(petsc_ksp);
> KSPGetPC(petsc_ksp, &petsc_pc);
> KSPSetType(petsc_ksp, KSPBCGS);
> PCSetType(petsc_pc, PCILU);
> PCFactorSetLevels(petsc_pc, 0);
> KSPSetInitialGuessNonzero(petsc_ksp, PETSC_TRUE);
> KSPSetTolerances(petsc_ksp, 1.e-12, 1.e-8, 1e14,1000);
> }
> void set_matrix(std::vector<int>& dsp,std::vector<int>& col,std::vector<double>& val){
> int *ptr_dsp = dsp.data();
> int *ptr_col = col.data();
> double *ptr_ele = val.data();
> n_rows=dsp.size()-1;
> std::cout<<"set petsc matrix, n_rows:"<<n_rows<<"\n";
> MatCreateSeqAIJWithArrays(PETSC_COMM_WORLD, n_rows,n_rows, ptr_dsp, ptr_col, ptr_ele,&petsc_A);
> MatSetOption(petsc_A, MAT_NEW_NONZERO_LOCATIONS, PETSC_FALSE);
> MatSetOption(petsc_A, MAT_NO_OFF_PROC_ZERO_ROWS, PETSC_TRUE);
> KSPSetOperators(petsc_ksp, petsc_A, petsc_A);
> }
> void set_rhs(std::vector<double>& val){
> double *ptr_ele = val.data();
> VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, n_rows, NULL,&petsc_b);
> VecPlaceArray(petsc_b, ptr_ele);
> }
> void set_sol(std::vector<double>& val){
> double *ptr_ele = val.data();
> VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, n_rows, NULL,&petsc_x);
> VecPlaceArray(petsc_x, ptr_ele);
> }
>
> int solve(bool force_rebuild) {
> int solver_stat = 0;
> KSPGetPC(petsc_ksp, &petsc_pc);
> int ierr;
> // ierr = PetscObjectStateIncrease((PetscObject)petsc_A);
> // ierr = PetscObjectStateIncrease((PetscObject)petsc_b);
> MatAssemblyBegin(petsc_A,MAT_FINAL_ASSEMBLY);
> MatAssemblyEnd(petsc_A,MAT_FINAL_ASSEMBLY);
> VecAssemblyBegin(petsc_b);
> VecAssemblyEnd(petsc_b);
>
> // KSPSetOperators(petsc_ksp, petsc_A, petsc_A);
> ierr = KSPSolve(petsc_ksp, petsc_b, petsc_x);
> KSPGetConvergedReason(petsc_ksp, &reason);
> if (reason < 0){
> KSPGetIterationNumber(petsc_ksp, &linear_iter);
> std::cout<<"NOT CONVERGED!\n";
> // PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason _aux: %D PCREUSE: %D (%D=False %D=True) IERR:%i ITERS:%i\n",reason, reuse, PETSC_FALSE, PETSC_TRUE, ierr,linear_iter);
> return -1;
> }
> KSPGetIterationNumber(petsc_ksp, &linear_iter);
> return linear_iter;
> }
> };
> void change_rhs(int i, int n_rows,std::vector<double>&rhs){
> for(int row=0;row<n_rows;row++)rhs[row]=sin(1.*i+row)+1.*i*exp(row*1.e-3);//pseduo random something
> }
> void change_matrix(int i, int n_rows,std::vector<double>& vals){
> int nnz = n_rows*3-2;
> for(int row=0;row<n_rows;row++){
> if(row==0) {
> vals[0]=3+cos(i+row);//pseduo random something
> vals[1]=-1+0.3*cos(i+row);//pseduo random something
> }else if(row==n_rows-1){
> vals[nnz-1]=3+cos(i+row);//pseduo random something
> vals[nnz-2]=-1+0.3*cos(i+row);//pseduo random something
> }else{
> vals[2+(row-1)*3] =-1+0.1*cos(i+row);//pseduo random something
> vals[2+(row-1)*3+1] = 4+0.3*cos(i+row);//pseduo random something
> vals[3+(row-1)*3+2] =-1+0.2*cos(i+row);//pseduo random something
> }
> }
> }
> void set_crs_structure(int n_rows,std::vector<int>& dsp,std::vector<int>& col,std::vector<double>& val ){
> int nnz = n_rows*3-2;
> std::cout<<"SETCRS ROWS:"<<n_rows<<" NNZ:"<<nnz<<"\n";
> dsp.resize(n_rows+1);
> col.resize(nnz);
> val.resize(nnz);
> for(int row=0;row<n_rows;row++){
> if(row==0) {
> col[0]=0;
> col[1]=1;
> dsp[row+1]=dsp[row]+2;
> }else if(row==n_rows-1){
> col[2+(row-1)*3+0]=row-1;
> col[2+(row-1)*3+1]=row;
> dsp[row+1]=dsp[row]+2;
> }
> else{
> dsp[row+1]=dsp[row]+3;
> col[2+(row-1)*3+0]=row-1;
> col[2+(row-1)*3+1]=row;
> col[2+(row-1)*3+2]=row+1;
> }
> }
> }
> double check_res(std::vector<int>& dsp,std::vector<int>& col,std::vector<double>& val ,std::vector<double>& sol,std::vector<double> rhs){
> int n_rows=dsp.size()-1;
> double res=0;
> for (int row=0;row<n_rows;row++){
> for(int entry=dsp[row];entry<dsp[row+1];entry++){
> int c=col[entry];
> rhs[row]-=val[entry]*sol[c];
> }
> res+=rhs[row]*rhs[row];
> }
> return sqrt(res);
> }
> int main(int argc, char **argv) {
>
> int n_rows = 20;
> std::vector<double> rhs(n_rows);
> std::vector<double> sol(n_rows);
> std::vector<int> dsp;
> std::vector<int> cols;
> std::vector<double> vals;
> set_crs_structure(n_rows,dsp,cols,vals);
> PetscInitializeNoArguments();
> petsc_solver p;
> p.set_matrix(dsp,cols,vals);
> p.set_rhs(rhs);
> p.set_sol(sol);
> for (int i=0;i<100;i++){
> change_rhs(i,n_rows,rhs);
> change_matrix(i,n_rows,vals);
> // std::cout<<"RES BEFORE:"<<check_res(dsp,cols,vals,sol,rhs)<<"\n";
> int iter = p.solve(false);
> std::cout<<"SOL:"<<i<<" ITERS:"<<iter<<" RES:"<<check_res(dsp,cols,vals,sol,rhs)<<"\n";
> }
> PetscFinalize();
> return -1;
> }
> //##########################################################################################
>
> This is a full working minimal example
> When I callgrind this, it shows me that the vecduplicate is called as often as the solve process itself,
> I hope this clarifies my issue,
>
> Best regards,
>
> Franz
>
> From: Stefano Zampini <stefano.zampini at gmail.com <mailto:stefano.zampini at gmail.com>>
> Sent: Tuesday, June 6, 2023 4:40 PM
> To: Pichler, Franz <franz.pichler at v2c2.at <mailto:franz.pichler at v2c2.at>>
> Cc: petsc-users at mcs.anl.gov <mailto:petsc-users at mcs.anl.gov>
> Subject: Re: [petsc-users] Petsc using VecDuplicate in solution process
>
>
>
> Il giorno mar 6 giu 2023 alle ore 09:24 Pichler, Franz <franz.pichler at v2c2.at <mailto:franz.pichler at v2c2.at>> ha scritto:
> Hello,
> I was just investigating my KSP_Solve_BCGS Routine with algrandcallgrind,
> I see there that petsc is using a vecduplicate (onvolvin malloc and copy) every time it is called.
>
> Do you mean KSPSolve_BCGS?
>
> There's only one VecDuplicate in there and it is called only once. An example code showing the problem would help
>
>
>
> I call it several thousand times (time evolution problem with rather small matrices)
>
> I am not quite sure which vector is copied there but I guess is the initial guess or the rhs,
> Is there a tool in petsc to avoid any vecduplication by providing a fixed memory for this vector?
>
> Some corner facts of my routine:
> I assemble the matrices(crs,serial) and vectors myself and then use
> MatCreateSeqAIJWithArrays and VecCreateSeqWithArray
> To wrap petsc around it,
>
> I use a ILU preconditioner and the sparsity patterns between the calls to not change, the values do,
>
> Thank you for any hint how to avoid the vecduplicate,
>
> Best regards
>
> Franz
>
>
> Dr. Franz Pichler
> Lead Researcher Area E
>
>
> Virtual Vehicle Research GmbH
>
> Inffeldgasse 21a, 8010 Graz, Austria
> Phone: +43 316 873 9818
> franz.pichler at v2c2.at <mailto:franz.pichler at v2c2.at>
> www.v2c2.at <http://www.v2c2.at/>
>
> Firmenname: Virtual Vehicle Research GmbH
> Rechtsform: Gesellschaft mit beschränkter Haftung
> Firmensitz: Austria, 8010 Graz, Inffeldgasse 21/A
> Firmenbuchnummer: 224755y
> Firmenbuchgericht: Landesgericht für ZRS Graz
> UID: ATU54713500
>
>
>
>
> --
> Stefano
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230607/c82d85eb/attachment-0001.html>
More information about the petsc-users
mailing list