如何使用dsyev例程计算特征值? [英] How to use dsyev routine to calculate eigenvalues?
问题描述
int main() {
char jobz='V',uplo='U';
int lda=3,n=3,info=8,lwork=9;
// lapack_int lda=3,n=3,info=8;
int i;
double w[3],work[3];
double a[9] = {
3,2,4,
2,0,2,
4,2,3
};
info=LAPACKE_dsyev(LAPACK_ROW_MAJOR,jobz,uplo, n ,a, lda , w);
//dsyev_( &jobz,&uplo,&n, a, &lda, w,work , &lwork, &info );
if( info > 0 ) {
printf( "The algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
for(i=0;i<3;i++)
{
printf("%f\n",w[i]);
}
for(i=0;i<9;i++)
{
printf("%f\n",a[i]);
}
exit( 0 );
}
输出: -1.000000 -1.000000 8.000000
0.617945 1.999713 -0.016938 0.010468 0.033876 0.999857 1.381966 0.618034 0.000000
我希望 k = -1 :[1，-2,0]，[4,2，-5]和 k = 8 :[2,1,2] 在输出中的某个位置！
我不正确地使用api还是我不正确地读取了输出? 还请建议我如何使用fortran api执行相同的任务? 与fortran一样，我无法获得适当的永恒值！ 即我通过fortran获得的永恒价值: -0.134742 0.050742 0.523036
外来向量: 0.617945 1.999713 -0.016938 0.010468 0.033876 0.999857 1.381966 0.618034 0.000000
如@francis在注释中所建议，如果将work[3]修改为work[9]，该程序将起作用.得到的结果是
Eigenvalues: w[0],w[1],w[2] => -1.000000 -1.000000 8.000000
1st eigenvector: a[0], a[1], a[2] => -0.494101 -0.472019 0.730111
2nd eigenvector: a[3], a[4], a[5] => -0.558050 0.816142 0.149979
3rd eigenvector: a[6], a[7], a[8] => 0.666667 0.333333 0.666667
为了进行比较，让我们用不同的程序对角相同的矩阵.例如，Python/Numpy给出结果
>>> import numpy as np
>>> a = np.array([[3,2,4], [2,0,2], [4,2,3]], dtype=np.float )
>>> np.linalg.eig( a )
(array([-1., 8., -1.]),
array([[-0.74535599, 0.66666667, -0.09414024],
[ 0.2981424 , 0.33333333, -0.84960833],
[ 0.59628479, 0.66666667, 0.5189444 ]]))
朱莉娅奉献时
julia> a = Float64[ 3 2 4 ; 2 0 2 ; 4 2 3 ]
3x3 Array{Float64,2}:
3.0 2.0 4.0
2.0 0.0 2.0
4.0 2.0 3.0
julia> eig( a )
([-0.9999999999999996,-0.9999999999999947,8.0],
3x3 Array{Float64,2}:
0.447214 -0.596285 -0.666667
-0.894427 -0.298142 -0.333333
0.0 0.745356 -0.666667)
在两种情况下，前三个数字均为特征值{-1,-1,8}，随后的矩阵为对应的特征向量(在其列中).您将看到所有程序对本征向量给出不同的结果.由于存在两个简并的特征值(-1)，因此对应特征向量的任何线性组合也是具有相同特征值的特征向量，因此结果不是唯一的.我们可以确定，通过2x2正交变换(或旋转" 101.6度)，从C获得的简并特征向量与Julia简并的特征向量有关.
有趣的是，您期望的特征向量[-1,2,0]和[4,2,-5]完全对应于从Julia获得的特征向量(在归一化之后)，但是这可能是偶然的，并且人们无法期望这样的简并特征向量的一致性.
i am trying to write code to calculate eign vector and eign values for a symmetric matrix. I understand how to calculate evalues using pen & paper but i am slightly confused with the api!. I am a beginner so i may be wrong in interpreting the api parameters.
int main() {
char jobz='V',uplo='U';
int lda=3,n=3,info=8,lwork=9;
// lapack_int lda=3,n=3,info=8;
int i;
double w[3],work[3];
double a[9] = {
3,2,4,
2,0,2,
4,2,3
};
info=LAPACKE_dsyev(LAPACK_ROW_MAJOR,jobz,uplo, n ,a, lda , w);
//dsyev_( &jobz,&uplo,&n, a, &lda, w,work , &lwork, &info );
if( info > 0 ) {
printf( "The algorithm failed to compute eigenvalues.\n" );
exit( 1 );
}
for(i=0;i<3;i++)
{
printf("%f\n",w[i]);
}
for(i=0;i<9;i++)
{
printf("%f\n",a[i]);
}
exit( 0 );
}
output: -1.000000 -1.000000 8.000000
0.617945 1.999713 -0.016938 0.010468 0.033876 0.999857 1.381966 0.618034 0.000000
whereas i expected k=-1: [1,-2,0] ,[4,2,-5] and k=8: [2,1,2] somewhere in the output!
am i using api incorrectly or am i reading the output incorrectly? also please suggest how do i do the same task with fortran api ? as with fortran i am unable to get proper eign values !. i.e. eign values i get with fortran: -0.134742 0.050742 0.523036
eign vectors: 0.617945 1.999713 -0.016938 0.010468 0.033876 0.999857 1.381966 0.618034 0.000000
As suggested by @francis in the comment, the program works if you modify work[3] to work[9]. The obtained result is
Eigenvalues: w[0],w[1],w[2] => -1.000000 -1.000000 8.000000
1st eigenvector: a[0], a[1], a[2] => -0.494101 -0.472019 0.730111
2nd eigenvector: a[3], a[4], a[5] => -0.558050 0.816142 0.149979
3rd eigenvector: a[6], a[7], a[8] => 0.666667 0.333333 0.666667
For comparison, let's diagonalize the same matrix with different programs. For example, Python/Numpy gives the result
>>> import numpy as np
>>> a = np.array([[3,2,4], [2,0,2], [4,2,3]], dtype=np.float )
>>> np.linalg.eig( a )
(array([-1., 8., -1.]),
array([[-0.74535599, 0.66666667, -0.09414024],
[ 0.2981424 , 0.33333333, -0.84960833],
[ 0.59628479, 0.66666667, 0.5189444 ]]))
while Julia gives
julia> a = Float64[ 3 2 4 ; 2 0 2 ; 4 2 3 ]
3x3 Array{Float64,2}:
3.0 2.0 4.0
2.0 0.0 2.0
4.0 2.0 3.0
julia> eig( a )
([-0.9999999999999996,-0.9999999999999947,8.0],
3x3 Array{Float64,2}:
0.447214 -0.596285 -0.666667
-0.894427 -0.298142 -0.333333
0.0 0.745356 -0.666667)
In both cases, the first three numbers are eigenvalues {-1,-1,8} and the following matrix are the corresponding eigenvectors (in their columns). You will see that all the programs give different results for eigenvectors. Because there are two degenerate eigenvalues (-1), any linear combination of the corresponding eigenvectors is also an eigenvector with the same eivenvalue, so the result is not unique. We can confirm that the degenerate eigenvectors obtained from C are related to those of Julia by a 2x2 orthogonal transformation (or a "rotation" by 101.6 degrees).
Interestingly, your expected eigenvectors [-1,2,0] and [4,2,-5] correspond exactly to the eigenvectors obtained from Julia (after normalization), but this is probably accidental and one cannot expect such agreement of degenerate eigenvectors.
这篇关于如何使用dsyev例程计算特征值?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
