banded jacobian works

src/forwarddiff_const.f90

@@ -1,3 +1,4 @@
 module forwarddiff_const
   use iso_fortran_env, only: wp => real64
+  implicit none
 end module

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

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

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
 

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