[petsc-users] KSPSolve takes more in each iteration

[Barry Smith]

  I don't see any for () loop ?

  solveAxb() seems to create KSP objects etc but where are they destroyed? If you call solveAxb repeatedly you are going to be creating more and more PETSc objects and never destroying them

  Barry

    
> On Jul 20, 2015, at 11:29 PM, Orxan Shibliyev <orxan.shibli at gmail.com> wrote:
> 
> The structure of matrix remains the same. Before for loop I use MatMPIAIJSetPreallocation to allocate for once. But I set d_nz and o_nz more than enough since I don't know their exact values but minimums. Number of iterations does not change too much and it is 32-33. This problem happens even for sequential code. I also don't explicitly allocate memory in the loop. All allocations are done in the constructor of user defined struct Petsc:
> 
> struct Solver
> {
>     struct Petsc
>     {
>         Vec x, b;
>         Mat A;
>         KSP ksp;
>         PC pc;
>         MPI_Comm world;
>         PetscInt n;
>         PetscInt bs;
>         PetscInt xGlobalSize;
>         PetscInt xLocalSize;
>         PetscInt vecFirst, vecLast;
>         PetscInt first, last;
>         double* DX;
>         
>         Petsc (Grid& gr);
>         void solveAxb (Grid& gr);
>         void finalize();
>     } petsc;
> 
>     // ...
> }
> 
> Solver::Petsc::Petsc (Grid& gr)
> {
>     world = PETSC_COMM_WORLD;
>     n = gr.n_in_elm;
>     bs = N_VAR;
>     xGlobalSize = n*bs;
>     DX = (double*) malloc (xGlobalSize*sizeof(double));
>     
>     // set x
>     VecCreate (world, &x);
>     VecSetType (x, VECSTANDARD);
>     VecSetSizes (x, PETSC_DECIDE, xGlobalSize);
>     VecGetLocalSize (x, &xLocalSize);
>     VecGetOwnershipRange (x, &vecFirst, &vecLast);
> 
>     // set b
>     VecDuplicate (x, &b);
>     
>     // set A
>     MatCreate (world, &A);        
>     MatSetType (A, MATMPIAIJ);
>     MatSetSizes (A, PETSC_DECIDE, PETSC_DECIDE, xGlobalSize, xGlobalSize);
>     MatMPIAIJSetPreallocation (A, 4*bs, NULL, 4*bs, NULL);
>     MatGetOwnershipRange (A, &first, &last);
> }
> 
> void Solver::Petsc::solveAxb (Grid& gr)
> {
>     // set values of b
>     // set values of A
> 
>     KSPCreate (world, &ksp);
>     KSPGetPC (ksp, &pc);
>     KSPSetOperators (ksp, A, A);
>     PCSetType (pc, PCSOR);    
>     KSPSetType (ksp, KSPGMRES); 
>     KSPSetFromOptions (ksp);    
>     KSPSolve (ksp, b, x);
>     
>     double *dx = NULL;
>     VecGetArray (x, &dx);
>     MPI_Allgatherv (...); // dx -> DX
>     VecRestoreArray (x, &dx);
> }
> 
> On Mon, Jul 20, 2015 at 8:20 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
>   Is the matrix changing in the for loop?
> 
>   Are the number of iterations increasing for each KSP? run with -ksp_monitor
> 
> 
> 
> > On Jul 20, 2015, at 8:01 PM, Matthew Knepley <knepley at gmail.com> wrote:
> >
> > It sounds like you are allocating memory each time. You can check using -malloc_dump
> > after a few iterates.
> >
> >    Matt
> >
> > On Mon, Jul 20, 2015 at 7:35 PM, Orxan Shibliyev <orxan.shibli at gmail.com> wrote:
> > I am solving a non-linear problem. Basically I do the following:
> >
> > ... Allocate A, x, b, ksp ...
> >
> > for_loop
> >     ... call myFunction ...
> >
> > where,
> >
> > myFunction:
> >     ...
> >     startTime
> >     KSPSolve
> >     stopTime
> >     ....
> >
> > I record time for KSPSolve. It takes more each time until it becomes impossible to wait. I put the code on a cluster with (14 x 1) processors. Any idea?
> >
> >
> >
> > --
> > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> > -- Norbert Wiener
> 
>