Record propagator join normalization factors

NEO-2-QL/join_diagnostics_mod.f90

@@ -3,6 +3,10 @@ module join_diagnostics_mod
     implicit none
     private
 
+    character(len=:), allocatable, save :: normalization_output_filename
+    integer, save :: normalization_sequence = 0
+    logical, save :: normalization_output_initialized = .false.
+
     public :: report_join_failure, validate_join_normalization
 
 contains
@@ -33,7 +37,8 @@ subroutine validate_join_normalization(facnorm, direction, column, old_tags, &
         integer, intent(out) :: ierr
 
         if (ieee_is_finite(facnorm) .and. facnorm /= 0.0d0) then
-            ierr = 0
+            call record_join_normalization(facnorm, direction, column, &
+                old_tags, new_tags, npass_l, npass_r, ierr)
             return
         end if
 
@@ -45,4 +50,67 @@ subroutine validate_join_normalization(facnorm, direction, column, old_tags, &
         write(*, '(a,2(1x,i0))') 'join_ripples: pass dimensions=', npass_l, npass_r
         ierr = 5
     end subroutine validate_join_normalization
+
+    subroutine record_join_normalization(facnorm, direction, column, old_tags, &
+            new_tags, npass_l, npass_r, ierr)
+        real(kind(1.0d0)), intent(in) :: facnorm
+        character(len=1), intent(in) :: direction
+        integer, intent(in) :: column, old_tags(2), new_tags(2), npass_l, npass_r
+        integer, intent(out) :: ierr
+
+        integer :: iunit, status
+
+        ierr = 0
+        if (.not. normalization_output_initialized) &
+            call initialize_normalization_output(ierr)
+        if (ierr /= 0 .or. .not. allocated(normalization_output_filename)) return
+        open(newunit=iunit, file=normalization_output_filename, status='old', &
+            position='append', action='write', iostat=status)
+        if (status /= 0) then
+            write(*, '(a,i0)') 'join_ripples: normalization output error=', status
+            ierr = 6
+            return
+        end if
+        normalization_sequence = normalization_sequence + 1
+        write(iunit, &
+            '(i0,",",a,",",7(i0,","),es25.16e3)', iostat=status) &
+            normalization_sequence, direction, column, old_tags, new_tags, &
+            npass_l, npass_r, facnorm
+        close(iunit)
+        if (status /= 0) then
+            write(*, '(a,i0)') 'join_ripples: normalization output error=', status
+            ierr = 6
+        end if
+    end subroutine record_join_normalization
+
+    subroutine initialize_normalization_output(ierr)
+        integer, intent(out) :: ierr
+
+        integer :: iunit, length, status
+
+        ierr = 0
+        normalization_output_initialized = .true.
+        call get_environment_variable('NEO2_JOIN_NORMALIZATION_FILE', &
+            length=length, status=status)
+        if (status /= 0 .or. length == 0) return
+        allocate(character(len=length) :: normalization_output_filename)
+        call get_environment_variable('NEO2_JOIN_NORMALIZATION_FILE', &
+            value=normalization_output_filename, status=status)
+        if (status == 0) then
+            open(newunit=iunit, file=normalization_output_filename, &
+                status='replace', action='write', iostat=status)
+        end if
+        if (status == 0) then
+            write(iunit, '(a)', iostat=status) &
+                'sequence,direction,column,old_start,old_end,new_start,' // &
+                'new_end,npass_l,npass_r,factor'
+            close(iunit)
+        end if
+        if (status /= 0) then
+            write(*, '(a,i0)') 'join_ripples: normalization output error=', status
+            if (allocated(normalization_output_filename)) &
+                deallocate(normalization_output_filename)
+            ierr = 6
+        end if
+    end subroutine initialize_normalization_output
 end module join_diagnostics_mod

TEST/CMakeLists.txt

@@ -133,6 +133,7 @@ add_test(NAME join_failure_diagnostic_test
          WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
 
 set_tests_properties(join_failure_diagnostic_test PROPERTIES
+    ENVIRONMENT "NEO2_JOIN_NORMALIZATION_FILE=${CMAKE_CURRENT_BINARY_DIR}/join_normalization.csv"
     PASS_REGULAR_EXPRESSION "join_ripples: DGBSV info=2"
     FAIL_REGULAR_EXPRESSION "FAIL"
     TIMEOUT 30

TEST/test_join_failure_diagnostic.f90

@@ -5,7 +5,12 @@ program test_join_failure_diagnostic
     implicit none
 
     real(kind(1.0d0)) :: matrix(2, 2), rhs(2, 1)
-    integer :: info, ierr, pivot(2)
+    real(kind(1.0d0)) :: factor
+    character(len=1024) :: filename
+    character(len=128) :: header
+    character(len=1) :: direction
+    integer :: column, info, ierr, iunit, new_tags(2), old_tags(2), &
+        npass(2), pivot(2), sequence, status
 
     matrix = reshape([1.0d0, 2.0d0, 2.0d0, 4.0d0], shape(matrix))
     rhs(:, 1) = [1.0d0, 2.0d0]
@@ -23,4 +28,19 @@ program test_join_failure_diagnostic
     call validate_join_normalization(1.0d0, 'p', 2, [12, 13], [13, 13], &
         24, 24, ierr)
     if (ierr /= 0) error stop 'FAIL: finite normalization was rejected'
+
+    call get_environment_variable('NEO2_JOIN_NORMALIZATION_FILE', filename)
+    open(newunit=iunit, file=trim(filename), status='old', action='read', &
+        iostat=status)
+    if (status /= 0) error stop 'FAIL: join normalization output is absent'
+    read(iunit, '(a)', iostat=status) header
+    if (status /= 0) error stop 'FAIL: join normalization header is absent'
+    read(iunit, *, iostat=status) sequence, direction, column, old_tags, &
+        new_tags, npass, factor
+    close(iunit)
+    if (status /= 0 .or. sequence /= 1 .or. direction /= 'p' &
+        .or. column /= 2 .or. any(old_tags /= [12, 13]) &
+        .or. any(new_tags /= [13, 13]) .or. any(npass /= [24, 24]) &
+        .or. factor /= 1.0d0) &
+        error stop 'FAIL: join normalization record is wrong'
 end program test_join_failure_diagnostic