SLEPcWrappers::SolverKrylovSchur Zero pivot problem

[Léonhard YU]

Hello Everyone:

I am dealing with a eigenvalue problem with SLEPc, but the program reports this error:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Zero pivot in LU factorization: https://petsc.org/release/faq/#zeropivot
[0]PETSC ERROR: Zero pivot row 1 value 0. tolerance 2.22045e-14

I have reviewed the documentation. It says that (Frequently Asked Questions (FAQ) — PETSc 3.16.4 documentation)

A zero pivot in LU, ILU, Cholesky, or ICC sparse factorization does not always mean that the matrix is singular. You can use

-pc_factor_shift_type nonzero -pc_factor_shift_amount [amount]

or

-pc_factor_shift_type positive_definite -[level]_pc_factor_shift_type nonzero -pc_factor_shift_amount [amount]

or

-[level]_pc_factor_shift_type positive_definite

to prevent the zero pivot. [level] is “sub” when lu, ilu, cholesky, or icc are employed in each individual block of the bjacobi or ASM preconditioner. [level] is “mg_levels” or “mg_coarse” when lu, ilu, cholesky, or icc are used inside multigrid smoothers or to the coarse grid solver. See PCFactorSetShiftType(), PCFactorSetShiftAmount().

This error can also happen if your matrix is singular, see MatSetNullSpace() for how to handle this. If this error occurs in the zeroth row of the matrix, it is likely you have an error in the code that generates the matrix.

So I write the code as follows:

   PETScWrappers::PreconditionBlockJacobi::AdditionalData data;
    PETScWrappers::PreconditionBlockJacobi  preconditioner(mpi_communicator, data);
    SolverControl linear_solver_control (dof_handler.n_dofs(),
                                         1e-9, false, false);
    PETScWrappers::SolverCG linear_solver(linear_solver_control,
                                          mpi_communicator);

    PCFactorSetShiftType(preconditioner.get_pc(), MAT_SHIFT_NONZERO);
//    PCFactorSetShiftAmount(preconditioner.get_pc(), 3e-14);

    linear_solver.initialize(preconditioner);

    SolverControl solver_control (100, 1e-9,false,false);
    SLEPcWrappers::SolverKrylovSchur eigensolver(solver_control,
                                                 mpi_communicator);

    eigensolver.set_which_eigenpairs(EPS_SMALLEST_REAL);

    eigensolver.set_problem_type(EPS_GHEP);

    SLEPcWrappers::TransformationShift spectral_transformation(mpi_communicator);
    spectral_transformation.set_solver(linear_solver);
    eigensolver.set_transformation(spectral_transformation);

    eigensolver.solve(stiffness_matrix,
                      mass_matrix,
                      eigenvalues,
                      eigenfunctions,
                      eigenfunctions.size());

But the error occurs as the same. I don't know how to solve it.

Thanks for your help.

Best Wishes.

[Wolfgang Bangerth]

Leonhard:
have you checked that the matrix really is invertible? We use the wrappers as
is for step-36, if I recall correctly. A good starting point for asking why
your approach doesn't work would be to ask what is different between your
program and step-36.

Best
W.


On 2/25/22 07:10, Léonhard YU wrote:
> *** Caution: EXTERNAL Sender ***

>
> Hello Everyone:
>
> I am dealing with a eigenvalue problem with SLEPc, but the program reports
> this error:
>
> [0]PETSC ERROR: --------------------- Error Message
> --------------------------------------------------------------
> [0]PETSC ERROR: Zero pivot in LU factorization:
> https://petsc.org/release/faq/#zeropivot
> [0]PETSC ERROR: Zero pivot row 1 value 0. tolerance 2.22045e-14
>
> I have reviewed the documentation. It says that (Frequently Asked Questions
> (FAQ) — PETSc 3.16.4 documentation

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Frelease%2Ffaq%2F%23zeropivot&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cff03e4917646439d68e908d9f868a339%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637813950590384693%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=UcAsQKuwWtUlwHjMjn8IGvLRvUbs8O592okHh6M5fuY%3D&reserved=0>)
>
> *A zero pivot in LU, ILU, Cholesky, or ICC sparse factorization does not
> always mean that the matrix is singular. You can use*
>
> *-pc_factor_shift_type nonzero -pc_factor_shift_amount [amount] *
>
> *or*
>
> *-pc_factor_shift_type positive_definite -[level]_pc_factor_shift_type nonzero
> -pc_factor_shift_amount [amount] *
>
> *or*
>
> *-[level]_pc_factor_shift_type positive_definite *
>
> *to prevent the zero pivot. [level] is “sub” when lu, ilu, cholesky, or icc

