Eigenvector issues in 3.12.1

[tweber-ill]

Description
Beginning with Lapack version 3.12.1 I obtain wrong eigenvalues and vectors. So far I've only tested on MacOS Tahoe using an M4 processor and compiling Lapack(e) with Clang 21.1.3 and GFortran 15.2.0, so other systems might not be affected.

• A small C++ test program (see below for the code) yields the following false results using Lapack 3.12.1 (testing matrix * eigenvector gives a different result than eigenvalue * eigenvector, which is obviously wrong):

M = [3,3](((1,0),(3,1.5),(5,2)),((3,-1.5),(2,0),(2.2,-1.7)),((5,-2),(2.2,1.7),(4,0)))

eval0: (7.93151,4.96788), evec0: [3]((0.545641,-0.186362),(0.429787,0.0385125),(0.693791,0)).
eval0 * evec0: [3]((5.25358,1.23254),(3.21754,2.44059),(5.50281,3.44666))
M     * evec0: [3]((5.24619,1.96144),(3.74329,-2.47997),(6.0107,-1.20772))

eval1: (-2.49846,2.29741), evec1: [3]((0.753772,0),(-0.155329,0.396544),(-0.218818,0.450079)).
eval1 * evec1: [3]((-1.88327,1.73173),(-0.522943,-1.34761),(-0.487309,-1.62722))
M     * evec1: [3]((-2.30128,2.7694),(2.23439,1.0246),(1.87774,0.901111))

eval2: (0.396735,0.0246189), evec2: [3]((-0.164256,-0.406065),(0.760621,0),(0.457211,-0.143336)).
eval2 * evec2: [3]((-0.0551693,-0.165144),(0.301765,0.0187256),(0.184921,-0.0456103))
M     * evec2: [3]((4.69033,0.932609),(1.18157,-2.06441),(1.8688,-0.982102))

• Under Lapack 3.12.0 the program gives the following correct output (matrix * eigenvector = eigenvalue * eigenvector):

M = [3,3](((1,0),(3,1.5),(5,2)),((3,-1.5),(2,0),(2.2,-1.7)),((5,-2),(2.2,1.7),(4,0)))

eval0: (10.1646,3.74461e-16), evec0: [3]((0.542088,0.157607),(0.413968,-0.185275),(0.68964,0)).
eval0 * evec0: [3]((5.51011,1.60201),(4.20782,-1.88325),(7.00991,2.58243e-16))
M     * evec0: [3]((5.51011,1.60201),(4.20782,-1.88325),(7.00991,-4.44089e-16))

eval1: (-3.65999,7.47669e-17), evec1: [3]((0.787496,0),(-0.339972,-0.031168),(-0.423307,0.290015)).
eval1 * evec1: [3]((-2.88223,5.88787e-17),(1.2443,0.114075),(1.5493,-1.06145))
M     * evec1: [3]((-2.88223,-2.22045e-16),(1.2443,0.114075),(1.5493,-1.06145))

eval2: (0.495399,-5.13887e-18), evec2: [3]((0.0880892,-0.231064),(0.823256,0),(-0.510609,-0.0194141)).
eval2 * evec2: [3]((0.0436393,-0.114469),(0.40784,-4.23061e-18),(-0.252955,-0.00961772))
M     * evec2: [3]((0.0436393,-0.114469),(0.40784,-4.44089e-16),(-0.252955,-0.00961772))

As far as I could reproduce the error, only complex calculations seem to be affected, not real ones. I've tried compiling different Lapack versions from its git history, and the error seems to have been introduced with pull request #1020.

Here's the source code for my minimal C++ test program (requires Boost):

#include <complex>
#include <iostream>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

namespace ublas = boost::numeric::ublas;
using t_cplx = std::complex<double>;

extern "C"
{
	#define lapack_complex_double t_cplx
	#define lapack_complex_double_real(c) (c.real())
	#define lapack_complex_double_imag(c) (c.imag())

	#include <lapacke.h>
}

bool eigenvec_cplx(const ublas::matrix<t_cplx>& mat,
	std::vector<ublas::vector<t_cplx>>& evecs,
	std::vector<t_cplx>& evals)
{
	const std::size_t order = mat.size1();
	evecs.resize(order);
	evals.resize(order);
	for(std::size_t i = 0; i < order; ++i)
		evecs[i].resize(order);

	// allocate buffer for matrix and eigenvectors
	std::unique_ptr<t_cplx, std::default_delete<t_cplx[]>>
		buf(new t_cplx[2*order*order]);
	t_cplx *buf_matrix = buf.get();
	t_cplx *buf_evecs = buf.get() + order*order;

	// copy matrix to buffer
	for(std::size_t i = 0; i < order; ++i)
		for(std::size_t j = 0; j < order; ++j)
			buf_matrix[i*order + j] = mat(i, j);

	// calculate eigenvectors and -values
	int iInfo = LAPACKE_zgeev(LAPACK_ROW_MAJOR, 'N', 'V', order,
		buf_matrix, order, evals.data(), nullptr, order, buf_evecs, order);
	if(iInfo != 0)
	{
		std::cerr << "Lapack error: " << iInfo << "." << std::endl;
		return false;
	}

	// copy buffer to eigenvectors
	for(std::size_t i = 0; i < order; ++i)
		for(std::size_t j = 0; j < order; ++j)
			evecs[i][j] = buf_evecs[j*order + i];

	return true;
}

int main()
{
	// example matrix
	ublas::matrix<t_cplx> M(3, 3);
	M(0, 0) = t_cplx(1., 0.);
	M(0, 1) = t_cplx(3., 1.5);
	M(0, 2) = t_cplx(5., 2.);
	M(1, 0) = t_cplx(3., -1.5);
	M(1, 1) = t_cplx(2., 0.);
	M(1, 2) = t_cplx(2.2, -1.7);
	M(2, 0) = t_cplx(5., -2.);
	M(2, 1) = t_cplx(2.2, 1.7);
	M(2, 2) = t_cplx(4., 0.);	
	std::cout << "M = " << M << "\n" << std::endl;

	std::vector<ublas::vector<t_cplx>> evecs;
	std::vector<t_cplx> evals;

	if(!eigenvec_cplx(M, evecs, evals))
		return -1;

	for(std::size_t i = 0; i < evals.size(); ++i)
	{
		std::cout << "eval" << i << ": " << evals[i]
			<< ", evec" << i << ": " << evecs[i] << ".\n"
			<< "eval" << i << " * evec" << i << ": " << evals[i]*evecs[i]
			<< "\nM     * evec" << i << ": " << ublas::prod(M, evecs[i])
			<< "\n" << std::endl;
	}

	return 0;
}

Checklist

• I've included a minimal example to reproduce the issue
• I'd be willing to make a PR to solve this issue