Giles method to differentiate solvers and minor edits to docs

BLAS/test/test_dtrsm.f90

@@ -0,0 +1,201 @@
+! Test program for DTRSM differentiation
+! Generated automatically by run_tapenade_blas.py
+! Using REAL*8 precision
+
+program test_dtrsm
+  implicit none
+
+  integer :: seed_array(33)
+
+  external :: dtrsm
+  external :: dtrsm_d
+
+  ! Test parameters
+  integer, parameter :: n = 4  ! Matrix/vector size for test
+  integer, parameter :: max_size = n  ! Maximum array dimension (rows/cols of matrices)
+  integer, parameter :: lda = max_size, ldb = max_size, ldc = max_size  ! Leading dimensions
+
+  character :: side
+  character :: uplo
+  character :: transa
+  character :: diag
+  integer :: msize
+  integer :: nsize
+  real(8) :: alpha
+  real(8), dimension(max_size,max_size) :: a
+  integer :: lda_val
+  real(8), dimension(max_size,max_size) :: b
+  integer :: ldb_val
+
+  ! Derivative variables
+  real(8) :: alpha_d
+  real(8), dimension(max_size,max_size) :: a_d
+  real(8), dimension(max_size,max_size) :: b_d
+
+  ! Storage variables for inout parameters
+  real(8), dimension(max_size,max_size) :: b_output
+
+  ! Array restoration variables for numerical differentiation
+  real(8) :: alpha_orig
+  real(8), dimension(max_size,max_size) :: b_orig
+  real(8), dimension(max_size,max_size) :: a_orig
+
+  ! Variables for central difference computation
+  real(8), dimension(max_size,max_size) :: b_forward, b_backward
+  ! Scalar variables for central difference computation
+  real(8) :: central_diff, ad_result
+  logical :: has_large_errors
+
+  ! Variables for storing original derivative values
+  real(8) :: alpha_d_orig
+  real(8), dimension(max_size,max_size) :: b_d_orig
+  real(8), dimension(max_size,max_size) :: a_d_orig
+
+  ! Temporary variables for matrix initialization
+  real(4) :: temp_real, temp_imag
+  integer :: i, j
+
+  ! Initialize random seed for reproducible results
+  seed_array = 42
+  call random_seed(put=seed_array)
+
+  side = 'L'
+  uplo = 'U'
+  transa = 'N'
+  diag = 'N'
+  msize = n
+  nsize = n
+  call random_number(alpha)
+  alpha = alpha * 2.0d0 - 1.0d0  ! Scale to [-1,1]
+  call random_number(a)
+  a = a * 2.0d0 - 1.0d0  ! Scale to [-1,1]
+  lda_val = lda  ! LDA must be at least max( 1
+  call random_number(b)
+  b = b * 2.0d0 - 1.0d0  ! Scale to [-1,1]
+  ldb_val = ldb
+
+  ! Initialize input derivatives to random values
+  call random_number(alpha_d)
+  alpha_d = alpha_d * 2.0e0 - 1.0e0  ! Scale to [-1,1]
+  call random_number(b_d)
+  b_d = b_d * 2.0e0 - 1.0e0  ! Scale to [-1,1]
+  call random_number(a_d)
+  a_d = a_d * 2.0e0 - 1.0e0  ! Scale to [-1,1]
+
+  ! Store initial derivative values after random initialization
+  alpha_d_orig = alpha_d
+  b_d_orig = b_d
+  a_d_orig = a_d
+
+  ! Store original values for central difference computation
+  alpha_orig = alpha
+  b_orig = b
+  a_orig = a
+
+  write(*,*) 'Testing DTRSM'
+  ! Store input values of inout parameters before first function call
+  b_orig = b
+
+  ! Re-initialize data for differentiated function
+  ! Only reinitialize inout parameters - keep input-only parameters unchanged
+
+  ! side already has correct value from original call
+  ! uplo already has correct value from original call
+  ! transa already has correct value from original call
+  ! diag already has correct value from original call
+  msize = n
+  nsize = n
+  ! alpha already has correct value from original call
+  ! a already has correct value from original call
+  lda_val = lda  ! LDA must be at least max( 1
+  b = b_orig
+  ldb_val = ldb
+
+  ! Call the differentiated function
+  call dtrsm_d(side, uplo, transa, diag, msize, nsize, alpha, alpha_d, a, a_d, lda_val, b, b_d, ldb_val)
+
+  ! Print results and compare
+  write(*,*) 'Function calls completed successfully'
+
+  ! Numerical differentiation check
+  call check_derivatives_numerically()
+
+  write(*,*) 'Test completed successfully'
+
+contains
+
+  subroutine check_derivatives_numerically()
+    implicit none
+    real(8), parameter :: h = 1.0e-6  ! Step size for finite differences
+    real(8) :: relative_error, max_error
+    real(8) :: output_orig, output_pert
+    real(8) :: numerical_result, analytical_result
+    real(8) :: abs_error, abs_reference, error_bound
+    integer :: i, j
+
+    max_error = 0.0e0
+    has_large_errors = .false.
+
+    write(*,*) 'Checking derivatives against numerical differentiation:'
+    write(*,*) 'Step size h =', h
+
+    ! Tolerance thresholds: rtol=1.0e-5, atol=1.0e-5
+
+    ! Original values already stored in main program
+
+    ! Central difference computation: f(x + h) - f(x - h) / (2h)
+    ! Forward perturbation: f(x + h)
+    alpha = alpha_orig + h * alpha_d_orig
+    b = b_orig + h * b_d_orig
+    a = a_orig + h * a_d_orig
+    call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val)
+    ! Store forward perturbation results
+    b_forward = b
+
+    ! Backward perturbation: f(x - h)
+    alpha = alpha_orig - h * alpha_d_orig
+    b = b_orig - h * b_d_orig
+    a = a_orig - h * a_d_orig
+    call dtrsm(side, uplo, transa, diag, msize, nsize, alpha, a, lda_val, b, ldb_val)
+    ! Store backward perturbation results
+    b_backward = b
+
+    ! Compute central differences and compare with AD results
+    ! Check derivatives for output B
+    do j = 1, min(2, n)  ! Check only first few elements
+      do i = 1, min(2, n)
+        ! Central difference: (f(x+h) - f(x-h)) / (2h)
+        central_diff = (b_forward(i,j) - b_backward(i,j)) / (2.0e0 * h)
+        ! AD result
+        ad_result = b_d(i,j)
+        ! Error check: |a - b| > atol + rtol * |b|
+        abs_error = abs(central_diff - ad_result)
+        abs_reference = abs(ad_result)
+        error_bound = 1.0e-5 + 1.0e-5 * abs_reference
+        if (abs_error > error_bound) then
+          has_large_errors = .true.
+          relative_error = abs_error / max(abs_reference, 1.0e-10)
+          write(*,*) 'Large error in output B(', i, ',', j, '):'
+          write(*,*) '  Central diff: ', central_diff
+          write(*,*) '  AD result:   ', ad_result
+          write(*,*) '  Absolute error:', abs_error
+          write(*,*) '  Error bound:', error_bound
+          write(*,*) '  Relative error:', relative_error
+        end if
+        ! Track max error for reporting (normalized)
+        relative_error = abs_error / max(abs_reference, 1.0e-10)
+        max_error = max(max_error, relative_error)
+      end do
+    end do
+
+    write(*,*) 'Maximum relative error:', max_error
+    write(*,*) 'Tolerance thresholds: rtol=1.0e-5, atol=1.0e-5'
+    if (has_large_errors) then
+      write(*,*) 'FAIL: Derivatives are outside tolerance'
+    else
+      write(*,*) 'PASS: Derivatives are within tolerance (rtol + atol)'
+    end if
+
+  end subroutine check_derivatives_numerically
+
+end program test_dtrsm