LAPACK error in syevx

[Tao Jin]

Dear all,

I received the message "LAPACK error in syevx" when using the compute_eigenvalues_symmetric of lapack_full_matrix (the interface to LAPACK in dealii) to compute the eigenvalues of a symmetric matrix with relatively small matrix entries. I have attached a test file (main.cc) in the end. 

To provide more information, I construct the following 2x2 symmetric matrix

  const double scalar = 1.0;
  for (unsigned int i = 0; i < dim; i++)
    for (unsigned int j = 0; j <= i; j++)
      {
symmetric_tensor[i][j] = (i + j + 1) * scalar;
        if (j != i)
          symmetric_tensor[j][i] = symmetric_tensor[i][j];
      }

In this case, the eigenvalue decomposition runs without error message:

Spectrum decomposition of a symmetric 2nd order tensor
Eigenvalues: -0.236068 4.23607
  Spectrum decomposition succeeded!

However, if I set 

  const double scalar = 1.0e-3; // to make the matrix entries relatively small

The eigenvalue decomposition is still successful. However, here is the output message:

Spectrum decomposition of a symmetric 2nd order tensor
LAPACK error in syevx
Eigenvalues: -0.000236068 0.00423607
  Spectrum decomposition succeeded!

This error message is from "lapack_full_matrix.cc"

// Negative return value implies a wrong argument. This should be internal.

2049 Assert(info >= 0, ExcInternalError());

2050 if (info != 0)

2051 std::cerr << "LAPACK error in syevx" << std::endl;

 

Since I need to run the spectrum decomposition at each Gauss point, it is annoying to receive tons of this error message, even though the simulation still finishes with no issue.

Any insights to avoid this issue?

Best,

Tao

[Wolfgang Bangerth]

Tao:
Yes, this is cumbersome indeed. I looked into this file, and it contains a
number of places where deal.II outputs to std::cerr. This is not the style we
favor today. I looked it up, and these lines were written in 2005, probably
with the very first version of that file.

I have opened an issue about this:
https://github.com/dealii/dealii/issues/15363
As with all things deal.II, we would gladly help any kind of help we can get
-- so if the error message bothers you, perhaps this is your project: Fix this
particular issue, and get it merged into deal.II!

Underlying your question is the issue that dependent on the value of the
scaling, you either do or do not get an error. I believe that that is because
you need to also scale the upper/lower bounds and the tolerance with which you
call the LAPACK function. When you multiply the 2x2 matrix by a scalar factor,
all of the eigenvalues also scale and consequently the tolerance should also
scale.

Best
W.

On 6/14/23 09:47, Tao Jin wrote:
> *** Caution: EXTERNAL Sender ***

>
> Dear all,
>
> I received the message "LAPACK error in syevx" when using the

> /*compute_eigenvalues_symmetric*/ of */lapack_full_matrix/* (the interface to

> LAPACK in dealii) to compute the eigenvalues of a symmetric matrix with
> relatively small matrix entries. I have attached a test file (main.cc) in the
> end.
>
> To provide more information, I construct the following 2x2 symmetric matrix

>   const double *scalar = 1.0;*

>   for (unsigned int i = 0; i < dim; i++)
>     for (unsigned int j = 0; j <= i; j++)
>       {

> symmetric_tensor[i][j] = (i + j + 1) * *scalar*;

>         if (j != i)
>           symmetric_tensor[j][i] = symmetric_tensor[i][j];
>       }
>

> *In this case, the eigenvalue decomposition runs without error message:*

> Spectrum decomposition of a symmetric 2nd order tensor
> Eigenvalues: -0.236068 4.23607
>   Spectrum decomposition succeeded!
>
> However, if I set

>   const double *scalar = 1.0e-3; // to make the matrix entries relatively small*

> The eigenvalue decomposition is still successful. However, here is the output
> message:
> Spectrum decomposition of a symmetric 2nd order tensor

> *LAPACK error in syevx
> *Eigenvalues: -0.000236068 0.00423607

>   Spectrum decomposition succeeded!
>
> This error message is from "lapack_full_matrix.cc"
> // Negative return value implies a wrong argument. This should be internal.
> 2049 Assert

> <https://www.dealii.org/current/doxygen/deal.II/group__Exceptions.html#ga70a0bb353656e704acf927945277bbc6>(info >= 0, ExcInternalError <https://www.dealii.org/current/doxygen/deal.II/group__Exceptions.html#gab7a0d88175320d08084b2f40f5e3380b>());

> 2050 if (info != 0)

> 2051 std::cerr << "*LAPACK error in syevx*" << std::endl;
> **

> Since I need to run the spectrum decomposition at each Gauss point, it is
> annoying to receive tons of this error message, even though the simulation
> still finishes with no issue.
>
> Any insights to avoid this issue?
>
> Best,
>
> Tao
>

> --
> The deal.II project is located at http://www.dealii.org/
> <http://www.dealii.org/>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://groups.google.com/d/forum/dealii?hl=en>
> ---
> 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/8adfe7d5-220c-4c4a-b3e5-60c2c9c6b278n%40googlegroups.com <https://groups.google.com/d/msgid/dealii/8adfe7d5-220c-4c4a-b3e5-60c2c9c6b278n%40googlegroups.com?utm_medium=email&utm_source=footer>.

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

[Tao Jin]

Hi Wolfgang,

Thanks for your input. I would be delighted to look into this issue and potentially make a contribution.

Also, I noticed that SymmetricTensor<2, dim> also has a function called eigenvectors, which could robustly provide the eigenvalues and eigenvectors for a symmetric 2nd order tensor. This function performs well regardless of the scaling. Below is my function to obtain the eigen system of a 2nd order symmetric tensor:

  template <int dim>
  void spectrum_decomposition(SymmetricTensor<2, dim> const & symmetric_tensor,
                                 Vector<double> & myEigenvalues,
                                 std::vector<Tensor<1, dim>> & myEigenvectors)
  {
    std::array< std::pair< double, Tensor< 1, dim > >, dim > myEigenSystem;

    myEigenSystem = eigenvectors(symmetric_tensor);

    for (int i = 0; i < dim; i++)
      {
        myEigenvalues[i] = myEigenSystem[i].first;
        myEigenvectors[i] = myEigenSystem[i].second;
      }
  }

Best,

Tao