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
