Merge pull request #81 from kousuke-nakano/devel-#12

Progress in refactoring code and implementing unit tests

CMakeLists.txt

@@ -577,8 +577,14 @@ if(EXT_MODTEST)
   add_target_exe_serial_wrapper(test_zgetrf test_tools)
   add_target_exe_serial_wrapper(test_dgetrfi test_tools)
   add_target_exe_serial_wrapper(test_zgetrfi test_tools)
+  add_target_exe_serial_wrapper(test_zsktri test_tools)
+  add_target_exe_serial_wrapper(test_zsktrs test_tools)
   add_target_exe_serial_wrapper(test_openmp_reduction_real test_tools)
   add_target_exe_serial_wrapper(test_openmp_reduction_complex test_tools)
+  add_target_exe_serial_wrapper(test_upwinvp test_tools)
+  add_target_exe_serial_wrapper(test_upwinvp_complex test_tools)
+  #add_target_exe_serial_wrapper(test_upwinvp_pfaff test_tools)
+  #add_target_exe_serial_wrapper(test_upwinvp_pfaff_complex test_tools)
 
   if(EXT_GPU)
 

src/a_module_tests/CMakeLists.txt

@@ -162,6 +162,44 @@ foreach( EXECUTABLE IN LISTS EXECUTABLES_S_L
       )
   endif()
 
+  ######################################################################
+  #
+  # Add sources and link libraries for zsktr[ri] wrapper tests
+  #
+
+  if(${EXECUTABLE} MATCHES test_zsktr[si])
+      set_target_properties( ${EXECUTABLE} PROPERTIES SUFFIX ".x")
+
+    target_sources( ${EXECUTABLE}
+      PRIVATE
+      ${EXECUTABLE}.f90
+      )
+
+    target_link_libraries( ${EXECUTABLE}
+      PRIVATE
+        common-serial
+      )
+  endif()
+
+  ######################################################################
+  #
+  # Add sources and link libraries for upwinvp wrapper tests
+  #
+
+  if(${EXECUTABLE} MATCHES test_upwinvp.*)
+      set_target_properties( ${EXECUTABLE} PROPERTIES SUFFIX ".x")
+
+    target_sources( ${EXECUTABLE}
+      PRIVATE
+      ${EXECUTABLE}.f90
+      )
+
+    target_link_libraries( ${EXECUTABLE}
+      PRIVATE
+        common-serial
+      )
+  endif()
+
   ######################################################################
   #
   #

src/a_module_tests/test_upwinvp.f90

@@ -0,0 +1,56 @@
+! Copyright (C) 2022 TurboRVB group
+!
+! 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 3 of the License, 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, see <http://www.gnu.org/licenses/>.
+
+program test_upwinvp
+    use constants, only: yes_ontarget
+    implicit none
+    integer, parameter :: nel = 3, indt = 4
+    real*8 :: psi(indt, nel, 2), ainv(nel), ainvn(nel), winv(nel, indt)
+    integer :: i, j, k
+    real*8, parameter :: eps = 1.0e-10 ! A small value for numerical comparison
+
+    yes_ontarget = .true.
+
+    ! sleep to avoid unexpected corruptions.
+    call sleep(2)
+
+    ! Initialize test data
+    winv = 0.0
+    do i = 1, indt
+        do j = 1, nel
+            do k = 1, 2
+                psi(i, j, k) = i*j*k
+            end do
+        end do
+    end do
+    ainv = (/1.0, 2.0, 3.0/)
+    ainvn = (/4.0, 5.0, 6.0/)
+
+    call upwinvp(nel, indt, winv, ainv, ainvn, psi)
+
+    ! Check the result
+    do i = 1, indt
+        do j = 1, nel
+            if (abs(winv(j, i) - (ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2))) > eps) then
+                print *, "Test failed at i=", i, ", j=", j
+                print *, "ERROR"
+                stop 1
+            end if
+        end do
+    end do
+
+    print *, "OK"
+    print *, "All tests passed."
+end program test_upwinvp

src/a_module_tests/test_upwinvp_complex.f90

