Petsc using VecDuplicate in solution process
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. 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
Il giorno mar 6 giu 2023 alle ore 09:24 Pichler, Franz < [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
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
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]> Sent: Tuesday, June 6, 2023 4:40 PM To: Pichler, Franz <[email protected]> Cc: [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
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
participants (3)
-
Barry Smith -
Pichler, Franz -
Stefano Zampini