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 <[email protected]> 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 <[email protected] > <mailto:[email protected]>> > Sent: Tuesday, June 6, 2023 4:40 PM > To: Pichler, Franz <[email protected] <mailto:[email protected]>> > Cc: [email protected] <mailto:[email protected]> > Subject: Re: [petsc-users] Petsc using VecDuplicate in solution process > > > > Il giorno mar 6 giu 2023 alle ore 09:24 Pichler, Franz <[email protected] > <mailto:[email protected]>> 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 > [email protected] <mailto:[email protected]> > 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
