Include matrix inversion utility.

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

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)") " ____    __    ____  ___      .__   __. .__   __.  __  .__   __. .___________. "

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 = "[email protected]"

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

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

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
 

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))