我使用 dstev 计算特征向量得到零值

Question

我使用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;
}

我希望我能清楚地解释这件事，让我知道有什么不清楚的地方。

Answer 1

这些 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 个元素。
因此，z 的最后一个 N-1 元素是没有意义的。在这种情况下，发现它们为零。

删除 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;
}