2017-01-24 5 views
0

Fortran Eispackルーチン(固有値と固有ベクトルを計算し、値がMatlabのものと異なるかどうかをチェックする)を使うコードを書いていましたが、プログラムがハングアップするqzhesサブルーチンを呼び出します。サブルーチンを呼び出してプログラムをクラッシュし、行列を渡します

私はファイルから行列を読み込みます。

呼び出しにコメントしましたが、問題なく動作します。

私はFortranのことを学んだし、インターネットの助けを借りて、私はこのコード(コンパイルおよび実行)を書いた:

program qz 
    IMPLICIT NONE 
INTEGER:: divm, i, divg 
INTEGER(kind=4) :: dimen 
LOGICAL :: matz 
REAL(kind = 8), DIMENSION(:,:), ALLOCATABLE:: ma 
REAL(kind = 8), DIMENSION(:), ALLOCATABLE:: tabm 
REAL(kind = 8), DIMENSION(:,:), ALLOCATABLE:: ga 
REAL(kind = 8), DIMENSION(:), ALLOCATABLE:: tabg 
REAL(kind = 8), DIMENSION(:,:), ALLOCATABLE:: zet 

divm = 1 
divg = 2 
dimen = 20 
matz = .TRUE. 
ALLOCATE(ma(1:dimen,1:dimen)) 
ALLOCATE(tabm(1:dimen)) 
ALLOCATE(ga(1:dimen,1:dimen)) 
ALLOCATE(tabg(1:dimen)) 

OPEN(divm, FILE='Em.txt') 

DO i=1,dimen 
    READ (divm,*) tabm 
    ma(1:dimen,i)=tabm 
END DO 

CLOSE(divm) 

OPEN(divg, FILE='Gje.txt') 

DO i=1,dimen 
    READ (divg,*) tabg 
    ga(1:dimen,i)=tabg 
END DO 

CLOSE(divg) 

call qzhes(dimen, ma, ga, matz, zet) 

OPEN(divm, FILE='Em2.txt') 

DO i=1,dimen 
    tabm = ma(1:dimen,i) 
    WRITE (divm,*) tabm 
END DO 

CLOSE(divm) 

OPEN(divg, FILE='Gje2.txt') 

DO i=1,dimen 
    tabg = ga(1:dimen,i) 
    WRITE (divg,*) tabg 
END DO 

CLOSE(divg) 

end program qz 

...//EISPACK subrotines//... 

マトリックス:

Gje.txt:https://drive.google.com/file/d/0BxH3QOkswLy_c2hmTGpGVUI3NzQ/view?usp=sharing

エム.txt:https://drive.google.com/file/d/0BxH3QOkswLy_OEtJUGQwN3ZXX2M/view?usp=sharing

編集:

subroutine qzhes (n, a, b, matz, z) 

!*****************************************************************************80 
! 
!! QZHES carries out transformations for a generalized eigenvalue problem. 
! 
! Discussion: 
! 
! This subroutine is the first step of the QZ algorithm 
! for solving generalized matrix eigenvalue problems. 
! 
! This subroutine accepts a pair of real general matrices and 
! reduces one of them to upper Hessenberg form and the other 
! to upper triangular form using orthogonal transformations. 
! it is usually followed by QZIT, QZVAL and, possibly, QZVEC. 
! 
! Licensing: 
! 
! This code is distributed under the GNU LGPL license. 
! 
! Modified: 
! 
! 18 October 2009 
! 
! Author: 
! 
! Original FORTRAN77 version by Smith, Boyle, Dongarra, Garbow, Ikebe, 
! Klema, Moler. 
! FORTRAN90 version by John Burkardt. 
! 
! Reference: 
! 
! James Wilkinson, Christian Reinsch, 
! Handbook for Automatic Computation, 
! Volume II, Linear Algebra, Part 2, 
! Springer, 1971, 
! ISBN: 0387054146, 
! LC: QA251.W67. 
! 
! Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, 
! Yasuhiko Ikebe, Virginia Klema, Cleve Moler, 
! Matrix Eigensystem Routines, EISPACK Guide, 
! Lecture Notes in Computer Science, Volume 6, 
! Springer Verlag, 1976, 
! ISBN13: 978-3540075462, 
! LC: QA193.M37. 
! 
! Parameters: 
! 
! Input, integer (kind = 4) N, the order of the matrices. 
! 
! Input/output, real (kind = 8) A(N,N). On input, the first real general 
! matrix. On output, A has been reduced to upper Hessenberg form. The 
! elements below the first subdiagonal have been set to zero. 
! 
! Input/output, real (kind = 8) B(N,N). On input, a real general matrix. 
! On output, B has been reduced to upper triangular form. The elements 
! below the main diagonal have been set to zero. 
! 
! Input, logical MATZ, should be TRUE if the right hand transformations 
! are to be accumulated for later use in computing eigenvectors. 
! 
! Output, real (kind = 8) Z(N,N), contains the product of the right hand 
! transformations if MATZ is TRUE. 
! 
    implicit none 

    integer (kind = 4) n 

    real (kind = 8) a(n,n) 
    real (kind = 8) b(n,n) 
    integer (kind = 4) i 
    integer (kind = 4) j 
    integer (kind = 4) k 
    integer (kind = 4) l 
    integer (kind = 4) l1 
    integer (kind = 4) lb 
    logical matz 
    integer (kind = 4) nk1 
    integer (kind = 4) nm1 
    real (kind = 8) r 
    real (kind = 8) rho 
    real (kind = 8) s 
    real (kind = 8) t 
    real (kind = 8) u1 
    real (kind = 8) u2 
    real (kind = 8) v1 
    real (kind = 8) v2 
    real (kind = 8) z(n,n) 