> are employed in each individual block of the bjacobi or ASM preconditioner.
> [level] is “mg_levels” or “mg_coarse” when lu, ilu, cholesky, or icc are used
> inside multigrid smoothers or to the coarse grid solver. See
> PCFactorSetShiftType

> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Frelease%2Fdocs%2Fmanualpages%2FPC%2FPCFactorSetShiftType.html%23PCFactorSetShiftType&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cff03e4917646439d68e908d9f868a339%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637813950590384693%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=zp3oMzA5JMPuOyVmMPXG3dGx4u1m0nGLcNVh4EmOmLs%3D&reserved=0>(),
> PCFactorSetShiftAmount
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Frelease%2Fdocs%2Fmanualpages%2FPC%2FPCFactorSetShiftAmount.html%23PCFactorSetShiftAmount&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cff03e4917646439d68e908d9f868a339%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637813950590384693%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=F85oRwsykILhaQjMn6lUfkDGQCXuHgl1ExpVllAVLAs%3D&reserved=0>().*
>
> *This error can also happen if your matrix is singular, see MatSetNullSpace
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpetsc.org%2Frelease%2Fdocs%2Fmanualpages%2FMat%2FMatSetNullSpace.html%23MatSetNullSpace&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cff03e4917646439d68e908d9f868a339%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637813950590384693%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=L3uPVPCwyJQF%2FqarLxETmnSrjdq5mCC2MXDXfR7bLCA%3D&reserved=0>() for

> how to handle this. If this error occurs in the zeroth row of the matrix, it

> is likely you have an error in the code that generates the matrix.*

>
> So I write the code as follows:
>

> *PETScWrappers::PreconditionBlockJacobi::AdditionalData data;

>     PETScWrappers::PreconditionBlockJacobi  preconditioner(mpi_communicator,
> data);
>     SolverControl linear_solver_control (dof_handler.n_dofs(),
>                                          1e-9, false, false);
>     PETScWrappers::SolverCG linear_solver(linear_solver_control,
>                                           mpi_communicator);
>
>     PCFactorSetShiftType(preconditioner.get_pc(), MAT_SHIFT_NONZERO);
> //    PCFactorSetShiftAmount(preconditioner.get_pc(), 3e-14);
>
>     linear_solver.initialize(preconditioner);
>
>     SolverControl solver_control (100, 1e-9,false,false);
>     SLEPcWrappers::SolverKrylovSchur eigensolver(solver_control,
>                                                  mpi_communicator);
>
>     eigensolver.set_which_eigenpairs(EPS_SMALLEST_REAL);
>
>     eigensolver.set_problem_type(EPS_GHEP);
>
>     SLEPcWrappers::TransformationShift spectral_transformation(mpi_communicator);
>     spectral_transformation.set_solver(linear_solver);
>     eigensolver.set_transformation(spectral_transformation);
>
>     eigensolver.solve(stiffness_matrix,
>                       mass_matrix,
>                       eigenvalues,
>                       eigenfunctions,

>                       eigenfunctions.size());*

>
> But the error occurs as the same. I don't know how to solve it.
>
> Thanks for your help.
>
> Best Wishes.
>
>
>

> --
> The deal.II project is located at http://www.dealii.org/
> <https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cff03e4917646439d68e908d9f868a339%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637813950590384693%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=AdEUEv%2BUs1X6IcZevlJmxNsSDy3Au7SW9R4Mhow8MAQ%3D&reserved=0>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cff03e4917646439d68e908d9f868a339%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637813950590384693%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=ehjgK4%2BGby6oq7nnOqm%2FqA712gY13A3%2FuBBP%2F%2BVNNQU%3D&reserved=0>
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/c7349e82-26bc-4014-a21d-c52311015921n%40googlegroups.com
> <https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2Fc7349e82-26bc-4014-a21d-c52311015921n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7Cff03e4917646439d68e908d9f868a339%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637813950590384693%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000&sdata=963v4J26aybDjICgrxJ6G%2B8jHiQwNH9TobTj5TIZfLI%3D&reserved=0>.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Léonhard YU]

Professor Bangerth:

I find that if FE_Q is replaced by FESystem in step-36 it reports that error.

