[petsc-users] Solving tridiagonal hermitian generalized eigenvalue problem with SLEPC

Mon Aug 18 15:24:51 CDT 2014

Toon :
> Sorry, I forgot to mention: I am looking for the Eigenvalues that are the
> largest, not in absolute value, but along the real axis.

Then, you do not need shift-invert, therefore, should not use LU
matrix factorization.

Hong
>
>
> On 18 August 2014 21:54, Toon Weyens <tweyens at fis.uc3m.es> wrote:
>>
>> Hi, thanks for the answers!
>>
>> I think I expressed myself wrong: I can indeed get it to work with just
>> using AIJ matrices, as in example 13. This is the way that I am currently
>> solving my problem. There are only two issues:
>> 1. memory is indeed important so I would certainly like to decrease it by
>> one third if possible :-) The goal is to make the simulations as fast and
>> light as possible to be able to perform parameter studies (on the stability
>> of MHD configurations).
>>
>> 2. I have played around a little bit with the different solvers but it
>> appears that the standard method and the Arnoldi with explicit restart
>> method are the best. Some of the others don't converge and some are slower.
>>
>> The thing is that in the end the matrices that I use are large but they
>> have a very easy structure: hermitian tri-diagonal. That's why, I think,
>> slepc usually converges in a few iterations (correct me if I'm wrong).
>>
>> The problem is that sometimes, when I consider more grid points, the
>> solver doesn't work any more because apparently it uses the LU decomposition
>> (not sure for the matrix A or B in A x = lambda B x) and there is a zero
>> pivot element (see below for error message). In other words: the matrices
>> become almost singular. This is a characteristic of the numerical method I
>> think. Is there any way to fix this by setting a more precise threshold or
>> something?
>>
>> Also, there is one more thing: is it possible and useful to use the
>> non-zero structure of A when defining matrix B? Does this affect perfomance?
>>
>> Thanks in advance!
>> Toon
>>
>>
>> Error messages: identical for KRYLOV-SCHUR and ARNOLDI
>>
>> [0]PETSC ERROR: --------------------- Error Message
>> ------------------------------------
>> [0]PETSC ERROR: Detected zero pivot in LU factorization:
>> see http://www.mcs.anl.gov/petsc/documentation/faq.html#ZeroPivot!
>> [0]PETSC ERROR: Zero pivot row 19 value 6.57877e-16 tolerance 2.22045e-14!
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: Petsc Release Version 3.4.5, Jun, 29, 2014
>> [0]PETSC ERROR: See docs/changes/index.html for recent updates.
>> [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.
>> [0]PETSC ERROR: See docs/index.html for manual pages.
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: ./PB3D on a debug-complex named toon-XPS-L501X by toon Mon
>> Aug 18 21:47:41 2014
>> [0]PETSC ERROR: Libraries linked from
>> /opt/petsc/petsc-3.4.5/debug-complex/lib
>> [0]PETSC ERROR: Configure run at Mon Aug 11 10:08:29 2014
>> [0]PETSC ERROR: Configure options PETSC_ARCH=debug-complex
>> --with-scalar-type=complex --with-debugging
>> [0]PETSC ERROR:
>> ------------------------------------------------------------------------
>> [0]PETSC ERROR: MatPivotCheck_none() line 589 in
>> /opt/petsc/petsc-3.4.5/include/petsc-private/matimpl.h
>> [0]PETSC ERROR: MatPivotCheck() line 608 in
>> /opt/petsc/petsc-3.4.5/include/petsc-private/matimpl.h
>> [0]PETSC ERROR: MatLUFactorNumeric_SeqAIJ_Inode() line 1837 in
>> /opt/petsc/petsc-3.4.5/src/mat/impls/aij/seq/inode.c
>> [0]PETSC ERROR: MatLUFactorNumeric() line 2889 in
>> /opt/petsc/petsc-3.4.5/src/mat/interface/matrix.c
>> [0]PETSC ERROR: PCSetUp_LU() line 152 in
>> /opt/petsc/petsc-3.4.5/src/ksp/pc/impls/factor/lu/lu.c
>> [0]PETSC ERROR: PCSetUp() line 890 in
>> /opt/petsc/petsc-3.4.5/src/ksp/pc/interface/precon.c
>> [0]PETSC ERROR: KSPSetUp() line 278 in
>> /opt/petsc/petsc-3.4.5/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: PCSetUp_Redundant() line 170 in
>> /opt/petsc/petsc-3.4.5/src/ksp/pc/impls/redundant/redundant.c
>> [0]PETSC ERROR: PCSetUp() line 890 in
>> /opt/petsc/petsc-3.4.5/src/ksp/pc/interface/precon.c
>> [0]PETSC ERROR: KSPSetUp() line 278 in
>> /opt/petsc/petsc-3.4.5/src/ksp/ksp/interface/itfunc.c
>> [0]PETSC ERROR: STSetUp_Shift() line 126 in
>> /opt/slepc/slepc-3.4.4/src/st/impls/shift/shift.c
>> [0]PETSC ERROR: STSetUp() line 294 in
>> /opt/slepc/slepc-3.4.4/src/st/interface/stsolve.c
>> [0]PETSC ERROR: EPSSetUp() line 215 in
>> /opt/slepc/slepc-3.4.4/src/eps/interface/setup.c
>> [0]PETSC ERROR: EPSSolve() line 90 in
>> /opt/slepc/slepc-3.4.4/src/eps/interface/solve.c
>>
>>
>>
>> On 18 August 2014 20:52, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>>
>>>
>>>    I would recommend just using AIJ matrices and not worry about SBAIJ at
>>> this stage. SBAIJ can save you memory but not operations and unless memory
>>> is really tight is not likely a useful optimization.
>>>
>>>    What eigensolvers have you tried and which eigenvalues do you want?
>>>
>>>    Barry
>>>
>>> On Aug 18, 2014, at 5:46 AM, Toon Weyens <tweyens at fis.uc3m.es> wrote:
>>>
>>> > Dear all,
>>> >
>>> > I am using PETSC and SLEPC to simulate a problem in MHD, described in
>>> > http://scitation.aip.org/content/aip/journal/pop/21/4/10.1063/1.4871859.
>>> >
>>> > I have a bit of experience with MPI but not too much with PETSC and
>>> > SLEPC. So after reading both user manuals and also the relevant chapters of
>>> > the PETSC developers manual, I still can't get it to work.
>>> >
>>> > The problem that I have to solve is a large generalized eigenvalue
>>> > system where the matrices are both Hermitian by blocks and tridiagonal,
>>> > e.g.:
>>> >
>>> > ( A11  A12    0     0       0   )                      ( B11  B12    0
>>> > 0       0   )
>>> > ( A12* A22  A23   0      0   )                       ( B12* B22  B23
>>> > 0      0   )
>>> > (   0    A23* A33  A34   0   )   =   lambda    (   0    B23* B33  B34
>>> > 0   )
>>> > (   0      0    A34* A44  A45)                       (   0      0
>>> > B34* B44  B45)
>>> > (   0      0    A45* A55    0  )                       (   0      0
>>> > B45* B55    0  )
>>> >
>>> > where Aii = Aii*, with * the Hermitian conjugate. I apologize for the
>>> > ugly representation.
>>> >
>>> > The dimensions of both A and B are around 50 to 100 blocks (as there is
>>> > a block per discretized point) and the blocks themselves can vary from 1 to
>>> > more than 100x100 as well (as they correspond to a spectral decomposition).
>>> >
>>> > Now, my question is: how to solve this economically?
>>> >
>>> > What I have been trying to do is to make use of the fact that the
>>> > matrices are Hermitian and by using matcreatesbaij and through the
>>> > recommended matcreate, matsettype(matsbaij), etc.
>>> >
>>> > Could someone help me out? All help would be greatly appreciated!
>>> >
>>> > Thank you in advance,
>>> > Toon
>>> > UC3M
>>>
>>
>