@@ -0,0 +1,63 @@
+! Copyright (C) 2022 TurboRVB group
+!
+! 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 3 of the License, 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, see <http://www.gnu.org/licenses/>.
+
+program test_upwinvp_complex
+    use constants, only: yes_ontarget
+    implicit none
+    integer, parameter :: nel = 3, indt = 4
+    complex*16 :: psi(indt, nel, 2), ainv(nel), ainvn(nel), winv(nel, indt), expected(nel, indt)
+    integer :: i, j, k
+    real*8, parameter :: eps = 1.0e-10 ! A small value for numerical comparison
+
+    yes_ontarget = .true.
+
+    ! sleep to avoid unexpected corruptions.
+    call sleep(1)
+
+    ! Initialize test data
+    winv = cmplx(0.0, 0.0)
+    do i = 1, indt
+        do j = 1, nel
+            do k = 1, 2
+                psi(i, j, k) = cmplx(i*j*k, i*j*k)
+            end do
+        end do
+    end do
+    ainv = (/cmplx(1.0, 1.0), cmplx(2.0, 2.0), cmplx(3.0, 3.0)/)
+    ainvn = (/cmplx(4.0, 4.0), cmplx(5.0, 5.0), cmplx(6.0, 6.0)/)
+
+    call upwinvp_complex(nel, indt, winv, ainv, ainvn, psi)
+
+    ! Calculate expected result
+    do i = 1, indt
+        do j = 1, nel
+            expected(j, i) = ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2)
+        end do
+    end do
+
+    ! Check the result
+    do i = 1, indt
+        do j = 1, nel
+            if (abs(winv(j, i) - expected(j, i)) > eps) then
+                print *, "Test failed at i=", i, ", j=", j
+                print *, "OK"
+                stop 1
+            end if
+        end do
+    end do
+
+    print *, "OK"
+    print *, "All tests passed."
+end program test_upwinvp_complex

src/a_module_tests/test_zsktri.f90

@@ -0,0 +1,96 @@
+! Copyright (C) 2022 TurboRVB group
+!
+! 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 3 of the License, 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, see <http://www.gnu.org/licenses/>.
+
+program test_zsktri
+    implicit none
+
+    complex*16, allocatable, dimension(:, :) :: A, A_inv, A_inv_orig
+    complex*16, allocatable, dimension(:) :: W
+    integer, allocatable, dimension(:) :: ipiv
+    real*8, allocatable, dimension(:, :) :: helper_r, helper_c
+    complex*16 :: one = 1.d0, zero = 0.d0
+    integer :: s, gen, ii, jj, info
+    character(len=1) :: uplo
+
+    ! s dimension of the test matrix, s x s.
+    ! uplo: U->upper triangular, L->lower triangular
+    ! gen = 0 : Compare matrices, gen = 1 : Generate matrices
+    write (*, *) 's, uplo, gen'
+    read (*, *) s, uplo, gen
+
+    allocate (ipiv(s))
+    allocate (A(s, s))
+    allocate (A_inv(s, s))
+    allocate (A_inv_orig(s, s))
+    allocate (W(s**2 + 12*s - 2))
+    allocate (helper_r(s, s))
+    allocate (helper_c(s, s))
+
+    if (gen .eq. 1) then
+
+        ! generate a skew_symmetric matrix A
+        call random_number(helper_r)
+        call random_number(helper_c)
+
+        do ii = 1, s
+            do jj = 1, ii - 1
+                helper_c(jj, ii) = -helper_c(ii, jj)
+                helper_r(jj, ii) = -helper_r(ii, jj)
+            end do
+        end do
+
+        do ii = 1, s
+            helper_r(ii, ii) = 0
+            helper_c(ii, ii) = 0
+        end do
+
+        A = cmplx(helper_r, helper_c)
+
+    else
+        open (unit=10, form="unformatted", file="A", action="read")
+        read (10) A
+        close (10)
+        open (unit=10, form="unformatted", file="A_inv_orig", action="read")
+        read (10) A_inv_orig
+        close (10)
+    end if
+
+    W = 0
+
+    do ii = 1, s
+        ipiv(ii) = ii
+    end do
+
+    call zsktri(uplo, s, A, s, A_inv, s, ipiv, W, info)
+
+    if (gen .eq. 1) then
+        open (unit=10, form="unformatted", file="A", action="write")
+        write (10) A
+        close (10)
+        open (unit=10, form="unformatted", file="A_inv_orig", action="write")
+        write (10) A_inv_orig
+        close (10)
+    else
+        A_inv = A_inv - A_inv_orig
+        if (maxval(abs(A_inv)) > 1.0d-10) then
+            print *, "ERROR"
+        else
+            print *, "OK"
+        end if
+    end if
+
+    deallocate (A, A_inv, A_inv_orig, W, helper_r, helper_c)
+
+end program test_zsktri

