GSL 复数矩阵 - eigenvalues/eigenvectors
GSL complex matrix - eigenvalues/eigenvectors
我编写了一个程序来计算厄密矩阵的特征值和特征向量。
有谁知道这在 GSL 中是如何正确完成的?这是我已有的。
//hermitian matrix
0 1 0 -i
1 0 -i 0
0 i 0 1
i 0 1 0
#include <iostream>
#include <stdio.h>
#include <cmath>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
using namespace std;
const int N = 4;
int main(){
gsl_eigen_hermv_workspace *workN = gsl_eigen_hermv_alloc(N);
gsl_matrix_complex *A = gsl_matrix_complex_alloc(N, N);
gsl_complex i = gsl_complex_rect(0.0,1.0);
gsl_complex ii = gsl_complex_rect(0.0,-1.0);
gsl_vector *eval = gsl_vector_alloc(N);
gsl_matrix_complex *evec = gsl_matrix_complex_alloc(N, N);
double mTab[] = {
0, 1, 0, 5,
1, 0, 5, 0,
0, 5, 0, 1,
5, 0, 1, 0
};
gsl_matrix_complex_view tmpM = gsl_matrix_complex_view_array(mTab, N, N);
gsl_matrix_complex_memcpy(A, &tmpM.matrix);
gsl_matrix_complex_set(A, 0, 3, ii);
gsl_matrix_complex_set(A, 1, 2, ii);
gsl_matrix_complex_set(A, 2, 1, i);
gsl_matrix_complex_set(A, 3, 0, i);
gsl_eigen_hermv(A, eval, evec, workN);
for(int i=0; i < N; i++){
for(int j=0; j < N; j++){
gsl_complex z = gsl_matrix_complex_get(A, i, j);
cout << GSL_REAL(z) << "+ i" << GSL_IMAG(z) << " ";
}
cout << "\n";
}
cout << "\n";
for(int i=0; i < N; i++){
cout << gsl_vector_get(eval, i) << " ";
}
return 0;
}
这就是我输出特征向量的方式
for(int i=0; i < N; i++){
for(int j=0; j < N; j++){
gsl_complex z = gsl_matrix_complex_get(A, i, j);
cout << GSL_REAL(z) << "+ i" << GSL_IMAG(z) << " ";
}
cout << "\n";
}
最后,这是我声明相关矩阵的方式。
double mTab[] = {
0, 1, 0, 5,
1, 0, 5, 0,
0, 5, 0, 1,
5, 0, 1, 0
};
后来,我加了复数。
我设法打印了特征向量,但我不知道如何打印特征值。对此有任何帮助吗?
您在将 double mTab
用于 gsl_matrix_complex_view_array
时出错。这假设一个复数数组表示为实部,后跟虚部,在一个 double
值的大数组中。您可以将定义更改为:
double mTab[] = {
0, 0, 1, 0, 0, 0, 5, 0,
1, 0, 0, 0, 5, 0, 0, 0,
0, 0, 5, 0, 0, 0, 1, 0,
5, 0, 0, 0, 1, 0, 0, 0,
};
(这也意味着你不需要使用 5 的 "dummy" 变量只是为了稍后用 ±i 重写它。)然后打印特征值的代码运行良好。
你在特征向量打印循环中也有错别字:它应该是
gsl_complex z = gsl_matrix_complex_get(evec, i, j);
没有
gsl_complex z = gsl_matrix_complex_get(A, i, j);
我编写了一个程序来计算厄密矩阵的特征值和特征向量。
有谁知道这在 GSL 中是如何正确完成的?这是我已有的。
//hermitian matrix
0 1 0 -i
1 0 -i 0
0 i 0 1
i 0 1 0
#include <iostream>
#include <stdio.h>
#include <cmath>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_eigen.h>
using namespace std;
const int N = 4;
int main(){
gsl_eigen_hermv_workspace *workN = gsl_eigen_hermv_alloc(N);
gsl_matrix_complex *A = gsl_matrix_complex_alloc(N, N);
gsl_complex i = gsl_complex_rect(0.0,1.0);
gsl_complex ii = gsl_complex_rect(0.0,-1.0);
gsl_vector *eval = gsl_vector_alloc(N);
gsl_matrix_complex *evec = gsl_matrix_complex_alloc(N, N);
double mTab[] = {
0, 1, 0, 5,
1, 0, 5, 0,
0, 5, 0, 1,
5, 0, 1, 0
};
gsl_matrix_complex_view tmpM = gsl_matrix_complex_view_array(mTab, N, N);
gsl_matrix_complex_memcpy(A, &tmpM.matrix);
gsl_matrix_complex_set(A, 0, 3, ii);
gsl_matrix_complex_set(A, 1, 2, ii);
gsl_matrix_complex_set(A, 2, 1, i);
gsl_matrix_complex_set(A, 3, 0, i);
gsl_eigen_hermv(A, eval, evec, workN);
for(int i=0; i < N; i++){
for(int j=0; j < N; j++){
gsl_complex z = gsl_matrix_complex_get(A, i, j);
cout << GSL_REAL(z) << "+ i" << GSL_IMAG(z) << " ";
}
cout << "\n";
}
cout << "\n";
for(int i=0; i < N; i++){
cout << gsl_vector_get(eval, i) << " ";
}
return 0;
}
这就是我输出特征向量的方式
for(int i=0; i < N; i++){
for(int j=0; j < N; j++){
gsl_complex z = gsl_matrix_complex_get(A, i, j);
cout << GSL_REAL(z) << "+ i" << GSL_IMAG(z) << " ";
}
cout << "\n";
}
最后,这是我声明相关矩阵的方式。
double mTab[] = {
0, 1, 0, 5,
1, 0, 5, 0,
0, 5, 0, 1,
5, 0, 1, 0
};
后来,我加了复数。
我设法打印了特征向量,但我不知道如何打印特征值。对此有任何帮助吗?
您在将 double mTab
用于 gsl_matrix_complex_view_array
时出错。这假设一个复数数组表示为实部,后跟虚部,在一个 double
值的大数组中。您可以将定义更改为:
double mTab[] = {
0, 0, 1, 0, 0, 0, 5, 0,
1, 0, 0, 0, 5, 0, 0, 0,
0, 0, 5, 0, 0, 0, 1, 0,
5, 0, 0, 0, 1, 0, 0, 0,
};
(这也意味着你不需要使用 5 的 "dummy" 变量只是为了稍后用 ±i 重写它。)然后打印特征值的代码运行良好。
你在特征向量打印循环中也有错别字:它应该是
gsl_complex z = gsl_matrix_complex_get(evec, i, j);
没有
gsl_complex z = gsl_matrix_complex_get(A, i, j);