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/
For mailing list/forum options, see
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/8adfe7d5-220c-4c4a-b3e5-60c2c9c6b278n%40googlegroups.com.
##
# CMake script for the linear thermoelastic monolithic approach:
##
# Set the name of the project and target:
SET(TARGET "main")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 3.3.0)
FIND_PACKAGE(deal.II 9.4.0
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
/* ---------------------------------------------------------------------
*
* Copyright (C) 2006 - 2020 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE.md at
* the top level directory of deal.II.
*
* ---------------------------------------------------------------------
*
* Author: Tao Jin
* University of Ottawa, Ottawa, Ontario, Canada
* May 2023
*/
// In this example, the spectrum decomposition of a symmetric 2nd order tensor
// and its relevant derivatives are calculated.
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/patterns.h>
#include <deal.II/physics/elasticity/standard_tensors.h>
#include <fstream>
#include <iostream>
#include<ctime>
template <int dim>
void spectrum_decomposition_test()
{
using namespace dealii;
SymmetricTensor<2, dim> symmetric_tensor;
const double scalar = 1.0e-3;
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];
}
Vector<double> eigenvalues(dim);
std::vector<Tensor<1, dim>> eigenvectors(dim);
FullMatrix<double> symmetric_matrix(dim, dim);
for (int i = 0; i < dim; i++)
for (int j = 0; j < dim; j++)
{
symmetric_matrix(i,j) = symmetric_tensor[i][j];
}
//std::ofstream myMatrix("mySymmetricMatrix.matrix");
//symmetric_matrix.print_formatted(myMatrix);
// create a LapackMatrix from full matrix to use the eigenvalue functions
// in Lapack
LAPACKFullMatrix<double> LapackMatrix(dim);
LapackMatrix = symmetric_matrix;
FullMatrix<double> eigenvectors_fullmatrix;
LapackMatrix.compute_eigenvalues_symmetric(Patterns::Double::min_double_value,
Patterns::Double::max_double_value,
1.0e-16,
eigenvalues,
eigenvectors_fullmatrix);
//std::ofstream output_eigenvalues("eigenvalues.matrix");
//eigenvalues.print(output_eigenvalues);
//std::ofstream output_eigenvectors("eigenvectors.matrix");
//eigenvectors_fullmatrix.print_formatted(output_eigenvectors);
for (int i = 0; i < dim; i++)
{
for (int j = 0; j < dim; j++)
eigenvectors[i][j] = eigenvectors_fullmatrix(j, i);
}
std::cout << "Eigenvalues: " << eigenvalues << std::endl;
// reconstruct original matrix via spectrum decomposition
Tensor<2, dim> reconstructed_tensor;
for (unsigned int i = 0; i < dim; i++)
reconstructed_tensor += eigenvalues[i] * outer_product(eigenvectors[i],
eigenvectors[i]);
bool spetrum_decomposition_success = true;
if ( (reconstructed_tensor - symmetric_tensor).norm() > 1.0e-12 )
spetrum_decomposition_success = false;
for (unsigned int i = 0; i < dim; i++)
if ( ( reconstructed_tensor * eigenvectors[0]
- eigenvalues[0] * eigenvectors[0]).norm() > 1.0e-12 )
spetrum_decomposition_success = false;
if (spetrum_decomposition_success)
std::cout << " Spectrum decomposition succeeded!" << std::endl;
else
Assert(false,
ExcMessage(" Spectrum decomposition failed!"));
}
int main()
{
std::cout << "Spectrum decomposition of a symmetric 2nd order tensor" << std::endl;
spectrum_decomposition_test<2>();
return 0;
}
