g95で実対称行列の固有値を求めるサブルーチンには,

  • dsyev (driver)
  • dsyevd (divide and conquer)
  • dsyevx (expert)
  • dsyevr (relatively robust representations)

といくつかある(LAPACKのユーザガイド及びmanページ参照).以下の簡単な例では,最初の3つの固有値は一致した.dsyevrはe-14〜e-15程度の差があった.

Macには,CBLAS, CLAPACKがAccelerate FrameworkとしてOSに添付されている.g95がインストールしてあれば,

g95 foo.f90 -framework accelerate

でライブラリをリンクできる.とっても簡単.

program eigen
  implicit none

  character, parameter :: jobz = "V", uplo = "L"
  integer, parameter :: n = 3, lda = n, lwork = 3*n-1
  integer :: info, i
  real*8, dimension(lda,n) :: a
  real*8, dimension(n) :: w
  real*8, dimension(lwork) :: work

  a(:,1) = (/1.0d0,  2.0d0,  3.0d0/)
  a(:,2) = (/2.0d0,  3.0d0, -4.0d0/)
  a(:,3) = (/3.0d0, -4.0d0,  5.0d0/)

  call dsyev(jobz, uplo, n, a, lda, w, work, lwork, info)
  print *, "info=", info
  do i=1, n
    print *, "mode:", i, " eigenvalue=", w(i)
    print *, "         eigenvector", a(:,i)
  end do
end program eigen