I'm getting zero values using dstev to calculate eigenvectors

Question

I use gcc to compile under mac os x, I have Intel's mkl_lapack.h library installed. In the program I have a NxN tridiagonal matrix, so I just use two vectors to store values of the matrix. "d" vector is the main diagonal, the values of subdiagonals are stored in "e". First of all I initialize values, then since the matrix is 16x16 (in input I'm giving the value 16 as argv[1]), I split the vector into two vectors (I could just use dstev once for all, but it's for experimental purposes), from d[0] to d[N/2-1] I have the first vector, from d[N/2] to d[N-1] the second one. So once initilized the values of "e" and "d" , I call two times dstev. But I don't bother writing all the values in "z" (z will contain eigenvectors), because I know that after calling dstev two times, in all the "z" vector I should have only two submatrixes of values, 8x8 of non-zero values. But if I try priting "z", some values are 0.0, and I can't explain why this happens.

#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;
}

I hope I explained this thing clearly, let me know is something isn't clear.

Try asking this at [SciComp](http://www.scicomp.stackexchange.com) — , Mar 08 '12 at 13:45
The documentation tells you that the value returned in `info` contains diagnostic information. You appear to be ignoring that. — David Heffernan, Jun 23 '14 at 10:05

francis · Answer 1 · 2015-02-20T17:31:29.213

These calls of dstev() are correct since info is 0 after.

There is a difference between the two calls of dstev() : dim is decremented by doing dim--.

The first matrix passed to dstev() is of size N/2 and the second matrix is of size N/2-1. The total size of z is N*N/4+(N-2)*(N-2)/4=N*N/2-N+1 while N*N/2 elements are printed. Hence, the last N-1 elements of z are meaningless. In this case, they are found to be zeros.

Revoming dim-- solves the problem : there are no more zeros in z, except if you change d or e.

Code compiled with 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;
}

To help distinguish between small numbers and zero recommend: `printf("(%e) ",z[i]);` e vs. f. — chux - Reinstate Monica, Feb 20 '15 at 17:25
@chux : right ! Scientific notations are much more practical for large or small numbers. I edited accordingly. — francis, Feb 20 '15 at 17:33

I'm getting zero values using dstev to calculate eigenvectors

1 Answers1