您在zheev() 上链接的示例已经使用了未打包的LDA*N=N*N 矩阵。实际上,厄米矩阵不需要打包:您可能不必复制矩阵。注意:zheev()修改矩阵A!
LAPACK 处理矩阵的其他存储模式:参见 LAPACK 的the naming scheme。例如:
-
zheev():内存占用N*N 和存储类似于一般解压的N*N 矩阵之一。根据参数UPLO 的值,使用或忽略上三角部分。无论如何,矩阵可以被填充,就好像它是一个大小为N*N 的一般解包矩阵。在这种情况下,参数UPLO 的值不应改变获得的特征值。
-
zhpev():打包格式。根据参数UPLO 的值,存储上对角线项或下对角线项。矩阵存储的内存占用为(N*(N+1))/2。
-
zhbev():专用于频段存储。
当您使用 C 或 C++ 时,这里是通过接口 LAPACKE 使用 zheev() 的示例代码。可以通过gcc main.c -o main -llapacke -llapack -lblas -lm -Wall编译。此外,此代码确保函数zheev() 返回右特征向量,而不是左特征向量。如here 所述,左侧特征向量是右侧特征向量的复共轭。
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#include <time.h>
#include <lapacke.h>
int main(void){
int n=200;
srand(time(NULL));
// allocate the matrix, unpacked storage
double complex** A=malloc(n*sizeof(double complex*));
if(A==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
A[0]=malloc(n*n*sizeof(double complex));
if(A[0]==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
int i;
for(i=1;i<n;i++){
A[i]=&A[0][n*i];
}
//populte the matrix, only the lower diagonal part
int j;
for(i=0;i<n;i++){
for(j=0;j<i;j++){
A[i][j]=rand()/((double)RAND_MAX)+rand()/((double)RAND_MAX)*I;
}
A[i][i]=rand()/((double)RAND_MAX)+42.;
}
//saving the matrix, because zheev() changes it.
double complex** As=malloc(n*sizeof(double complex*));
if(As==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
As[0]=malloc(n*n*sizeof(double complex));
if(As[0]==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
for(i=1;i<n;i++){
As[i]=&As[0][n*i];
}
for(i=0;i<n;i++){
for(j=0;j<i;j++){
As[i][j]=A[i][j];
}
As[i][i]=A[i][i];
}
//transpose part, conjugate
for(i=0;i<n;i++){
for(j=i+1;j<n;j++){
As[i][j]=conj(A[j][i]);
}
}
/*
for(i=0;i<n;i++){
for(j=0;j<n;j++){
printf("%g+I%g ",creal(As[i][j]),cimag(As[i][j]));
}
printf("\n");
}
*/
//let's get the eigenvalues and the eigenvectors!
double* w=malloc(n*sizeof(double));
if(w==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
int lda = n;
int info = LAPACKE_zheev(LAPACK_ROW_MAJOR, 'V', 'L', n, A[0], lda, w);
if(info<0){
fprintf(stderr,"argument %d has illegal value\n",-info);
}
if(info>0){
fprintf(stderr,"algorithm failed to converge... bad condition number ?\n");
}
//printing the eigenvalues...
for(i=0;i<n;i++){
printf("%d %g\n",i,w[i]);
}
//let's check that the column i of A is now a right eigenvectors corresponding to the eigenvalue w[i]...
int l=42;
double complex *res=malloc(n*sizeof(double complex));
if(res==NULL){fprintf(stderr,"malloc failed\n");exit(1);}
for(i=0;i<n;i++){
res[i]=0;
for(j=0;j<n;j++){
res[i]+=As[i][j]*A[j][l];
}
res[i]-=w[l]*A[i][l];
}
double norm2res=0;
double norm2vec=0;
for(i=0;i<n;i++){
norm2res+=creal(res[i])*creal(res[i])+cimag(res[i])*cimag(res[i]);
norm2vec+=creal(A[i][l])*creal(A[i][l])+cimag(A[i][l])*cimag(A[i][l]);
}
printf("the norm of the eigenvector is %g\n",sqrt(norm2vec));
printf("||Ax-\\lambda*x||_2/||x||_2 is about %g\n",sqrt(norm2res/norm2vec));
//free the matrix
free(A[0]);
free(A);
free(As[0]);
free(As);
free(w);
free(res);
return 0;
}
在上面的代码中,执行了矩阵的副本,但LAPACKE_zheev() 不需要这样做。处理一个2000*2000的矩阵,上面代码的内存占用大约是167MB。这是矩阵大小 (64MB) 的两倍多,因为执行了复制。但如果不进行复制,则不到两倍。因此,LAPACKE_zheev() 不执行任何矩阵副本。另请注意,LAPACKE_zheev() 为工作数组分配了一些空间。