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 <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 <[email protected]>
> > 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
>
>
