# Why does PETSc take unexpectedly long to set up its KSP solver with a custom preconditioner? [closed]

I am attempting to solve a large system, $\bf{Ax} = \bf{b}$ with the help of PETSc. Due to the size of the problem, I'm using a matrix-free approach, where $\bf{A}$ is just a shell. I'm also providing my own preconditioner (which is not a shell), and I'm using ilu(2) factorization on the preconditioner.

The issue: the setup stage of the solver (see relevant block of code below) is taking an extremely long time. I suspect that this would mostly be the ilu of the preconditioner that's taking a while. I know that it is expected that the ilu will take time, but I'm concerned for the following reason: when I try to solve the same problem using a direct solver (using MKL Lapacke to find the LU decomposition and then inversion of $\bf{A}$, outside of PETSc), the LU is 10x faster. I expect that PETSc's ilu should take a comparable time to the full LU factorization, but it is 10x slower, which seems very strange. (By the way, you may ask, why even use an iterative solver if I can do it with LU, but this example is not as large as the ones I actually want to run, at which point I won't be able to use a direct solver).

Here is the code snippet relevant to the issue:

MatCreateShell(comm, Nu, Nu, Nu, Nu, ctx, &A_shell);
MatShellSetOperation(A_shell, MATOP_MULT, (void(*)(void))usermult);

KSPCreate(comm, &solver);
KSPSetOperators(solver, A_shell, PreconditionerMatrix);
KSPSetInitialGuessNonzero(solver, PETSC_TRUE);
KSPSetNormType(solver, KSP_NORM_UNPRECONDITIONED);

KSPSetFromOptions(solver);
KSPSetUp(solver);


Things I know / have tried:

• The condition number of the matrix can be as large as $10^7$, but I don't think that my issue has anything to do with that, because again, the time sink is in the setup. If this was the problem, it would also manifest itself when I do the full LU decomposition, but it doesn't.
• I know that I have to provide a large enough fill factor guess for the post-ilu matrix, and I've set the option -pc_factor_fill to 3. After running the code with -info, I confirmed that this is sufficient to prevent any reallocation of memory. Interesting side note: when I do run it with -info, it reports the required fill factor fairly quickly. Does this mean that it actually does the ilu expectedly fast, but then gets stuck somewhere else? Am I barking up the wrong tree? Here's what it reports:

  [0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
[0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
[0] PCSetUp(): Setting up PC for first time
[0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
[0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
[0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
[0] PetscCommDuplicate(): Using internal PETSc communicator 7412512 20851120
[0] MatILUFactorSymbolic_SeqAIJ(): Reallocs 0 Fill ratio:given 3. needed 1.7385
[0] MatILUFactorSymbolic_SeqAIJ(): Run with -[sub_]pc_factor_fill 1.7385 or use
[0] MatILUFactorSymbolic_SeqAIJ(): PCFactorSetFill([sub]pc,1.7385);
[0] MatILUFactorSymbolic_SeqAIJ(): for best performance.
[0] MatSeqAIJCheckInode_FactorLU(): Found 2030 nodes of 6096. Limit used: 5. Using Inode routines


Then it gets stuck for a really long time...so maybe I set the fill factor too large? I tried the same thing again with a fill factor of 2 instead, but it made no difference.

• I have already made sure I'm not using the debugging installation of PETSc; when timing the code, I'm definitely using --with-debugging=0 in the PETSc configuration step.

• I am not using any parallelization.

Here is the output generated with -log_view:

Using Petsc Release Version 3.7.6, Apr, 24, 2017

Max       Max/Min        Avg      Total
Time (sec):           4.057e+02      1.00000   4.057e+02
Objects:              7.050e+02      1.00000   7.050e+02
Flops:                2.161e+11      1.00000   2.161e+11  2.161e+11
Flops/sec:            5.327e+08      1.00000   5.327e+08  5.327e+08
MPI Messages:         0.000e+00      0.00000   0.000e+00  0.000e+00
MPI Message Lengths:  0.000e+00      0.00000   0.000e+00  0.000e+00
MPI Reductions:       0.000e+00      0.00000

Flop counting convention: 1 flop = 1 real number operation of type (multiply/divide/add/subtract)
e.g., VecAXPY() for real vectors of length N --> 2N flops
and VecAXPY() for complex vectors of length N --> 8N flops

Summary of Stages:   ----- Time ------  ----- Flops -----  --- Messages ---  -- Message Lengths --  -- Reductions --
Avg     %Total     Avg     %Total   counts   %Total     Avg         %Total   counts   %Total
0:      Main Stage: 4.0567e+02 100.0%  2.1612e+11 100.0%  0.000e+00   0.0%  0.000e+00        0.0%  0.000e+00   0.0%

------------------------------------------------------------------------------------------------------------------------
See the 'Profiling' chapter of the users' manual for details on interpreting output.
Phase summary info:
Count: number of times phase was executed
Time and Flops: Max - maximum over all processors
Ratio - ratio of maximum to minimum over all processors
Mess: number of messages sent
Avg. len: average message length (bytes)
Reduct: number of global reductions
Global: entire computation
Stage: stages of a computation. Set stages with PetscLogStagePush() and PetscLogStagePop().
%T - percent time in this phase         %F - percent flops in this phase
%M - percent messages in this phase     %L - percent message lengths in this phase
%R - percent reductions in this phase
Total Mflop/s: 10e-6 * (sum of flops over all processors)/(max time over all processors)
------------------------------------------------------------------------------------------------------------------------
Event                Count      Time (sec)     Flops                             --- Global ---  --- Stage ---   Total
Max Ratio  Max     Ratio   Max  Ratio  Mess   Avg len Reduct  %T %F %M %L %R  %T %F %M %L %R Mflop/s
------------------------------------------------------------------------------------------------------------------------

--- Event Stage 0: Main Stage

MatMult              624 1.0 3.2336e+01 1.0 4.23e+10 1.0 0.0e+00 0.0e+00 0.0e+00  8 20  0  0  0   8 20  0  0  0  1307
MatMultAdd          2480 1.0 5.4122e-02 1.0 5.00e+07 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0   923
MatMultTranspose    3100 1.0 2.8153e+00 1.0 9.05e+09 1.0 0.0e+00 0.0e+00 0.0e+00  1  4  0  0  0   1  4  0  0  0  3215
MatSolve             608 1.0 4.6020e+01 1.0 7.41e+10 1.0 0.0e+00 0.0e+00 0.0e+00 11 34  0  0  0  11 34  0  0  0  1611
MatLUFactorNum         1 1.0 4.4004e+01 1.0 9.31e+10 1.0 0.0e+00 0.0e+00 0.0e+00 11 43  0  0  0  11 43  0  0  0  2115
MatILUFactorSym        1 1.0 3.4659e+01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  9  0  0  0  0   9  0  0  0  0     0
MatAssemblyBegin      27 1.0 1.5497e-05 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatAssemblyEnd        27 1.0 6.8570e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatGetRow        10650098 1.0 6.6575e+00 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  2  0  0  0  0   2  0  0  0  0     0
MatGetRowIJ            1 1.0 9.5367e-07 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatGetOrdering         1 1.0 1.3018e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatZeroEntries        39 1.0 1.5883e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatMatMult             2 1.0 3.5977e+00 1.0 5.97e+09 1.0 0.0e+00 0.0e+00 0.0e+00  1  3  0  0  0   1  3  0  0  0  1660
MatMatMultSym          2 1.0 7.2353e-01 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
MatMatMultNum          2 1.0 2.8741e+00 1.0 5.97e+09 1.0 0.0e+00 0.0e+00 0.0e+00  1  3  0  0  0   1  3  0  0  0  2078
VecMDot              302 1.0 4.2068e-02 1.0 2.06e+08 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  4907
VecNorm              621 1.0 1.4746e-02 1.0 3.03e+07 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  2054
VecScale             314 1.0 2.0843e-03 1.0 7.66e+06 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  3674
VecCopy              636 1.0 1.0492e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecSet              4096 1.0 4.1316e-02 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecAXPY             5278 1.0 9.9347e-02 1.0 1.96e+08 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  1977
VecAYPX              306 1.0 4.3933e-03 1.0 7.46e+06 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  1698
VecMAXPY             608 1.0 5.9476e-02 1.0 4.16e+08 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  6993
VecAssemblyBegin    2493 1.0 1.5733e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecAssemblyEnd      2493 1.0 1.5340e-03 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
VecNormalize         314 1.0 1.0375e-02 1.0 2.30e+07 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  2214
KSPGMRESOrthog       302 1.0 7.3104e-02 1.0 4.13e+08 1.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0  5647
KSPSetUp               1 1.0 1.2398e-04 1.0 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00  0  0  0  0  0   0  0  0  0  0     0
KSPSolve               4 1.0 7.8504e+01 1.0 1.17e+11 1.0 0.0e+00 0.0e+00 0.0e+00 19 54  0  0  0  19 54  0  0  0  1491
PCSetUp                1 1.0 7.8663e+01 1.0 9.31e+10 1.0 0.0e+00 0.0e+00 0.0e+00 19 43  0  0  0  19 43  0  0  0  1183
PCApply              608 1.0 4.6022e+01 1.0 7.41e+10 1.0 0.0e+00 0.0e+00 0.0e+00 11 34  0  0  0  11 34  0  0  0  1611
------------------------------------------------------------------------------------------------------------------------

Memory usage is given in bytes:

Object Type          Creations   Destructions     Memory  Descendants' Mem.
Reports information only for process 0.

--- Event Stage 0: Main Stage

Matrix    33             30   1306648980     0.
Vector   663            663     65247472     0.
Krylov Solver     1              1        35264     0.
Preconditioner     1              1         1008     0.
Viewer     2              0            0     0.
Index Set     5              5        87624     0.
========================================================================================================================
Average time to get PetscTime(): 5.96046e-07
#PETSc Option Table entries:
-ksp_atol 1e-8
-ksp_converged_reason
-ksp_monitor
-ksp_monitor_true_residual
-ksp_rtol 1e-8
-log_view
-pc_factor_fill 3
-pc_factor_levels 2
-pc_type ilu
#End of PETSc Option Table entries
Compiled with FORTRAN kernels
Compiled with full precision matrices (default)
sizeof(short) 2 sizeof(int) 4 sizeof(long) 8 sizeof(void*) 8 sizeof(PetscScalar) 16 sizeof(PetscInt) 4
Configure options: PETSC_ARCH=arch-linux2-cxx-nodebug --with-scalar-type=complex --with-fortran-kernels=1 --with-clanguage=c++ --with-debugging=0 --with-cxx=g++ CXXOPTFLAGS=-O3 COPTFLAGS=O3 FOPTFLAGS=-O3 --download-openmpi --with-blaslapack-dir=/opt/intel/mkl


Questions

• I know that there may not be an easy / obvious solution to this, but at least I'd like insight on why this happens, and what PETSc is doing internally that takes this long
• If there is no clear solution, what are some steps I can take to try and mitigate this or investigate this further?
• Is this expected / normal, and should I stop worrying about it and just suck it up?

• Please run with -log_view and attach the result. – Jed Brown Oct 26 '17 at 18:22