! 
! Set Z to the identity matrix. 
! 
    if (matz) then 

    z(1:n,1:n) = 0.0D+00 

    do i = 1, n 
     z(i,i) = 1.0D+00 
    end do 

    end if 
! 
! Reduce B to upper triangular form. 
! 
    if (n <= 1) then 
    return 
    end if 

    nm1 = n - 1 

    do l = 1, n - 1 

    l1 = l + 1 

    s = sum (abs (b(l+1:n,l))) 

    if (s /= 0.0D+00) then 

     s = s + abs (b(l,l)) 
     b(l:n,l) = b(l:n,l)/s 

     r = sqrt (sum (b(l:n,l)**2)) 
     r = sign (r, b(l,l)) 
     b(l,l) = b(l,l) + r 
     rho = r * b(l,l) 

     do j = l + 1, n 

     t = dot_product (b(l:n,l), b(l:n,j)) 

     b(l:n,j) = b(l:n,j) - t * b(l:n,l)/rho 

     end do 

     do j = 1, n 

     t = dot_product (b(l:n,l), a(l:n,j)) 

     a(l:n,j) = a(l:n,j) - t * b(l:n,l)/rho 

     end do 

     b(l,l) = - s * r 
     b(l+1:n,l) = 0.0D+00 

    end if 

    end do 
! 
! Reduce A to upper Hessenberg form, while keeping B triangular. 
! 
    if (n == 2) then 
    return 
    end if 

    do k = 1, n - 2 

    nk1 = nm1 - k 

    do lb = 1, nk1 

     l = n - lb 
     l1 = l + 1 
! 
! Zero A(l+1,k). 
! 
     s = abs (a(l,k)) + abs (a(l1,k)) 

     if (s /= 0.0D+00) then 

     u1 = a(l,k)/s 
     u2 = a(l1,k)/s 
     r = sign (sqrt (u1**2 + u2**2), u1) 
     v1 = - (u1 + r)/r 
     v2 = - u2/r 
     u2 = v2/v1 

     do j = k, n 
      t = a(l,j) + u2 * a(l1,j) 
      a(l,j) = a(l,j) + t * v1 
      a(l1,j) = a(l1,j) + t * v2 
     end do 

     a(l1,k) = 0.0D+00 

     do j = l, n 
      t = b(l,j) + u2 * b(l1,j) 
      b(l,j) = b(l,j) + t * v1 
      b(l1,j) = b(l1,j) + t * v2 
     end do 
! 
! Zero B(l+1,l). 
! 
     s = abs (b(l1,l1)) + abs (b(l1,l)) 

     if (s /= 0.0) then 

      u1 = b(l1,l1)/s 
      u2 = b(l1,l)/s 
      r = sign (sqrt (u1**2 + u2**2), u1) 
      v1 = -(u1 + r)/r 
      v2 = -u2/r 
      u2 = v2/v1 

      do i = 1, l1 
      t = b(i,l1) + u2 * b(i,l) 
      b(i,l1) = b(i,l1) + t * v1 
      b(i,l) = b(i,l) + t * v2 
      end do 

      b(l1,l) = 0.0D+00 

      do i = 1, n 
      t = a(i,l1) + u2 * a(i,l) 
      a(i,l1) = a(i,l1) + t * v1 
      a(i,l) = a(i,l) + t * v2 
      end do 

      if (matz) then 

      do i = 1, n 
       t = z(i,l1) + u2 * z(i,l) 
       z(i,l1) = z(i,l1) + t * v1 
       z(i,l) = z(i,l) + t * v2 
      end do 

      end if 

     end if 

     end if 

    end do 

    end do 

    return 
end 
+0

私のミス、と改めて確認することができます。サブルーチンコードを編集して追加します。 が動作を停止するようにプロンプ​​トを表示すると、プログラムがクラッシュします。 –

+0

それは唯一の出力ですか、それともその行の前に何かを印刷しましたか? –

+0

Googleドライブへのリンクを削除し、ここで入力データをコピーします。この質問は、削除可能なファイルストレージに依存してはいけません。将来の訪問者には意味をなさない必要があります。 –

答えて

0

あなたは多分今はmemoryleakを持っていない、ので、私は、割り当てプロセス

integer :: status1, status2, status3, status4, status5 
! check the allocation, returnvalue 0 means ok 

ALLOCATE(ma(1:dimen,1:dimen), stat=status1) 
ALLOCATE(tabm(1:dimen), stat=status2) 
ALLOCATE(ga(1:dimen,1:dimen), stat=status3) 
ALLOCATE(tabg(1:dimen), stat=status4) 
ALLOCATE(zet(1:dimen,1:dimen), stat=status5) 

とプログラムDEALLOCATEの終わりにすべてのアレイを展開するだろうが、あなたはサブルーチンにこのプログラムを置くと、いくつかのそれを使用する場合プログラム実行中に大規模なmatriciesで時間、プログラムは深刻な記憶を漏らす可能性があります。

.... 

DO i=1,dimen 
    tabg = ga(1:dimen,i) 
    WRITE (divg,*) tabg 
END DO 

CLOSE(divg) 

DEALLOCATE(ma, stat=status1) 
DEALLOCATE(tabm, stat=status2) 
DEALLOCATE(ga, stat=status3) 
DEALLOCATE(tabg, stat=status4) 
DEALLOCATE(zet, stat=status5) 

あなたは「QZ」私は意味「QZHES」によるステータス整数、割り当て解除はOKだった場合は、戻り値を再び0

+1

サブルーチンまたは関数に割り当てられた変数は、サブルーチンの最後に自動的に割り当て解除されます(それらに 'move_alloc'を使用しない限り)。 –

関連する問題