[petsc-users] SLEPc spectrum slicing

[Jeff Steward]

Thank you very much for your quick response Jose. I really appreciate it.
Your point about load balancing even with the same number of pairs in the
interval is well taken.

For now I can provide hints with EPSKrylovSchurSetSubintervals and that
should work fine. However, my problem is reaching the upper limits of
memory with plans to increase the size, so I'd like to work on an approach
to spectrum slicing with iterative/shell methods for the long run.

The idea for spectrum slicing as I discussed in my previous message is
based on the paper of Aktulga et al (2014):

https://math.berkeley.edu/~linlin/publications/ParallelEig.pdf

They address the issues involved with computing the (very expensive) matrix
inertia. These experiments use PARPACK and FEAST to do a parallel spectrum
slicing method. Of course PARPACK and FEAST have their own issues, and
SLEPc is much easier to use and more flexible, extensible, etc. I think the
MSIL approach they outline might be beneficial for my relatively large
eigenproblem where I need many eigenpairs.

Computing the matrix inertia for such an approach is likely to be the
bottleneck at least in terms of memory. How far do you think the incomplete
Cholesky with MUMPS would scale?

I will send a separate e-mail to slepc-maint regarding my problem.

Thank you again!

Best wishes,

Jeff

On Thu, May 19, 2016 at 2:19 PM, Jose E. Roman <jroman at dsic.upv.es> wrote:

>
> > El 19 may 2016, a las 22:28, Jeff Steward <jeffsteward at gmail.com>
> escribió:
> >
> > I have some questions regarding spectrum slicing in SLEPc, especially in
> the new version 3.7, as I see a line in the Change Notes that I don't quite
> understand ("in spectrum slicing in multi-communicator mode now it is
> possible to update the problem matrices directly on the sub-communicators".)
>
> This is intended for applications that require solving a sequence of
> eigenproblems, where the problem matrices in a given problem are a
> modification of the previous ones. Instead of updating the matrices in the
> parent communicator it is now possible to do the update in the partitioned
> communicators. This avoids lots of communication required for moving data
> to/from the parent communicator. This is certainly an "advanced feature".
>
> >
> > 1) In light of this statement, how should I divide my Mat to best work
> with eps_krylovschur_partitions? Let's say I have 384 processors and
> eps_krylovschur_partitions=48, so there will be 8 processors in each group.
> Should I distribute the matrix over all 384 processors (so let's say this
> gives 32 rows per processor) and have the entire communicator call EPS
> solve, or should I (can I?) distribute the matrix over each of the 8 groups
> (giving say 1536 rows per processor) and have the 48 subcommunicators call
> EPSSolve? Am I correct in thinking that distributing over all 384
> processors requires collecting and duplicating the matrix to the 48
> different groups?
>
> The idea is that you create the matrix in the initial communicator
> (typically PETSC_COMM_WORLD) and then EPS will automatically create the
> partitioned sub-communicators and replicate the matrices via
> MatCreateRedundantMatrix(). You have to choose the number of partitions and
> the number of processes in the original communicator that are best suited
> for you application. If your matrix is relatively small but you need to
> compute many eigenvalues, then you can consider setting as many partitions
> as processes in the original communicator (hence the size of each partition
> is 1). But it will be necessary to communicate a lot for setting up data in
> the sub-communicators. The main point of sub-communicators is to workaround
> the limited scalability of parallel direct linear solvers (e.g. MUMPS) when
> one wants to use may processes.
>
> > 2) If I'm understanding it correctly, the approach for Krylov-Schur
> spectrum slicing described in the manual for SLEPc seems like a wasteful
> default method. From what I gather, regions are divided by equal distance,
> and different regions are bound to end up with different (and potentially
> vastly different) numbers of eigenpairs. I understand the user can provide
> their own region intervals, but wouldn't it be better for SLEPc to first
> compute the matrix inertias at some given first guess regions, interpolate,
> then fix the spectra endpoints so they contain approximately the same
> number in each region? An option for logarithmically spaced points rather
> than linearly spaced points would be helpful as well, as for the problem I
> am looking at the spectrum decays in this way (few large eigenvalues with
> an exponential decay down to many smaller eigenvalues). I require
> eigenpairs with eigenvalues that vary by several orders of magnitude (1e-3
> to 1e3), so the linear equidistant strategy is hopeless.
>
> If you have an a priori knowledge of eigenvalue distribution then you can
> use EPSKrylovSchurSetSubintervals() to give hints. We do not compute
> inertia to get a rough estimation of the distribution because in the
> general case computing inertia is very costly (it requires a factorization
> in our implementation). The main intention of
> EPSKrylovSchurSetSubintervals() is also the case where a sequence of
> eigenproblems must be solved, so the solution of one problem provides an
> estimation of eigenvalue distribution for the next problem.
>
> Anyway, giving subintervals that roughly contain the same number of
> eigenvalues does not guarantee that the workload will be balanced, since
> convergence may be quite different from one subinterval to the other.
>
>
> > 3) The example given for spectrum slicing in the user manual is
> >
> >  mpiexec -n 20 ./ex25 -eps_interval 0.4,0.8 -eps_krylovschur_partitions 4
> >                   -st_type sinvert -st_ksp_type preonly -st_pc_type
> cholesky
> >                   -st_pc_factor_mat_solver_package mumps
> -mat_mumps_icntl_13 1
> >
> > which requires a direct solver. If I can compute the matrix inertias
> myself and come up with spectral regions and subcommunicators as described
> above, is there a way to efficiently use SLEPc with an iterative solver?
> How about with a matrix shell? (I'm getting greedy now ^_^).
>
> I don't know how you can possibly compute the "exact" inertia cheaply. In
> that case yes, it would be worth using a shell matrix but the solver is
> probably not prepared for this case. If you want us to have a look at this
> possibility, send more details to slepc-maint.
>
> Jose
>
> >
> > I would really appreciate any help on these questions. Thank you for
> your time.
> >
> > Best wishes,
> >
> > Jeff
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20160519/4dfb31a4/attachment-0001.html>