From ec51d8afb634bcd8ad8f27ef09c5e57abb0f6e9e Mon Sep 17 00:00:00 2001 From: Irukoa Date: Thu, 12 Sep 2024 15:57:29 +0200 Subject: [PATCH] Include matrix inversion utility. --- README.md | 10 +++ app/main.F90 | 2 +- fpm.toml | 2 +- src/WannInt.F90 | 3 +- src/WannInt_utilities.F90 | 83 +++++++++++++++++++++++ test/suites/Basic_Suite.F90 | 97 ++++++++++++++++++++++++++- test/suites/Extra_Utilities_Suite.F90 | 1 - 7 files changed, 193 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index 53bc804..6c18b8b 100644 --- a/README.md +++ b/README.md @@ -296,6 +296,16 @@ where - `real/complex(dp), optional, intent(out) :: D(n, n)` is the diagonal form $D$ of `matrix`. - `real(dp), optional, intent(out) :: eig(n)` are the eigenvalues of `matrix`. +## Matrix inverse utility + +The library includes a utility to invert matrices by employing singular value decomposition. It can be called by +```fortran +inv = inverse(matrix) +``` +where +- `real/complex(dp), intent(in) :: matrix(n, n)` is an invertible square matrix $A$. +- `real/complex(dp) :: inv(n, n)` is an invertible square matrix $B$ obeying $AB=I$. + ## Schur decomposition utility The library includes a wrapper of `zgees` routine from LAPACK. It can be called by diff --git a/app/main.F90 b/app/main.F90 index 1a5f0a5..630dc50 100644 --- a/app/main.F90 +++ b/app/main.F90 @@ -6,7 +6,7 @@ program main implicit none - character(len=10) :: ver = "1.0.0" + character(len=10) :: ver = "1.0.1" write (unit=output_unit, fmt="(A)") "" write (unit=output_unit, fmt="(A)") " ____ __ ____ ___ .__ __. .__ __. __ .__ __. .___________. " diff --git a/fpm.toml b/fpm.toml index 5f9a71c..e9e6794 100644 --- a/fpm.toml +++ b/fpm.toml @@ -1,5 +1,5 @@ name = "WannInt" -version = "1.0.0" +version = "1.0.1" license = "GNU General Public License v3.0" author = "Álvaro R. Puente-Uriona" maintainer = "alvaro.ruiz@ehu.eus" diff --git a/src/WannInt.F90 b/src/WannInt.F90 index cc1ec61..a5ae3c6 100644 --- a/src/WannInt.F90 +++ b/src/WannInt.F90 @@ -5,7 +5,7 @@ module WannInt use WannInt_kinds, only: wp => dp use WannInt_definitions, only: cmplx_0, cmplx_i, pi - use WannInt_utilities, only: diagonalize, dirac_delta, & + use WannInt_utilities, only: diagonalize, inverse, dirac_delta, & deg_list, schur, & SVD, expsh, logu use MAC, only: container_specifier, container @@ -15,6 +15,7 @@ module WannInt private public :: diagonalize + public :: inverse public :: dirac_delta public :: deg_list public :: schur diff --git a/src/WannInt_utilities.F90 b/src/WannInt_utilities.F90 index acdeea0..1f68535 100644 --- a/src/WannInt_utilities.F90 +++ b/src/WannInt_utilities.F90 @@ -9,12 +9,18 @@ module WannInt_utilities private public :: diagonalize + public :: inverse interface diagonalize module procedure :: diagonalize_symmetric module procedure :: diagonalize_hermitian end interface + interface inverse + module procedure :: inverse_real + module procedure :: inverse_complex + end interface + public :: dirac_delta public :: deg_list public :: schur @@ -395,4 +401,81 @@ function logu(matrix) end function logu + function inverse_real(matrix) result(inv_r) + !=========================================! + ! ! + ! Given a general real matrix, ! + ! try to compute its inverse. ! + ! ! + !=========================================! + + real(wp), intent(in) :: matrix(:, :) + real(wp) :: inv_r(size(matrix(:, 1)), size(matrix(1, :))) + + complex(wp) :: matrix_c(size(matrix(:, 1)), size(matrix(1, :))), & + U(size(matrix(:, 1)), size(matrix(1, :))), & + V(size(matrix(:, 1)), size(matrix(1, :))), & + S(size(matrix(:, 1)), size(matrix(1, :))) + + integer :: i + + character(len=1024) :: errormsg + + if (size(matrix(:, 1)) /= size(matrix(1, :))) error stop & + "WannInt: Error #1: matrix to invert is not square." + + matrix_c = cmplx(matrix, 0.0_wp, wp) + + call SVD(matrix=matrix_c, U=U, V=V, sigma=S) + + do i = 1, size(matrix(:, 1)) + if (abs(real(S(i, i), wp)) < 1.0E-8_wp) then + write (errormsg, "(i20)") i + errormsg = "WannInt: Error #2: the eigenvalue #"//trim(adjustl(errormsg))//" is 0, cannot invert." + error stop trim(errormsg) + endif + S(i, i) = cmplx(1.0_wp/real(S(i, i), wp), 0.0_wp, wp) + enddo + + inv_r = real(matmul(matmul(V, S), conjg(transpose(U))), wp) + + end function inverse_real + + function inverse_complex(matrix) result(inv_c) + !=========================================! + ! ! + ! Given a general real matrix, ! + ! try to compute its inverse. ! + ! ! + !=========================================! + + complex(wp), intent(in) :: matrix(:, :) + complex(wp) :: inv_c(size(matrix(:, 1)), size(matrix(1, :))) + + complex(wp) :: U(size(matrix(:, 1)), size(matrix(1, :))), & + V(size(matrix(:, 1)), size(matrix(1, :))), & + S(size(matrix(:, 1)), size(matrix(1, :))) + + integer :: i + + character(len=1024) :: errormsg + + if (size(matrix(:, 1)) /= size(matrix(1, :))) error stop & + "WannInt: Error #1: matrix to invert is not square." + + call SVD(matrix=matrix, U=U, V=V, sigma=S) + + do i = 1, size(matrix(:, 1)) + if (abs(real(S(i, i), wp)) < 1.0E-8_wp) then + write (errormsg, "(i20)") i + errormsg = "WannInt: Error #2: the eigenvalue #"//trim(adjustl(errormsg))//" is 0, cannot invert." + error stop trim(errormsg) + endif + S(i, i) = cmplx(1.0_wp/real(S(i, i), wp), 0.0_wp, wp) + enddo + + inv_c = matmul(matmul(V, S), conjg(transpose(U))) + + end function inverse_complex + end module WannInt_utilities diff --git a/test/suites/Basic_Suite.F90 b/test/suites/Basic_Suite.F90 index 517dbd3..b972196 100644 --- a/test/suites/Basic_Suite.F90 +++ b/test/suites/Basic_Suite.F90 @@ -1,7 +1,7 @@ module Basic_Suite use WannInt_kinds, only: wp => dp use WannInt_definitions, only: cmplx_0, cmplx_1 - use WannInt_utilities, only: diagonalize + use WannInt_utilities, only: diagonalize, inverse use WannInt, only: crystal use testdrive, only: new_unittest, unittest_type, error_type implicit none @@ -20,6 +20,8 @@ subroutine Collect_Basic_Tests(testsuite) new_unittest("Diagonalization of Hermitian matrices (2x2)", test_diagonalization_herm_22), & new_unittest("Diagonalization of symmetric matrices (random)", test_diagonalization_sym_rand), & new_unittest("Diagonalization of Hermitian matrices (random)", test_diagonalization_herm_rand), & + new_unittest("Inverse of real matrices (random)", test_inv_re_rand), & + new_unittest("Inverse of complex matrices (random)", test_inv_c_rand), & new_unittest("Initialization of system and return basic properties from input.", create_crystal_from_input), & new_unittest("Initialization of system and return basic properties from file.", create_crystal_from_file)] @@ -146,6 +148,99 @@ subroutine test_diagonalization_sym_rand(error) end subroutine test_diagonalization_sym_rand + subroutine test_inv_re_rand(error) + type(error_type), allocatable, intent(out) :: error + + real(wp), allocatable :: mat(:, :), prod(:, :) + real(wp) :: szr + real(wp) :: offdiag + integer :: sz, n, m + + call random_seed() + call random_number(szr) + + sz = nint(50.0_wp*szr) + 5 + + allocate (mat(sz, sz), prod(sz, sz)) + + call random_number(mat) + mat = 10.0_wp*(mat - 0.5_wp) + + !Since strictly diagonally dominant matrices are invertible, + !(see https://en.wikipedia.org/wiki/Diagonally_dominant_matrix), + !create a random diagonally dominant matrix. + do n = 1, sz + offdiag = 0.0_wp + do m = 1, sz + if (n == m) cycle + offdiag = offdiag + abs(mat(n, m)) + enddo + mat(n, n) = offdiag + 1.0_wp + enddo + + prod = matmul(mat, inverse(mat)) + + do n = 1, sz + if (abs(prod(n, n) - 1.0_wp) > tol*epsilon(1.0_wp)) then + allocate (error) + return + endif + enddo + + deallocate (mat, prod) + + end subroutine test_inv_re_rand + + subroutine test_inv_c_rand(error) + type(error_type), allocatable, intent(out) :: error + + complex(wp), allocatable :: mat(:, :), prod(:, :) + real(wp) :: szr, entryr, entryc + real(wp) :: offdiag + integer :: sz, n, m + + call random_seed() + call random_number(szr) + + sz = nint(50.0_wp*szr) + 5 + + allocate (mat(sz, sz), prod(sz, sz)) + + do n = 1, sz + do m = 1, sz + call random_number(entryr) + call random_number(entryc) + entryr = 10.0_wp*(entryr - 0.5_wp) + entryc = 10.0_wp*(entryc - 0.5_wp) + mat(n, m) = cmplx(entryr, entryc, wp) + enddo + enddo + + !Since strictly diagonally dominant matrices are invertible, + !(see https://en.wikipedia.org/wiki/Diagonally_dominant_matrix), + !create a random diagonally dominant matrix. + do n = 1, sz + offdiag = 0.0_wp + do m = 1, sz + if (n == m) cycle + offdiag = offdiag + abs(mat(n, m)) + enddo + mat(n, n) = cmplx(offdiag + 1.0_wp, 0.0_wp, wp) + enddo + + prod = matmul(mat, inverse(mat)) + + do n = 1, sz + if (abs(prod(n, n) - 1.0_wp) > tol*epsilon(1.0_wp)) then + allocate (error) + return + endif + enddo + + deallocate (mat, prod) + + end subroutine test_inv_c_rand + subroutine test_diagonalization_herm_rand(error) type(error_type), allocatable, intent(out) :: error diff --git a/test/suites/Extra_Utilities_Suite.F90 b/test/suites/Extra_Utilities_Suite.F90 index 33ce3db..e0e2edc 100644 --- a/test/suites/Extra_Utilities_Suite.F90 +++ b/test/suites/Extra_Utilities_Suite.F90 @@ -119,7 +119,6 @@ subroutine test_exp_and_log(error) call random_number(nr) n = nint(2.0_wp + 28.0_wp*nr) - n = 30 allocate (rnd1(n, n), rnd2(n, n), eig(n), m(n, n), res(n, n), p(n, n))