1
$\begingroup$

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?

Sorry about the length of this question. Thanks for your time!

$\endgroup$

closed as off-topic by Wolfgang Bangerth, Christian Clason, nicoguaro, Brian Borchers, Mauro Vanzetto Oct 29 '17 at 16:18

This question appears to be off-topic. The users who voted to close gave this specific reason:

If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Please run with -log_view and attach the result. $\endgroup$ – Jed Brown Oct 26 '17 at 18:22
  • $\begingroup$ @JedBrown: done! $\endgroup$ – EM_IE Oct 26 '17 at 19:32
  • 2
    $\begingroup$ This is a question you should ask on the PETSc specific forums. The people who can help you with your question are all there, not necessarily here. $\endgroup$ – Wolfgang Bangerth Oct 27 '17 at 0:25
  • 4
    $\begingroup$ I agree with @WolfgangBangerth that this question is more appropriate on petsc-users. Your circumstances just don't generalize that well. From the log above, you can see that your program takes a total of 405 seconds, with about 157 seconds of that spent in PETSc, almost exactly evenly split between solve time (KSPSolve consists of MatMult, MatSolve, and vector operations) and setup (PCSetUp consists of factorization). ILU(3) is quite expensive for this matrix, and the convergence isn't particularly fast. You might benefit a greatly from a preconditioner well suited to your application. $\endgroup$ – Jed Brown Oct 27 '17 at 1:05