我使用gcc在mac os x下编译,我安装了Intel的mkl_lapack.h库。
在程序中,我有一个NxN三对角矩阵,所以我只使用两个向量来存储矩阵的值。
“ d”向量是主要对角线,次对角线的值存储在“ e”中。
首先,我初始化值,然后由于矩阵为16x16(在输入中,我将值16设为argv [1]),因此我将向量分为两个向量(我只能一次使用dstev,但这是为了实验目的),从d [0]到d [N / 2-1]我有第一个向量,从d [N / 2]到d [N-1]我有第二个向量。
因此,一旦初始化“ e”和“ d”的值,我将调用两次dstev。
但是我不必费心在“ z”中写入所有值(z将包含本征向量),因为我知道在两次调用dstev之后,在所有“ z”向量中我应该只有两个值的子矩阵,即8x8的非子矩阵。 -零值。
但是,如果我尝试使用“ z”,则某些值为0.0,并且无法解释为什么会发生这种情况。
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <mpi.h>
#include "mkl_lapack.h"
int main(int argc, char **argv)
{
int N,dim,info;
double *d,*e,*z,*work;
char jobz='V';
switch(argc)
{
case 2:
N=atoi(argv[1]);
break;
default:
fprintf(stderr,"Errore nell' inserimento degli argomenti\n");
exit(EXIT_FAILURE);
break;
}
if(N%2!=0)
{
fprintf(stderr,"La dimensione della matrice deve essere pari\n");
exit(EXIT_FAILURE);
}
dim=N/2;
d=(double*)malloc(N*sizeof(double));
e=(double*)malloc((N-1)*sizeof(double));
z=(double*)malloc(N*N/2*sizeof(double));
work=(double*)malloc((N-1)*2*sizeof(double));
for(int i=0;i<N-1;i++)
{
d[i]=(double)(i+3);
e[i]=1.0;
}
dstev(&jobz,&dim,d,e,z,&dim,work,&info);
dim--;
dstev(&jobz,&dim,&d[N/2],&e[N/2],&z[N*N/4],&dim,&work[N-1],&info);
for(int i=0;i<(N*N/2);i++)
printf("(%f) ",z[i]);
return 0;
}
我希望我清楚地解释了这件事,让我知道有些事情不清楚。
最佳答案
这些dstev()
的调用是正确的,因为info
之后为0。dstev()
的两个调用之间有区别:dim
通过执行dim--
递减。
传递给dstev()
的第一个矩阵的大小为N/2
,第二个矩阵的大小为N/2-1
。打印z
元素时,N*N/4+(N-2)*(N-2)/4=N*N/2-N+1
的总大小为N*N/2
。
因此,N-1
的最后z
个元素是没有意义的。在这种情况下,发现它们为零。
重新整理dim--
可解决问题:z
中不再有零,除非您更改d
或e
。
用gcc main.c -o main -llapack -lblas -lm
编译的代码:
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
//#include <mpi.h>
extern dstev_(char* jobz,int* n,double* d,double* e,double* z,int* ldz,double* work,int* info);
int main(int argc, char **argv)
{
int N,dim,info;
double *d,*e,*z,*work;
char jobz='V';
switch(argc)
{
case 2:
N=atoi(argv[1]);
break;
default:
fprintf(stderr,"Errore nell' inserimento degli argomenti\n");
exit(EXIT_FAILURE);
break;
}
if(N%2!=0)
{
fprintf(stderr,"La dimensione della matrice deve essere pari\n");
exit(EXIT_FAILURE);
}
dim=N/2;
d=malloc(N*sizeof(double));
e=malloc((N-1)*sizeof(double));
z=malloc(N*N/2*sizeof(double));
work=malloc((N-1)*2*sizeof(double));
int i;
for(i=0;i<N-1;i++)
{
d[i]=(double)(i+3);
e[i]=1.0;
}
dstev_(&jobz,&dim,d,e,z,&dim,work,&info);
printf("info %d\n",info);
//dim--;
dstev_(&jobz,&dim,&d[N/2],&e[N/2],&z[N*N/4],&dim,&work[N-1],&info);
printf("info %d\n",info);
for( i=0;i<(N*N/2);i++)
printf("(%e) ",z[i]);
free(e);
free(d);
free(z);
free(work);
return 0;
}
关于c - 我正在使用dstev计算特征向量得到零值,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/8840125/