Re: [petsc-users] Petsc using VecDuplicate in solution process

[Pichler, Franz]

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 
#include 
#include 
#include 
#include 
#include 
#include  /*I "petscksp.h" I*/
#include 
#include 
#include 
#include 

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, _ksp);
KSPSetFromOptions(petsc_ksp);
KSPGetPC(petsc_ksp, _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& dsp,std::vector& 
col,std::vector& 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:"<& val){
double *ptr_ele = val.data();
VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, n_rows, NULL,_b);
VecPlaceArray(petsc_b, ptr_ele);
}
void set_sol(std::vector& val){
double *ptr_ele = val.data();
VecCreateSeqWithArray(PETSC_COMM_WORLD, 1, n_rows, NULL,_x);
VecPlaceArray(petsc_x, ptr_ele);
}

int solve(bool force_rebuild) {
int solver_stat = 0;
KSPGetPC(petsc_ksp, _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, );
if (reason < 0){
KSPGetIterationNumber(petsc_ksp, _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, _iter);
return linear_iter;
}
};
void change_rhs(int i, int n_rows,std::vector){
for(int 
row=0;row& vals){
int nnz= n_rows*3-2;
for(int row=0;row& dsp,std::vector& 
col,std::vector& val ){
int nnz= n_rows*3-2;
std::cout<<"SETCRS ROWS:"<& dsp,std::vector& 
col,std::vector& val ,std::vector& sol,std::vector rhs){
int n_rows=dsp.size()-1;
double res=0;
for (int row=0;row rhs(n_rows);
std::vector sol(n_rows);
std::vector dsp;
std::vector cols;
std::vector 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:"<
Sent: Tuesday, June 6, 2023 4:40 PM
To: Pichler, Franz 
Cc: petsc-users@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 
mailto:franz.pich...@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.pich..

[petsc-users] Petsc using VecDuplicate in solution process

[Pichler, Franz]

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
franz.pich...@v2c2.at
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

Re: [petsc-users] Petsc ObjectStateIncrease without proivate header

[Pichler, Franz]

Hello, very sorry for the very late reply, but thank you even more for the ver 
helpful suggestion!
Using valgrinds callgrind I cann see that matassemblyBegin/End takes some 
cycles, but I guess I can take take,
a perfect solution would have no overhead,

Ii still changed the code to get rid of the private dependency,

Thank you!

From: Stefano Zampini 
Sent: Monday, May 8, 2023 3:31 PM
To: Pichler, Franz 
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] Petsc ObjectStateIncrease without proivate header

You can achieve the same effect by calling MatAssemblyBegin/End

Il giorno lun 8 mag 2023 alle ore 15:54 Pichler, Franz 
mailto:franz.pich...@v2c2.at>> ha scritto:
Hello,
i am using petsc in a single cpu setup where I have preassembled crs matrices 
that I wrap via PetSC’s
MatCreateSeqAIJWithArrays  Functionality.

Now I manipulate the values of these matrices (wohtout changing the sparsity) 
without using petsc,

When I now want to solve again I have to call
PetscObjectStateIncrease((PetscObject)petsc_A);

So that oetsc actually solves again (otherwise thinking nothing hs changed ,
This means I have to include the private header
#include 

Which makes a seamingless implementation of petsc into a cmake process more 
complicate (This guy has to be stated explicitly in the cmake process at the 
moment)

I would like to resolve that by “going” around the private header,
My first intuition was to increase the state by hand
((PetscObject)petsc_A_aux[the_sys])->state++;
This is the definition of petscstateincrease in the header. This throws me an
error: invalid use of incomplete type ‘struct _p_PetscObject’

compilation error.

Is there any elegeant way around this?
This is the first time I use the petsc mailing list so apologies for any 
beginners mistake I did in formatting or anything else.

Best regards

Franz Pichler



--
Stefano

[petsc-users] Petsc ObjectStateIncrease without proivate header

[Pichler, Franz]

Hello,
i am using petsc in a single cpu setup where I have preassembled crs matrices 
that I wrap via PetSC's
MatCreateSeqAIJWithArrays  Functionality.

Now I manipulate the values of these matrices (wohtout changing the sparsity) 
without using petsc,

When I now want to solve again I have to call
PetscObjectStateIncrease((PetscObject)petsc_A);

So that oetsc actually solves again (otherwise thinking nothing hs changed ,
This means I have to include the private header
#include 

Which makes a seamingless implementation of petsc into a cmake process more 
complicate (This guy has to be stated explicitly in the cmake process at the 
moment)

I would like to resolve that by "going" around the private header,
My first intuition was to increase the state by hand
((PetscObject)petsc_A_aux[the_sys])->state++;

This is the definition of petscstateincrease in the header. This throws me an
error: invalid use of incomplete type 'struct _p_PetscObject'

compilation error.

Is there any elegeant way around this?
This is the first time I use the petsc mailing list so apologies for any 
beginners mistake I did in formatting or anything else.

Best regards

Franz Pichler