Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse code

Implement eigh() and a test

  • Loading branch information...
commit 9a0b34e22ffee3822da7747ce62efbcaf3e7d8a4 1 parent 67df0d1
Ondřej Čertík authored
2  src/CMakeLists.txt
@@ -5,7 +5,7 @@ set(SRC
5 5 )
6 6
7 7 if(WITH_LAPACK)
8   - set(SRC ${SRC} lapack.f90 splines.f90)
  8 + set(SRC ${SRC} lapack.f90 splines.f90 linalg.f90)
9 9 endif()
10 10
11 11 if(WITH_HDF5)
50 src/linalg.f90
... ... @@ -0,0 +1,50 @@
  1 +module linalg
  2 +use types, only: dp
  3 +use lapack, only: dsygvd
  4 +use utils, only: stop_error
  5 +implicit none
  6 +private
  7 +public eigh
  8 +
  9 +contains
  10 +
  11 +subroutine eigh(Am, Bm, lam, c)
  12 +! solves generalized eigen value problem for all eigenvalues and eigenvectors
  13 +! Am must by symmetric, Bm symmetric positive definite.
  14 +! Only the lower triangular part of Am and Bm is used.
  15 +real(dp), intent(in) :: Am(:,:) ! LHS matrix: Am c = lam Bm c
  16 +real(dp), intent(in) :: Bm(:,:) ! RHS matrix: Am c = lam Bm c
  17 +real(dp), intent(out) :: lam(:) ! eigenvalues: Am c = lam Bm c
  18 +real(dp), intent(out) :: c(:,:) ! eigenvectors: Am c = lam Bm c; c(i,j) = ith component of jth vec.
  19 +integer n
  20 +! lapack variables
  21 +integer lwork, liwork, info
  22 +integer, allocatable:: iwork(:)
  23 +real(dp), allocatable:: Amt(:,:), Bmt(:,:), work(:)
  24 +
  25 +! solve
  26 +n = size(Am,1)
  27 +lwork = 1 + 6*n + 2*n**2
  28 +liwork = 3 + 5*n
  29 +allocate(Amt(n,n), Bmt(n,n), work(lwork), iwork(liwork))
  30 +Amt = Am; Bmt = Bm ! Amt,Bmt temporaries overwritten by dsygvd
  31 +call dsygvd(1,'V','L',n,Amt,n,Bmt,n,lam,work,lwork,iwork,liwork,info)
  32 +if (info /= 0) then
  33 + print *, "dsygvd returned info =", info
  34 + if (info < 0) then
  35 + print *, "the", -info, "-th argument had an illegal value"
  36 + else if (info <= n) then
  37 + print *, "the algorithm failed to compute an eigenvalue while working"
  38 + print *, "on the submatrix lying in rows and columns", 1.0_dp*info/(n+1)
  39 + print *, "through", mod(info, n+1)
  40 + else
  41 + print *, "The leading minor of order ", info-n, &
  42 + "of B is not positive definite. The factorization of B could ", &
  43 + "not be completed and no eigenvalues or eigenvectors were computed."
  44 + end if
  45 + call stop_error('eigh: dsygvd error')
  46 +end if
  47 +c = Amt
  48 +end subroutine
  49 +
  50 +end module
1  tests/CMakeLists.txt
@@ -6,6 +6,7 @@ add_subdirectory(strings)
6 6 add_subdirectory(mesh)
7 7 if(WITH_LAPACK)
8 8 add_subdirectory(splines)
  9 + add_subdirectory(linalg)
9 10 endif()
10 11 if(WITH_HDF5)
11 12 add_subdirectory(hdf5)
12 tests/linalg/CMakeLists.txt
... ... @@ -0,0 +1,12 @@
  1 +include_directories(${PROJECT_BINARY_DIR}/src)
  2 +
  3 +project(linalg)
  4 +
  5 +add_executable(test_eig test_eig.f90)
  6 +set(LAPACK_LIBS
  7 + lapack
  8 + blas
  9 + )
  10 +target_link_libraries(test_eig fortran_utils ${LAPACK_LIBS})
  11 +
  12 +add_test(test_eig ${PROJECT_BINARY_DIR}/test_eig)
38 tests/linalg/test_eig.f90
... ... @@ -0,0 +1,38 @@
  1 +program test_mesh
  2 +use types, only: dp
  3 +use utils, only: assert
  4 +use linalg, only: eigh
  5 +implicit none
  6 +
  7 +real(dp), parameter :: eps = 1e-9_dp
  8 +real(dp) :: A(2, 2), B(2, 2), lam(2), c(2, 2), r(2), n
  9 +integer :: i
  10 +A = reshape([1, 0, 0, -1], [2, 2])
  11 +B = reshape([1, 0, 0, 1], [2, 2])
  12 +call eigh(A, B, lam, c)
  13 +! Test eigenvalues:
  14 +call assert(all(abs(lam - [-1, 1]) < eps))
  15 +do i = 1, 2
  16 + ! Test that c(:, i) are eigenvectors:
  17 + r = matmul(A-lam(i)*B, c(:, i))
  18 + call assert(sqrt(dot_product(r, r)) < eps)
  19 + ! Test that eigenvectors are properly normalized:
  20 + n = dot_product(c(:, i), matmul(B, c(:, i)))
  21 + call assert(abs(n - 1) < eps)
  22 +end do
  23 +
  24 +A = reshape([2, -4, -4, 2], [2, 2])
  25 +B = reshape([2, 1, 1, 2], [2, 2])
  26 +call eigh(A, B, lam, c)
  27 +! Test eigenvalues:
  28 +call assert(all(abs(lam - [-2._dp/3, 6._dp]) < eps))
  29 +do i = 1, 2
  30 + ! Test that c(:, i) are eigenvectors:
  31 + r = matmul(A-lam(i)*B, c(:, i))
  32 + call assert(sqrt(dot_product(r, r)) < eps)
  33 + ! Test that eigenvectors are properly normalized:
  34 + n = dot_product(c(:, i), matmul(B, c(:, i)))
  35 + call assert(abs(n - 1) < eps)
  36 +end do
  37 +
  38 +end program

0 comments on commit 9a0b34e

Please sign in to comment.
Something went wrong with that request. Please try again.