Add tests

Author: jchristopherson

tests/nonlin_test.f90

@@ -18,6 +18,9 @@ program main
     rst = test_jacobian_1()
     if (.not.rst) overall = .false.
 
+    rst = test_jacobian_2()
+    if (.not.rst) overall = .false.
+
     rst = test_quasinewton_1()
     if (.not.rst) overall = .false.
 

tests/nonlin_test_jacobian.f90

@@ -8,6 +8,7 @@ module nonlin_test_jacobian
     implicit none
     private
     public :: test_jacobian_1
+    public :: test_jacobian_2
 
 contains
 ! ******************************************************************************
@@ -44,6 +45,46 @@ subroutine jac1(x, jac, args)
             [2, 2])
     end subroutine
 
+! ------------------------------------------------------------------------------
+    ! FUNCTIONS:
+    ! x = Y * r * cos(theta)
+    ! y = Y * r * sin(theta)
+    subroutine fcn2(x, f, args)
+        real(real64), intent(in), dimension(:) :: x
+        real(real64), intent(out), dimension(:) :: f
+        class(*), intent(inout), optional :: args
+        real(real64) :: r, theta, y
+        r = x(1)
+        theta = x(2)
+        select type (args)
+        type is (real(real64))
+            y = args
+        end select
+        f(1) = y * r * cos(theta) 
+        f(2) = y * r * sin(theta)
+    end subroutine
+
+    ! JACOBIAN:
+    !     |Y * cos  -Y * r * sin|
+    ! J = |                     |
+    !     |Y * sin   Y * r * cos|
+    subroutine jac2(x, jac, args)
+        real(real64), intent(in), dimension(:) :: x
+        real(real64), intent(out), dimension(:,:) :: jac
+        class(*), intent(inout), optional :: args
+        real(real64) :: r, theta, y
+        r = x(1)
+        theta = x(2)
+        select type (args)
+        type is (real(real64))
+            y = args
+        end select
+        jac = reshape(&
+            [y * cos(theta), y * sin(theta), &
+            -y * r * sin(theta), y * r * cos(theta)], &
+            [2, 2])
+    end subroutine
+
 ! ------------------------------------------------------------------------------
     function test_jacobian_1() result(rst)
         ! Local Variables
@@ -135,5 +176,97 @@ function test_jacobian_1() result(rst)
 100     format(A)
     end function
 
+! ------------------------------------------------------------------------------
+    function test_jacobian_2() result(rst)
+        ! Local Variables
+        type(vecfcn_helper) :: obj
+        procedure(vecfcn), pointer :: fcn
+        procedure(jacobianfcn), pointer :: jac
+        real(real64) :: numjac(2, 2), exact(2, 2), x(2), y
+        logical :: rst
+        integer(int32) :: i
+
+        ! Parameters
+        real(real64), parameter :: tol = 1.0d-4
+
+        ! Initialization
+        rst = .true.
+        fcn => fcn2
+        jac => jac2
+        call obj%set_fcn(fcn, 2, 2)
+        call random_number(y)
+
+        ! Test point 1 (0, 0)
+        x = 0.0d0
+        call obj%jacobian(x, numjac, args = y)
+        call jac(x, exact, args = y)
+        if (.not.assert(numjac, exact, tol)) then
+            rst = .false.
+            print 100, "Test Failed: Jacobian Test 2A"
+            print 100, "Expected:"
+            do i = 1, size(exact, 1)
+                print *, exact(i,:)
+            end do
+            print 100, "Found:"
+            do i = 1, size(exact, 1)
+                print *, numjac(i,:)
+            end do
+        end if
+
+        ! Test point 2 (1, 0)
+        x = [1.0d0, 0.0d0]
+        call obj%jacobian(x, numjac, args = y)
+        call jac(x, exact, args = y)
+        if (.not.assert(numjac, exact, tol)) then
+            rst = .false.
+            print 100, "Test Failed: Jacobian Test 2B"
+            print 100, "Expected:"
+            do i = 1, size(exact, 1)
+                print *, exact(i,:)
+            end do
+            print 100, "Found:"
+            do i = 1, size(exact, 1)
+                print *, numjac(i,:)
+            end do
+        end if
+
+        ! Test point 3 (0, 1)
+        x = [0.0d0, 1.0d0]
+        call obj%jacobian(x, numjac, args = y)
+        call jac(x, exact, args = y)
+        if (.not.assert(numjac, exact, tol)) then
+            rst = .false.
+            print 100, "Test Failed: Jacobian Test 2C"
+            print 100, "Expected:"
+            do i = 1, size(exact, 1)
+                print *, exact(i,:)
+            end do
+            print 100, "Found:"
+            do i = 1, size(exact, 1)
+                print *, numjac(i,:)
+            end do
+        end if
+
+        ! Test point 4 (0.5, -0.5)
+        x = [0.5d0, -0.5d0]
+        call obj%jacobian(x, numjac, args = y)
+        call jac(x, exact, args = y)
+        if (.not.assert(numjac, exact, tol)) then
+            rst = .false.
+            print 100, "Test Failed: Jacobian Test 2D"
+            print 100, "Expected:"
+            do i = 1, size(exact, 1)
+                print *, exact(i,:)
+            end do
+            print 100, "Found:"
+            do i = 1, size(exact, 1)
+                print *, numjac(i,:)
+            end do
+        end if
+
+        ! Formatting
+100     format(A)
+    end function
+
 ! ------------------------------------------------------------------------------
 end module