zheev дает неправильные собственные значения (проверено на zgeev и numpy.linalg.eig)
Работая с библиотеками linalg, я попытался настроить подпрограмму диагонализации для эрмитовых матриц и запустить ее на C++, используя LAPACKE.
Я следовал этому примеру, чтобы использовать ZHEEV, а затем проверил некоторые другие методы, в частности, numpy's eig и LAPACK(E) zgeev. Я не хотел использовать материал от собственной марки Intel, поэтому я избегал MKL и просто выбрал LAPACKE, но большая часть кода такая же, как в примере.
Для ясности, я не вижу причин, по которым общий ge ge ev не сможет обрабатывать конкретный случай гермитовой матрицы, даже если z he ev оптимизирован.
Вот с ++
#include <stdlib.h>
#include <stdio.h>
#include <lapacke.h>
//Parameters
#define N 4
#define LDA N
#define lint lapack_int
#define ldcmplex lapack_complex_double
//Auxiliary routines prototypes
extern void print_matrix( char* desc, lint m, lint n, ldcmplex* a, lint lda );
extern void print_rmatrix( char* desc, lint m, lint n, double* a, lint lda );
//Main program
int main()
{
//Locals
lint n = N, lda = LDA, info;
;
//Local arrays
double wr[N];
ldcmplex ah[LDA*N] = {
{ 9.14, 0.00}, { 0.00, 0.00}, { 0.00, 0.00}, { 0.00, 0.00},
{-4.37, 9.22}, {-3.35, 0.00}, { 0.00, 0.00}, { 0.00, 0.00},
{-1.98, 1.72}, { 2.25, 9.51}, {-4.82, 0.00}, { 0.00, 0.00},
{-8.96, 9.50}, { 2.57, -2.40}, {-3.24, -2.04}, { 8.44, 0.00}
};
;
//Executable statements
printf( "LAPACKE_zheev (row-major, high-level) Example Program Results\n" ) ;
;
//Print martix
print_matrix( "Input Matrix", n, n, ah, lda );
;
//Solve eigenproblem
info = LAPACKE_zheev( LAPACK_ROW_MAJOR, 'V', 'L', n, ah, lda, wr );
;
//Check for convergence
if( info > 0 ) {
printf( "zheev algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
;
//Print eigenvalues
print_rmatrix( "zheev Eigenvalues", 1, n, wr, 1 );
;
//Print eigenvectors
print_matrix( "Eigenvectors (stored columnwise)", n, n, ah, lda );
;
//Local arrays
ldcmplex wc[N];
ldcmplex ag[LDA*N] = {
{ 9.14, 0.00}, {-4.37, -9.22}, {-1.98, -1.72}, {-8.96, -9.50},
{-4.37, 9.22}, {-3.35, 0.00}, { 2.25, -9.51}, { 2.57, 2.40},
{-1.98, 1.72}, { 2.25, 9.51}, {-4.82, 0.00}, {-3.24, 2.04},
{-8.96, 9.50}, { 2.57, -2.40}, {-3.24, -2.04}, { 8.44, 0.00},
};
;
//Executable statements
printf( "LAPACKE_zgeev (row-major, high-level) Example Program Results\n" );
;
//Print martix
print_matrix( "Input Matrix", n, n, ag, lda );
;
//Solve eigenproblem
info = LAPACKE_zgeev( LAPACK_ROW_MAJOR, 'N', 'V', n, ag, lda, wc, 0, lda, 0, lda);
;
//Check for convergence
if( info > 0 ) {
printf( "zgeev algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
;
//Print eigenvalues
print_matrix( "zgeev Eigenvalues", 1, n, wc, 1);
;
//Print eigenvectors
print_matrix( "Eigenvectors (stored columnwise)", n, n, ag, lda );
exit( 0 );
}
//Auxiliary routine: printing a matrix
void print_matrix( char* desc, lint m, lint n, ldcmplex* a, lint lda ) {
lint i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ )
printf( " (%6.2f,%6.2f)", creal(a[i*lda+j]), cimag(a[i*lda+j]) );
printf( "\n" );
}
}
//Auxiliary routine: printing a real matrix
void print_rmatrix( char* desc, lint m, lint n, double* a, lint lda ) {
lint i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %6.2f", a[i*lda+j] );
printf( "\n" );
}
}
составлено с
g++ diag.cc -L /usr/lib/lapack/ -llapacke -lcblas -o diag.out
Единственные нестандартные liblapacke-dev, а также libcblas-dev устанавливается через apt-get install, Что возможно могло пойти не так?
Выход
LAPACKE_zheev (row-major, high-level) Example Program Results
Input Matrix
( 9.14, 0.00) ( 0.00, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( -4.37, 9.22) ( -3.35, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( -1.98, 1.72) ( 2.25, 9.51) ( -4.82, 0.00) ( 0.00, 0.00)
( -8.96, 9.50) ( 2.57, -2.40) ( -3.24, -2.04) ( 8.44, 0.00)
zheev Eigenvalues
-18.96 -12.85 18.78 30.71
Eigenvectors (stored columnwise)
( 0.16, 0.00) ( 0.57, 0.00) ( -0.73, 0.00) ( 0.35, 0.00)
( 0.26, -0.81) ( 0.17, -0.25) ( 0.22, -0.38) ( 0.06, -0.02)
( 0.29, 0.27) ( -0.11, -0.30) ( -0.26, -0.42) ( -0.50, -0.50)
( -0.21, 0.23) ( 0.50, -0.49) ( 0.18, -0.09) ( -0.33, 0.51)
LAPACKE_zgeev (row-major, high-level) Example Program Results
Input Matrix
( 9.14, 0.00) ( -4.37, -9.22) ( -1.98, -1.72) ( -8.96, -9.50)
( -4.37, 9.22) ( -3.35, 0.00) ( 2.25, -9.51) ( 2.57, 2.40)
( -1.98, 1.72) ( 2.25, 9.51) ( -4.82, 0.00) ( -3.24, 2.04)
( -8.96, 9.50) ( 2.57, -2.40) ( -3.24, -2.04) ( 8.44, 0.00)
zgeev Eigenvalues
( 25.51, 0.00) (-16.00, -0.00) ( -6.76, 0.00) ( 6.67, 0.00)
Eigenvectors (stored columnwise)
( 25.51, 0.00) ( -0.00, 0.00) ( 0.00, 0.00) ( 0.00, -0.00)
( 0.00, 0.00) (-16.00, -0.00) ( 0.00, -0.00) ( 0.00, -0.00)
( 0.00, 0.00) ( 0.00, 0.00) ( -6.76, 0.00) ( -0.00, -0.00)
( 0.00, 0.00) ( 0.00, 0.00) ( 0.00, 0.00) ( 6.67, 0.00)
Я попытался использовать верхний треугольник, заполнить матрицу и различные другие исправления. Те же результаты каждый раз.
Я с подозрением отношусь к #define ldcmplex lapack_complex_double макрос, но вся документация, которую я могу найти, говорит, что я должен использовать двойной комплекс, так что я немного растерялся. Во всяком случае, если бы это была проблема, почему бы работал zgeev?
Во всяком случае, вот скрипт проверки Python:
#!/usr/bin/env python
from numpy import linalg as li
import numpy as np
mat=np.array([
[ 9.14 + 0.00j, 0.00 + 0.00j, 0.00 + 0.00j, 0.00 +0.00j],
[ -4.37 + 9.22j, -3.35 + 0.00j, 0.00 + 0.00j, 0.00 +0.00j],
[ -1.98 + 1.72j, 2.25 + 9.51j, -4.82 + 0.00j, 0.00 +0.00j],
[ -8.96 + 9.50j, 2.57 - 2.40j, -3.24 - 2.04j, 8.44 +0.00j]])
mat[0]=np.conj(mat[:,0])
mat[1]=np.conj(mat[:,1])
mat[2]=np.conj(mat[:,2])
mat[3]=np.conj(mat[:,3])
mat=np.matrix(mat)
w, v = li.eig(mat)
print w
print v
Согласен с zgeev (до некоторых ошибок округления / машин). Результат также подтверждается указанным выше учебником Intel. Метод Жеева явно в меньшинстве одного, я просто не знаю почему.
Я пробовал это на нескольких машинах:
Linux parabox 4.8.0-52-generic #55-Ubuntu SMP Fri Apr 28 13:28:50 UTC 2017 x86_64 x86_64 x86_64 GNU/Linux
Linux glass 4.10.0-21-generic #23-Ubuntu SMP Fri Apr 28 16:14:22 UTC 2017 x86_64 x86_64 x86_64 GNU/Linux
Любая помощь приветствуется.
3 ответа
Замена -cblas с -blas в строке компиляции решает проблему.
В пакете cblas должна быть ошибка.
Вот что я получаю, когда запускаю ваш скрипт на python:
$ ./diag.py
[ 25.51400517 +1.20330583e-15j -16.00474647 -2.91871119e-15j
-6.76497015 -6.59630730e-16j 6.66571145 +1.54590036e-16j]
[[ 0.69747489+0.j 0.21857873+0.26662122j 0.47736933+0.26449375j -0.02829581-0.30988089j]
[-0.21578745+0.28003172j 0.69688890+0.j -0.14143627-0.2852389j 0.24437193-0.47778739j]
[-0.14607303-0.08302697j -0.01445974-0.60818924j 0.44678181+0.26546077j 0.57583205+0.j ]
[-0.34133591+0.49376693j 0.15930699-0.00061647j 0.57507627+0.j -0.45823952+0.2713093j ]]
Я не знаю, что должно соответствовать. Собственные значения совпадают, но не собственные векторы.
Мне просто пришлось иметь дело с
ZGEEV()и наткнулся на этот вопрос. В этой довольно старой системе (Ubuntu 14.04) с LAPACK, поддерживаемым BLAS 3, результаты этого образца также не соответствовали ожидаемым.
Вот и оказалось, что в моем случае виновата инициализация матрицы. Вместо того, чтобы делать его встроенным, как в исходном примере, я решил прочитать его в
double[2] введите, а затем повторно инициализируйте его в
lapack_complex_double используя добавленные
set_matrix() функция.
Кроме того, для правильного вывода собственных векторов из
wgr в него должен быть передан массив, а затем распечатан.
Ниже представлен обновленный исходный код и его вывод. Надеюсь, это будет полезно для кого-то еще, исследующего
LAPACKE_zgeev а также
LAPACKE_zheev функции (еще не так много документации).
#include <stdlib.h>
#include <stdio.h>
#include <lapacke.h>
//Parameters
#define N 4
#define LDA N
#define lint lapack_int
#define ldcmplex lapack_complex_double
typedef struct double2 {
double v[2];
} double2_t;
//Auxiliary routines prototypes
extern void print_matrix( char* desc, lint m, lint n, ldcmplex* a, lint lda);
extern void print_rmatrix( char* desc, lint m, lint n, double* a, lint lda);
extern void set_matrix(lapack_int n, lapack_complex_double* a, lapack_int lda, double2_t *a2);
//Main program
int main()
{
//Locals
lint n = N, lda = LDA, info;
;
//Local arrays
double wr[N];
ldcmplex ah[LDA*N] = {0};
double2_t ah2[LDA*N] = {
{ 9.14, 0.00}, { 0.00, 0.00}, { 0.00, 0.00}, { 0.00, 0.00},
{-4.37, 9.22}, {-3.35, 0.00}, { 0.00, 0.00}, { 0.00, 0.00},
{-1.98, 1.72}, { 2.25, 9.51}, {-4.82, 0.00}, { 0.00, 0.00},
{-8.96, 9.50}, { 2.57, -2.40}, {-3.24, -2.04}, { 8.44, 0.00}
};
;
//Executable statements
set_matrix(n, ah, lda, ah2);
printf( "LAPACKE_zheev (row-major, high-level) Example Program Results\n" ) ;
;
//Print martix
print_matrix( "Input Matrix", n, n, ah, lda );
;
//Solve eigenproblem
info = LAPACKE_zheev( LAPACK_ROW_MAJOR, 'V', 'L', n, ah, lda, wr );
;
//Check for convergence
if( info > 0 ) {
printf( "zheev algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
;
//Print eigenvalues
print_rmatrix( "zheev Eigenvalues", 1, n, wr, 1 );
;
//Print eigenvectors
print_matrix( "Eigenvectors (stored columnwise)", n, n, ah, lda );
;
//Local arrays
ldcmplex wc[N];
ldcmplex ag[LDA*N] = {0};
ldcmplex wgr[LDA*N] = {0};
double2_t ag2[LDA*N] = {
{ 9.14, 0.00}, {-4.37, -9.22}, {-1.98, -1.72}, {-8.96, -9.50},
{-4.37, 9.22}, {-3.35, 0.00}, { 2.25, -9.51}, { 2.57, 2.40},
{-1.98, 1.72}, { 2.25, 9.51}, {-4.82, 0.00}, {-3.24, 2.04},
{-8.96, 9.50}, { 2.57, -2.40}, {-3.24, -2.04}, { 8.44, 0.00},
};
;
printf("\n\n");
//Executable statements
set_matrix(n, ag, lda, ag2);
printf( "LAPACKE_zgeev (row-major, high-level) Example Program Results\n" );
;
//Print martix
print_matrix( "Input Matrix", n, n, ag, lda );
;
//Solve eigenproblem
info = LAPACKE_zgeev( LAPACK_ROW_MAJOR, 'N', 'V', n, ag, lda, wc, 0, lda, wgr, lda);
;
//Check for convergence
if( info > 0 ) {
printf( "zgeev algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
;
//Print eigenvalues
print_matrix( "zgeev Eigenvalues", 1, n, wc, 1);
;
//Print eigenvectors
print_matrix( "Eigenvectors (stored columnwise)", n, n, wgr, lda );
exit( 0 );
}
//Auxiliary routine: printing a matrix
void print_matrix( char* desc, lint m, lint n, ldcmplex* a, lint lda ) {
lint i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ )
printf( " (%6.2f,%6.2f)", creal(a[i*lda+j]), cimag(a[i*lda+j]) );
printf( "\n" );
}
}
//Auxiliary routine: printing a real matrix
void print_rmatrix( char* desc, lint m, lint n, double* a, lint lda ) {
lint i, j;
printf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) printf( " %6.2f", a[i*lda+j] );
printf( "\n" );
}
}
//Auxiliary routine: set a complex matrix from a double[2] type matrix
void set_matrix(lapack_int n, lapack_complex_double* a, lapack_int lda, double2_t *a2) {
lapack_int i, j;
for( i = 0; i < n; i++ ) {
for( j = 0; j < n; j++ )
a[i*lda+j] = lapack_make_complex_double(a2[i*lda+j].v[0], a2[i*lda+j].v[1]);
}
}
Сборка и запуск:
gcc sample.c -llapacke && ./a.out. Выход:
LAPACKE_zheev (row-major, high-level) Example Program Results
Input Matrix
( 9.14, 0.00) ( 0.00, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( -4.37, 9.22) ( -3.35, 0.00) ( 0.00, 0.00) ( 0.00, 0.00)
( -1.98, 1.72) ( 2.25, 9.51) ( -4.82, 0.00) ( 0.00, 0.00)
( -8.96, 9.50) ( 2.57, -2.40) ( -3.24, -2.04) ( 8.44, 0.00)
zheev Eigenvalues
-16.00 -6.76 6.67 25.51
Eigenvectors (stored columnwise)
( 0.34, -0.00) ( -0.55, 0.00) ( 0.31, 0.00) ( -0.70, 0.00)
( 0.44, -0.54) ( 0.26, 0.18) ( 0.45, 0.29) ( 0.22, -0.28)
( -0.48, -0.37) ( -0.52, -0.02) ( -0.05, 0.57) ( 0.15, 0.08)
( 0.10, -0.12) ( -0.50, 0.28) ( -0.23, -0.48) ( 0.34, -0.49)
LAPACKE_zgeev (row-major, high-level) Example Program Results
Input Matrix
( 9.14, 0.00) ( -4.37, -9.22) ( -1.98, -1.72) ( -8.96, -9.50)
( -4.37, 9.22) ( -3.35, 0.00) ( 2.25, -9.51) ( 2.57, 2.40)
( -1.98, 1.72) ( 2.25, 9.51) ( -4.82, 0.00) ( -3.24, 2.04)
( -8.96, 9.50) ( 2.57, -2.40) ( -3.24, -2.04) ( 8.44, 0.00)
zgeev Eigenvalues
( 25.51, -0.00) (-16.00, -0.00) ( -6.76, -0.00) ( 6.67, -0.00)
Eigenvectors (stored columnwise)
( 0.70, 0.00) ( 0.22, 0.27) ( 0.48, 0.26) ( -0.03, -0.31)
( -0.22, 0.28) ( 0.70, 0.00) ( -0.14, -0.29) ( 0.24, -0.48)
( -0.15, -0.08) ( -0.01, -0.61) ( 0.45, 0.27) ( 0.58, 0.00)
( -0.34, 0.49) ( 0.16, -0.00) ( 0.58, 0.00) ( -0.46, 0.27)
Как можно заметить, собственные значения и собственные векторы от обоих
ZHEEV а также
ZGEEVдействительно одинаковы, только их порядок разный. Поэтому важно убедиться, что исходная матрица сформирована правильно, поскольку данные в конечном итоге передаются в Фортран по ссылке ('мусор на входе -> мусор на выходе').