src/a_module_tests/test_zsktrs.f90

@@ -0,0 +1,77 @@
+! Copyright (C) 2022 TurboRVB group
+!
+! 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 3 of the License, 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, see <http://www.gnu.org/licenses/>.
+
+program test_zsktrs
+    implicit none
+
+    complex*16, allocatable, dimension(:, :) :: B, B_orig, X
+    complex*16, allocatable, dimension(:) :: A
+    real*8, allocatable, dimension(:) :: helper_r, helper_c
+    complex*16 :: one = 1.d0, zero = 0.d0
+    integer :: s, gen, ii, jj, info
+    character(len=1) :: uplo
+
+    ! s dimension of the test matrix, s x s.
+    ! uplo: U->upper triangular, L->lower triangular
+    ! gen = 0 : Compare matrices, gen = 1 : Generate matrices
+    write (*, *) 's, uplo, gen'
+    read (*, *) s, uplo, gen
+
+    allocate (A(s - 1))
+    allocate (B(s, s))
+    allocate (B_orig(s, s))
+    allocate (X(s, s))
+    allocate (helper_r(s - 1))
+    allocate (helper_c(s - 1))
+
+    if (gen .eq. 1) then
+
+        ! generate tridiagonal skewsymmetric matrix A
+        call random_number(helper_r)
+        call random_number(helper_c)
+        A = cmplx(helper_r, helper_c)
+
+    else
+        open (unit=10, form="unformatted", file="A", action="read")
+        read (10) A
+        close (10)
+        open (unit=10, form="unformatted", file="B_orig", action="read")
+        read (10) B_orig
+        close (10)
+    end if
+
+    X = 0
+
+    call zsktrs(uplo, s, s, A, B, s, X, info)
+
+    if (gen .eq. 1) then
+        open (unit=10, form="unformatted", file="A", action="write")
+        write (10) A
+        close (10)
+        open (unit=10, form="unformatted", file="B_orig", action="write")
+        write (10) B_orig
+        close (10)
+    else
+        B = B - B_orig
+        if (maxval(abs(B)) > 1.0d-10) then
+            print *, "ERROR"
+        else
+            print *, "OK"
+        end if
+    end if
+
+    deallocate (A, B, B_orig, X, helper_r, helper_c)
+
+end program test_zsktrs

src/m_common/upwinv.f90

@@ -16,8 +16,15 @@
 subroutine upwinv(nel, jel, indt, nelorb, winv, v, psi)
     use constants, only: yes_ontarget
     implicit none
-    integer nel, indt, i, j, jel, nelorb
-    real*8 winv(nel, indt), psi(indt, nel), v(nel)
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt, jel, nelorb
+    real*8, intent(in) :: psi(indt, nel), v(nel)
+    real*8, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer i, j
+
 #ifdef _OFFLOAD
 !$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
 #endif
@@ -43,11 +50,18 @@ subroutine upwinv(nel, jel, indt, nelorb, winv, v, psi)
     end do
     return
 end
+
 subroutine upwinv_complex(nel, jel, indt, nelorb, winv, v, psi)
     use constants, only: yes_ontarget
     implicit none
-    integer nel, indt, i, j, jel, nelorb
-    complex*16 winv(nel, indt), psi(indt, nel), v(nel)
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt, jel, nelorb
+    complex*16, intent(in) :: psi(indt, nel), v(nel)
+    complex*16, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer i, j
 #ifdef _OFFLOAD
 !$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
 #endif