I printed all diagonal element of two matrices and they are all non-zero, so the matrices are invertable.

My problem is a 3-dim beam eigenvalue calculation and I use FESystem just like step-18.

And if I use FE_Q rather than FESystem, the function shape_grad_component will report error.

I use that function to calculate stiffness_matrix, called by get_strain() like step-18.

And the theoretical formulae are as follows:

\begin{array}{l}
K = \int {{D^{\rm{T}}}ED{\rm{d}}{V_e}} \\
\varepsilon  = Du\\
\sigma  = \varepsilon E
\end{array}%

I guess that here exits conflict between FESystem and SolverKrylovSchur or just because it is a vector-value problem?

Is it a vector-value the 3-dim eigenvalue calculation?

Thanks for your help.

Best wishes

[Wolfgang Bangerth]

Leonhard,

> I find that if FE_Q is replaced by FESystem in step-36 it reports that
> error.
> I printed all diagonal element of two matrices and they are all
> non-zero, so the matrices are invertable.

That is not the right criterion. The following matrix also has all
nonzero diagonal entries and it still isn't invertible:
[1 1]
[1 1]

> My problem is a 3-dim beam eigenvalue calculation and I use FESystem
> just like step-18.
> And if I use FE_Q rather than FESystem, the function
> shape_grad_component will report error.
> I use that function to calculate stiffness_matrix, called by
> get_strain() like step-18.
> And the theoretical formulae are as follows:
>
> \begin{array}{l}
> K = \int {{D^{\rm{T}}}ED{\rm{d}}{V_e}} \\
> \varepsilon  = Du\\
> \sigma  = \varepsilon E
> \end{array}%
>
> I guess that here exits conflict between FESystem and SolverKrylovSchur
> or just because it is a vector-value problem?
> Is it a vector-value the 3-dim eigenvalue calculation?

There are many things that can go wrong, of course. You need to start
with a simple problem and make incremental changes to find out where the
problem is. I would start with using a bilinear form for the stiffness
matrix that matches what step-36 does, so
(grad u, grad v)
where now u,v are vector valued. Make sure you have zero Dirichlet
boundary conditions for all components. That should work. If it does,
move on to more complicated bilinear forms; with each modification, make
sure that it continues to work, and if it doesn't, you know what
modification was the problem.

Best
W.

[Sebastian D.-R.]

I wonder if there was any follow up replies about this. 

I started to work a few days ago on an eigensolver for linear elasticity with the use of the SLEPc wrappers in deal.II  (based on step-8 and step-36), and I came across the same issue Leonhard described above. After some debugging I realised the issue was in the fact that I was not assembling the mass matrix correctly. Essentially, if following step-8, then you need to only add terms to the local mass matrix whenever comp(i) == comp(j) (following the notation in step-8 documentation), as it is done for the \nabla\phi_i\dot\nabla\phi_j term. 

Another way to go around this issue is to use "FEValuesExtractors::Vector" for the displacement (the only unknown in the system). With this way the assemble looks a lot more closer to how one would write a primal formulation for linear elasticity. So far, both ways give the same eigenvalues, which also make sense (in 2 and 3D, and on different domains).

[Wolfgang Bangerth]

On 10/31/22 22:12, Sebastian D.-R. wrote:
>
> I started to work a few days ago on an eigensolver for linear elasticity
> with the use of the SLEPc wrappers in deal.II  (based on step-8 and
> step-36), and I came across the same issue Leonhard described above.
> After some debugging I realised the issue was in the fact that I was not
> assembling the mass matrix correctly. Essentially, if following step-8,
> then you need to only add terms to the local mass matrix whenever
> comp(i) == comp(j) (following the notation in step-8 documentation), as
> it is done for the \nabla\phi_i\dot\nabla\phi_j term.
>
> Another way to go around this issue is to use
> "FEValuesExtractors::Vector" for the displacement (the only unknown in
> the system). With this way the assemble looks a lot more closer to how
> one would write a primal formulation for linear elasticity. So far, both
> ways give the same eigenvalues, which also make sense (in 2 and 3D, and
> on different domains).

Yes, and that makes sense.

The observation about having to have a check of the form
comp(i)==comp(j), in any form, is correct. If you don't have this, then
the matrix you build is of the form
[M M]
[M M]
which is rank deficient because rows (columns) in the first block of
rows (columns) are identical to rows (columns) in the second block of
rows (columns). So it makes sense that you get a zero pivot: the matrix
really is singular.