From 5a9cb37d226ad22d12ab3c18ba1c2558fa6e1aeb Mon Sep 17 00:00:00 2001 From: Nick Wogan Date: Fri, 29 Mar 2024 12:30:19 -0700 Subject: [PATCH] banded jacobian works --- src/forwarddiff_const.f90 | 1 + src/forwarddiff_derivative.f90 | 121 +++++++++++++++++++++++++++++---- test/CMakeLists.txt | 6 +- test/test_forwarddiff.f90 | 38 +++++------ test/test_sparse.f90 | 86 +++++++++++++++++++++++ 5 files changed, 217 insertions(+), 35 deletions(-) create mode 100644 test/test_sparse.f90 diff --git a/src/forwarddiff_const.f90 b/src/forwarddiff_const.f90 index 1b196f5..7083d18 100644 --- a/src/forwarddiff_const.f90 +++ b/src/forwarddiff_const.f90 @@ -1,3 +1,4 @@ module forwarddiff_const use iso_fortran_env, only: wp => real64 + implicit none end module \ No newline at end of file diff --git a/src/forwarddiff_derivative.f90 b/src/forwarddiff_derivative.f90 index 22fc7d2..8ed64e2 100644 --- a/src/forwarddiff_derivative.f90 +++ b/src/forwarddiff_derivative.f90 @@ -53,7 +53,7 @@ subroutine grad(fcn, x, f, dfdx, err) integer :: i if (size(x) /= size(dfdx)) then - err = 'Output dfdx array is not the right size.' + err = 'Output `dfdx` array is not the right size.' return endif @@ -78,47 +78,138 @@ subroutine grad(fcn, x, f, dfdx, err) end subroutine - subroutine jacobian(fcn, x, f, dfdx, err) + subroutine jacobian(fcn, x, f, dfdx, bandwidth, err) procedure(jacobian_sig) :: fcn real(wp), intent(in) :: x(:) real(wp), intent(out) :: f(:) real(wp), intent(out) :: dfdx(:,:) + integer, optional, intent(in) :: bandwidth character(:), allocatable :: err + ! real(wp), allocatable :: JacT(:,:) + integer :: bandwidth_, hbw type(dual) :: xx(size(x)) type(dual) :: ff(size(x)) - integer :: i + integer :: i, j, ii, jj, kk + + bandwidth_ = size(x) + if (present(bandwidth)) then + bandwidth_ = bandwidth + if (bandwidth_ > size(x)) then + err = '`bandwidth` can not be > size(x).' + return + endif + if (bandwidth_ < 1) then + err = '`bandwidth` can not be < 1.' + return + endif + if (mod(bandwidth_,2) == 0) then + err = '`bandwidth` must be odd.' + return + endif + endif if (size(x) /= size(f)) then - err = 'Output f array is not the right size.' + err = 'Output `f` array is not the right size.' return endif - if (size(x) /= size(dfdx,1) .or. size(x) /= size(dfdx,2)) then - err = 'Output dfdx array is not the right size.' + if (bandwidth_ /= size(dfdx,1) .or. size(x) /= size(dfdx,2)) then + err = 'Output `dfdx` array is not the right size.' return endif + ! allocate(JacT(size(x),bandwidth_)) + + ! Set x xx%val = x + + ! Allocate dual number do i = 1,size(x) - allocate(xx(i)%der(size(x))) + allocate(xx(i)%der(bandwidth_)) xx(i)%der = 0.0_wp - xx(i)%der(i) = 1.0_wp enddo + ! Seed the dual number + j = 1 + outer : do + do i = 1,bandwidth_ + if (j > size(x)) exit outer + xx(j)%der(i) = 1.0_wp + j = j + 1 + enddo + enddo outer + + do i = 1,size(x) + print*,xx(i)%der(:) + enddo + print*,"hi" + ! Do differentiation call fcn(xx, ff) ! Unpack f(x) f = ff%val - + ! Unpack jacobian - do i = 1,size(x) - dfdx(i,:) = ff(i)%der(:) - enddo + ! do j = 1,bandwidth_ + ! do i = 1,size(x) + ! JacT(i,j) = ff(i)%der(j) + ! ! dfdx(j,i) = ff(i)%der(j) + ! ! dfdx(1,:) is the first column of the jacobian + ! enddo + ! enddo + ! JacT(:,1) is the first column of the jacobian + + ! We now want to disentangle the jacobian + ! dfdx(1,:) is a row. It should contain diagonals. Row 1 is the highest diagonal + + ! First row of dfdx + ! dfdx(1,1:hbw) = 0.0_dp, unused + ! JacT(1,hbw+1) should be at dfdx(1,hbw+1) + ! JacT(2,hbw+2) should be at dfdx(1,hbw+2) + ! JacT(3,hbw) should be at dfdx(1,hbw+3) + ! JacT(4,hbw+1) should be at dfdx(1,hbw+4) + ! JacT(5,hbw+2) should be at dfdx(1,hbw+5) + ! JacT(6,hbw-1) should be at dfdx(1,hbw+6) + + ! Second row + + ! JacT(1,bw) is first diagonal entry + + if (present(bandwidth)) then + + hbw = (bandwidth_ - 1)/2 ! halfbandwidth + j = 1 + outer1 : do + do i = 1,bandwidth_ + if (j > size(x)) exit outer1 + do jj = -hbw,hbw + kk = jj + hbw + 1 + ii = j + jj + if (ii < 1) then + dfdx(kk,j) = 0.0_wp + elseif (ii > size(x)) then + dfdx(kk,j) = 0.0_wp + else + ! dfdx(kk,j) = JacT(ii,i) + dfdx(kk,j) = ff(ii)%der(i) + endif + enddo + j = j + 1 + enddo + enddo outer1 + + else + + ! Dense matrix + do j = 1,bandwidth_ + do i = 1,size(x) + dfdx(i,j) = ff(i)%der(j) + enddo + enddo + + endif end subroutine - - - end module \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index af4930c..663ffbe 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -1,3 +1,7 @@ add_executable(test_forwarddiff test_forwarddiff.f90) target_link_libraries(test_forwarddiff forwarddiff) -target_include_directories(test_forwarddiff PUBLIC ${CMAKE_Fortran_MODULE_DIRECTORY}) \ No newline at end of file +target_include_directories(test_forwarddiff PUBLIC ${CMAKE_Fortran_MODULE_DIRECTORY}) + +add_executable(test_sparse test_sparse.f90) +target_link_libraries(test_sparse forwarddiff) +target_include_directories(test_sparse PUBLIC ${CMAKE_Fortran_MODULE_DIRECTORY}) \ No newline at end of file diff --git a/test/test_forwarddiff.f90 b/test/test_forwarddiff.f90 index 98ec7d4..fb12b87 100644 --- a/test/test_forwarddiff.f90 +++ b/test/test_forwarddiff.f90 @@ -1,8 +1,9 @@ program test_forwarddiff use forwarddiff, only: wp, derivative, grad, jacobian implicit none + call test_dual() - call test_jacobian1() + call test_jacobian() contains @@ -43,7 +44,7 @@ subroutine test_dual() end subroutine function func_operators(x) result(res) - use forwarddiff_dual + use forwarddiff type(dual), intent(in) :: x type(dual) :: res @@ -56,7 +57,7 @@ function func_operators(x) result(res) end function function func_intrinsics1(x) result(res) - use forwarddiff_dual + use forwarddiff type(dual), intent(in) :: x type(dual) :: res @@ -72,7 +73,7 @@ function func_intrinsics1(x) result(res) end function function func_intrinsics2(x) result(res) - use forwarddiff_dual + use forwarddiff type(dual), intent(in) :: x type(dual) :: res @@ -89,48 +90,48 @@ function func_intrinsics2(x) result(res) end function function func_grad1(x) result(res) - use forwarddiff_dual + use forwarddiff type(dual), intent(in) :: x(:) type(dual) :: res res = x(1)*x(1)*x(2) + x(1) + x(2) end function function func_grad2(x) result(res) - use forwarddiff_dual + use forwarddiff type(dual), intent(in) :: x(:) type(dual) :: res res = sum(x*3.14_wp) end function - subroutine test_jacobian1() + subroutine test_jacobian() real(wp) :: u(3), f(3), dfdu(3,3) real(wp) :: f1(3), dfdu1(3,3) character(:), allocatable :: err u = [1.0_wp, 2.0_wp, 3.0_wp] - call jacobian(rhs_rober_dual, u, f, dfdu, err) + call jacobian(rhs_rober_dual, u, f, dfdu, err=err) - print*,f + ! print*,f print*,'' - print*,dfdu(:,1) - print*,dfdu(:,2) - print*,dfdu(:,3) + print*,dfdu(1,:) + print*,dfdu(2,:) + print*,dfdu(3,:) ! Check against analytical solution call rhs_rober(u, f1) call jac_rober(u, dfdu1) if (.not. all(f == f1)) then - error stop "test_jacobian1 failed" + error stop "test_jacobian failed" endif if (.not. all(dfdu == dfdu1)) then - error stop "test_jacobian1 failed" + error stop "test_jacobian failed" endif end subroutine subroutine rhs_rober_dual(u, du) - use forwarddiff_dual + use forwarddiff type(dual), intent(in) :: u(:) type(dual), intent(out) :: du(:) @@ -145,7 +146,6 @@ subroutine rhs_rober_dual(u, du) end subroutine subroutine rhs_rober(u, du) - use forwarddiff_dual real(wp), intent(in) :: u(:) real(wp), intent(out) :: du(:) @@ -167,9 +167,9 @@ subroutine jac_rober(u, pd) k2 = 3.0e7_wp, & k3 = 1.0e4_wp - pd(:,1) = [-k1, k1, 0.0_wp] - pd(:,2) = [k3*u(3), -2.0_wp*k2*u(2) - k3*u(3), 2.0_wp*k2*u(2)] - pd(:,3) = [k3*u(2), -k3*u(2), 0.0_wp] + pd(:,1) = [-k1, k1, 0.0_wp] ! column 1 + pd(:,2) = [k3*u(3), -2.0_wp*k2*u(2) - k3*u(3), 2.0_wp*k2*u(2)] ! column 2 + pd(:,3) = [k3*u(2), -k3*u(2), 0.0_wp] ! column 3 end subroutine diff --git a/test/test_sparse.f90 b/test/test_sparse.f90 new file mode 100644 index 0000000..f5a2190 --- /dev/null +++ b/test/test_sparse.f90 @@ -0,0 +1,86 @@ + +program main + use forwarddiff, only: wp, jacobian + implicit none + integer, parameter :: nz = 7 + + call test_banded() + +contains + + subroutine test_banded() + integer, parameter :: bandwidth = 3 + real(wp) :: u(nz), f(nz), dfdu(bandwidth,nz) + real(wp) :: f1(nz), dfdu1(nz,nz) + character(:), allocatable :: err + integer :: i + + do i = 1,nz + u(i) = i + enddo + call jacobian(rhs_banded_dual, u, f, dfdu, bandwidth=bandwidth, err=err) + if (allocated(err)) then + print*,err + stop + endif + + call rhs_banded(u, f1) + call jac_banded(u, dfdu1) + + do i = 1,bandwidth + print*,dfdu(i,:) + enddo + print*,'' + + do i = 1,nz + print*,dfdu1(i,:) + enddo + + end subroutine + + subroutine rhs_banded_dual(u, du) + use forwarddiff + type(dual), intent(in) :: u(:) + type(dual), intent(out) :: du(:) + integer :: i + + du(1) = 3*u(2) - u(1) + do i = 2,nz-1 + du(i) = 3*u(i+1) - 2.0_wp*u(i) + u(i-1) + enddo + du(nz) = - u(nz) + u(nz-1) + + end subroutine + + subroutine rhs_banded(u, du) + real(wp), intent(in) :: u(:) + real(wp), intent(out) :: du(:) + integer :: i + + du(1) = 3*u(2) - u(1) + do i = 2,nz-1 + du(i) = 3*u(i+1) - 2.0_wp*u(i) + u(i-1) + enddo + du(nz) = - u(nz) + u(nz-1) + + end subroutine + + subroutine jac_banded(u, pd) + real(wp), intent(in) :: u(:) + real(wp), intent(out) :: pd(:,:) + integer :: i + + pd = 0.0_wp + pd(1,1) = -1.0_wp + pd(2,1) = 1.0_wp + do i = 2,nz-1 + pd(i,i) = -2.0_wp + pd(i+1,i) = 1.0_wp + pd(i-1,i) = 3.0_wp + enddo + pd(nz,nz) = -1.0_wp + pd(nz-1,nz) = 3.0_wp + + end subroutine + +end program \ No newline at end of file