Я получаю нулевые значения, используя dstev для вычисления собственных векторов - PullRequest
Новая
       9

Я получаю нулевые значения, используя dstev для вычисления собственных векторов

3 голосов
/ 12 января 2012

Я использую 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] второй. Поэтому, как только инициализируются значения «е» и «д», я два раза называю дстев. Но я не пишу все значения в «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;
}

Я надеюсь, что объяснил это ясно, дайте мне знать, что-то не ясно.

1 Ответ

1 голос
/ 20 февраля 2015

Эти вызовы 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 не имеют смысла. В этом случае они оказываются нулями.

Revoming 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;
}
...