2012-01-12 3 views
3

を計算するdstevを使用してゼロ値を取得しています、私はインテルのmkl_lapack.hライブラリがインストールされています。 プログラムにはNxNの三重対角行列がありますので、行列の値を格納するために2つのベクトルを使用します。 "d"ベクトルが主対角線である場合、副対角線の値は "e"に格納されます。私はマトリックスが16×16であるため、その後、私は2つのベクトル(私はすべてのために一度dstev使用することができるベクターを分割し、それはです(入力にI [1] argvのように値16を与えている)、値を初期化し、全ての ファーストd [0]からd [N/2-1]には、d [N/2]からd [N-1]までの最初のベクトルがあります。 "e"と "d"の値を初期化すると、私は2回dstevを呼び出します。 私はすべての「Z」ベクトルで2回、dstev呼び出した後、私は、8×8の値の2つのだけ部分行列を持つべきであることを知っているので、しかし、私は、(zは固有ベクトルが含まれています)、「Z」のすべての値を書く気にしないでくださいゼロ以外の値 しかし、私が "z"を整理しようとすると、いくつかの値は0.0であり、なぜこれが起こるのか説明できません。は、私はマックOS Xの下コンパイルにはgccを使用固有ベクトル

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

私はこのことを明確に説明したいと思っていますが、何かが分かりません。 infoが0の後にあるので、

+0

[SciComp](http://www.scicomp.stackexchange.com)でこれを尋ねてみ –

+0

ドキュメントはinfo' 'で返された値が診断含まれていることを示しています情報。あなたはそれを無視しているようです。 –

答えて

1

dstev()のこれらの呼び出しは正しいです。

dstev()の2つの呼び出しの間に違いがあります。は、dim--を実行することによって減分されます。

dstev()に渡される最初の行列のサイズはN/2であり、2番目の行列のサイズはN/2-1です。 N*N/2要素が印刷されている間zの合計サイズがN*N/4+(N-2)*(N-2)/4=N*N/2-N+1あります。 はしたがって、zの最後N-1要素は無意味です。この場合、それらはゼロであることがわかります。

Revoming dim--は、問題を解決:これ以上のゼロは、あなたがdeを変更した場合を除き、zではありません。 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; 
} 
+1

小さい数字とゼロを区別するために、 'printf("(%e) "、z [i]);' eとfを比較してください。 – chux

+0

@chux:右!科学的表記法は、大規模または小規模の方がはるかに実用的です。私はそれに応じて編集した。 – francis

関連する問題