2016-06-26 17 views
0

複雑なエルミート行列を対角化するためにlapackのzheevdを使用しようとしています。私は、コードに結果を実行するときlapack zheevdが間違った結果を返しました

program test 
    implicit none 

    INTEGER, PARAMETER :: N=4 
    INTEGER, PARAMETER :: LDA = N 
    INTEGER, PARAMETER :: LWMAX = 1000 
    INTEGER :: INFO, LWORK, LIWORK, LRWORK,i,j 

    INTEGER ::  IWORK(LWMAX) 
    REAL(8) :: W(N), RWORK(LWMAX) 
    COMPLEX(16) :: A(LDA, N), WORK(LWMAX), zero 
    character(len=1) :: job,uplo 



! the matrix I want to diagonalize is: 
!  ( 3.40, 0.00) (-2.36, -1.93) (-4.68, 9.55) ( 5.37, -1.23) 
! A= (-2.36, 1.93) ( 6.94, 0.00) ( 8.13, -1.47) ( 2.07, -5.78) 
!  (-4.68, -9.55) ( 8.13, 1.47) (-2.14, 0.00) ( 4.68, 7.44) 
!  ( 5.37, 1.23) ( 2.07, 5.78) ( 4.68, -7.44) (-7.42, 0.00) 

    zero=dcmplx(0.0d0,0.0d0) 

    A=zero 
    A(1,1)= dcmplx(3.40d0, 0.0d0); A(1,2)=dcmplx(-2.36d0, -1.93d0); A(1,3)= dcmplx(-4.68d0,9.55d0) 
    A(1,4)= dcmplx(5.37d0, -1.23d0) 
    A(2,2)= dcmplx(6.94d0, 0.0d0); A(2,3)=dcmplx(8.13d0, -1.47d0); A(2,4)= dcmplx(2.07d0, -5.78d0) 
    A(3,3)= dcmplx(-2.14d0, 0.0d0); A(3,4)=dcmplx(4.68d0, 7.44d0); A(4,4)= dcmplx(-7.42d0, 0.0d0) 


    job='V'; uplo='U' 

    LWORK= N**2 + 2*N; LRWORK= 2*N**2 + 5*N + 1; LIWORK= 5*N+3 

    CALL ZHEEVD(job, uplo, N, A, LDA, W, WORK, LWORK, RWORK,LRWORK,IWORK,LIWORK, INFO) 

    IF(INFO > 0) THEN 
    WRITE(*,*)'The algorithm failed to compute eigenvalues.' 
    STOP 
    END IF 


    print*, 'eigenvalues found' 
    do i=1,N 
    print*, W(i) 
    end do 

    open(1, file='eigenvectors.dat') 

    write(1,10) ((A(i,j),j=1,N),i=1,N) 
10 format(4(F10.5,2X,F10.5))  


    end program test 

は私が取得:私はここでは、コードの...任意のコンパイルや実行時エラーを生成しますが、固有値のために間違った結果を示していない小さな例を書かれまし」固有値は以下のとおりです。 -2.8413、0、0、2.8413

実際の固有値がありながら:-21.968、16.3387、6.45946、-0.0501069

私は、ルーチンのリファレンスガイドを見ておくと、私が正しいすべてを持っているようですそれは正常に動作する必要があります期待していない...私のコードで何が間違っているのか誰かが考えている?

おかげ

+0

'COMPLEX(16)'は 'COMPLEX * 16'と同じではありません。 – talonmies

+0

これを解決できましたか?私が提供した答えは助けになりましたか? – talonmies

答えて

0

私が見ることができることをここに三つの主要な問題があります。

  1. 最も深刻な問題は、あなたがCOMPLEX(16)ようにあなたのコードをベースとしているMKLの例でCOMPLEX*16種類を翻訳しているということですが。それは間違っている。 COMPLEX(8)を使用してください。あなたのツールチェーンが実際に拡張された精度の複合型を持っているかどうかはわかりませんが、あなたのコードとLAPACK呼び出しの間にサイズの不一致があるかもしれません。
  2. コードには、渡す行列の値LAPACKはあなたのコメントと同じではない(おそらく、あなたが固有値を計算した行列と同じではない)
  3. 最後に、重要なことに、ZHEEVDのインターフェイスを定義していない)。これにより、暗黙のインターフェイスがコンパイラによって推測されるようになります。コード内で引き渡される引数とLAPACKが期待するものとの間に矛盾がある可能性は非常に高いです。特に複合引数の型不一致がある。

私は3つの組み合わせによって結果を修正する必要があると思います。

関連する問題