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_fillto 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 routinesThen 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=0in 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!
-log_viewand attach the result. $\endgroup$