lalg_adv_tests.F90

Go to the documentation of this file.

!! Copyright (C) 2025 N. Tancogne-Dejean
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include "global.h"
 
module lalg_adv_tests_oct_m
  use debug_oct_m
  use global_oct_m
  use iso_c_binding
  use, intrinsic :: iso_fortran_env
  use lalg_adv_oct_m
  use profiling_oct_m
  use messages_oct_m
  use namespace_oct_m
 
  implicit none
 
  private
  public ::                    &
    test_exponential_matrix
 
contains
 
  !---------------------------------------------------------------------------
  subroutine test_exponential_matrix(namespace)
    type(namespace_t), intent(in) :: namespace
 
    ! Tolerance for floating-point comparisons
    real(real64), parameter :: tol = 1.0e-12_real64
    integer, parameter :: n = 2
    real(real64) :: factor
    real(real64), dimension(n, n) :: A, funA, expected
    complex(real64), dimension(n, n) :: zA, zfunA, zexpected
    logical :: hermitian
    integer, parameter :: n3 = 3
    complex(real64), dimension(n3, n3) :: A3, funA3, ident3, product
    complex(real64) :: zfactor
 
    push_sub(test_exponential_matrix)
 
    !-----------------------------------------------------------------------
    ! Test 1: Diagonal matrices
    ! A = [ [1, 0], [0, 2] ]
    ! so that exp(A) should be diag(exp(1), exp(2))
    ! and exp(iA) should be [ [exp(i), 1], [1, exp(2i)] ]
    !-----------------------------------------------------------------------
    a = reshape([1.0_real64, 0.0_real64, 0.0_real64, 2.0_real64], shape(a))
    factor = 1.0_real64
    hermitian = .true.
    call lalg_matrix_function(n, factor, a, funa, expfun, hermitian)
 
    ! Build expected result: only diagonal entries nonzero
    expected = 0.0_real64
    expected(1, 1) = exp(1.0_real64)
    expected(2, 2) = exp(2.0_real64)
 
    message(1) = "Testing real diagonal matrix exponential..."
    if ( maxval( abs(funa - expected) ) < tol ) then
      message(2) = "Diagonal matrix test passed."
    else
      message(2) = "Diagonal matrix test failed."
    endif
    call messages_info(2, namespace=namespace)
 
 
    za = reshape([m_z1, m_z0, m_z0, m_two*m_z1], shape(a))
    zfactor = m_zi
 
    ! Build expected result: only diagonal entries nonzero
    zexpected = m_zero
    zexpected(1, 1) = exp(m_zi)
    zexpected(2, 2) = exp(m_two*m_zi)
 
    ! Assuming general matrix first
    hermitian = .false.
    call lalg_matrix_function(n, zfactor, za, zfuna, zexpfun, hermitian)
 
    message(1) = "Testing complex diagonal matrix exponential..."
    if ( maxval( abs(zfuna - zexpected) ) < tol ) then
      message(2) = "Diagonal matrix test passed."
    else
      message(2) = "Diagonal matrix test failed."
    endif
    call messages_info(2, namespace=namespace)
 
    ! Going via the Hermitian patch
    hermitian = .true.
    call lalg_matrix_function(n, zfactor, za, zfuna, zexpfun, hermitian)
 
    message(1) = "Testing complex diagonal Hermitian matrix exponential..."
    if ( maxval( abs(zfuna - zexpected) ) < tol ) then
      message(2) = "Diagonal matrix test passed."
    else
      message(2) = "Diagonal matrix test failed."
    endif
    call messages_info(2, namespace=namespace)
 
    !-----------------------------------------------------------------------
    ! Test 2: Nilpotent matrix
    ! A = [ [0, 1+i], [0, 0] ] where A**2 = 0 and exp(A) = I + A
    !-----------------------------------------------------------------------
    za = m_z0
    za(1,2) = m_one + m_zi
    ! Use a factor of 1 (i.e. compute exp(1*A))
    zfactor = m_one
    hermitian = .false.
    call lalg_matrix_function(n, zfactor, za, zfuna, zexpfun, hermitian)
 
    ! Build expected result: I + A
    zexpected = (0.0_real64, 0.0_real64)
    zexpected(1,1) = (1.0_real64, 0.0_real64)
    zexpected(2,2) = (1.0_real64, 0.0_real64)
    zexpected(1,2) = (1.0_real64, 1.0_real64)
 
    message(1) = "Testing nilpotent matrix exponential..."
    if ( maxval( abs(zfuna - zexpected) ) < tol ) then
      message(2) = "Nilpotent matrix test passed."
    else
      message(2) = "Nilpotent matrix test failed."
    endif
    call messages_info(2, namespace=namespace)
 
    !-----------------------------------------------------------------------
    ! Test 3: Tridiagonal matrix with complex factor: exp(iA)
    ! We use a 3x3 symmetric tridiagonal matrix:
    !   A = [ 2  1  0 ]
    !       [ 1  2  1 ]
    !       [ 0  1  2 ]
    ! With factor = i, we compute exp(iA). For hermitian A, exp(iA) is unitary.
    ! We verify unitarity by checking that M * M^dagger = I.
    !-----------------------------------------------------------------------
    hermitian = .true.
    zfactor = m_zi
 
    a3 = (0.0_real64, 0.0_real64)
    ! Set the main diagonal to 2
    a3(1,1) = (2.0_real64, 0.0_real64)
    a3(2,2) = (2.0_real64, 0.0_real64)
    a3(3,3) = (2.0_real64, 0.0_real64)
    ! Set the off-diagonals to 1
    a3(1,2) = (1.0_real64, 0.0_real64); a3(2,1) = (1.0_real64, 0.0_real64)
    a3(2,3) = (1.0_real64, 0.0_real64); a3(3,2) = (1.0_real64, 0.0_real64)
 
    call lalg_matrix_function(n3, zfactor, a3, funa3, zexpfun, hermitian, tridiagonal=.true.)
 
    ! Build the 3x3 identity matrix
    ident3 = (0.0_real64, 0.0_real64)
    ident3(1,1) = (1.0_real64, 0.0_real64)
    ident3(2,2) = (1.0_real64, 0.0_real64)
    ident3(3,3) = (1.0_real64, 0.0_real64)
 
    ! Check unitarity: funA3 * funA3^dagger should equal identity.
    product = matmul(funa3, transpose(conjg(funa3)))
 
    message(1) = "Testing tridiagonal matrix exp(iA) for unitarity..."
    if ( maxval( abs(product - ident3) ) < tol ) then
      message(2) = "Tridiagonal matrix exp(iA) test passed."
    else
      message(2) = "Tridiagonal matrix exp(iA) test failed."
    endif
    call messages_info(2, namespace=namespace)
 
    pop_sub(test_exponential_matrix)
  end subroutine test_exponential_matrix
 
  !-----------------------------------------------------------------------
  real(real64) function expfun(z)
    implicit none
    real(real64), intent(in) :: z
    expfun = exp(z)
  end function expfun
 
  !-----------------------------------------------------------------------
  complex(real64) function zexpfun(z)
    implicit none
    complex(real64), intent(in) :: z
    zexpfun = exp(z)
  end function zexpfun
 
end module lalg_adv_tests_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: