From 643b219ed2acfbcbfc1c590ff22add750a40fb58 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 16 Jul 2025 17:28:59 +0200 Subject: [PATCH 1/5] Fast default splines --- COMMON/spline_cof.f90 | 79 ++++++++++++++++++++++++++++++++++++------- 1 file changed, 66 insertions(+), 13 deletions(-) diff --git a/COMMON/spline_cof.f90 b/COMMON/spline_cof.f90 index f51bb48e..db4e9a27 100644 --- a/COMMON/spline_cof.f90 +++ b/COMMON/spline_cof.f90 @@ -1,6 +1,6 @@ !*********************************************************************** -! +! ! routines for calculating spline coefficients ! drivers ! @@ -13,11 +13,47 @@ !*********************************************************************** -! +! ! routines for third order spline ! !*********************************************************************** +module fastspline + +contains + +SUBROUTINE splinecof3_fast(x, y, a, b, c, d) + use nrtype, only : I4B, DP + real(DP), dimension(:), intent(in) :: x, y + real(DP), dimension(:), intent(out) :: a, b, c, d + + integer(I4B) :: info + real(DP) :: r(size(x)-1), h(size(x)-1), dl(size(x)-3), ds(size(x)-2), cs(size(x)-2) + integer(I4B) :: n + + n = size(x) + + h = x(2:) - x(1:n-1) + r = y(2:) - y(1:n-1) + + dl = h(2:n-2) + ds = 2d0*(h(1:n-2)+h(2:)) + + cs = 3d0*(r(2:)/h(2:)-r(1:n-2)/h(1:n-2)) + + call dptsv(n-2, 1, ds, dl, cs, n-2, info) + a(1:n-1) = y(1:n-1) + b(1) = r(1)/h(1) - h(1)/3d0*cs(1) + b(2:n-2) = r(2:n-2)/h(2:n-2)-h(2:n-2)/3d0*(cs(2:n-2) + 2d0*cs(1:n-3)) + b(n-1) = r(n-1)/h(n-1)-h(n-1)/3d0*(2d0*cs(n-2)) + c(4) = 0 + c(2:n-1) = cs + d(1) = 1d0/(3d0*h(1))*cs(1) + d(2:n-2) = 1d0/(3d0*h(2:n-2))*(cs(2:n-2)-cs(1:n-3)) + d(n-1) = 1d0/(3d0*h(n-1))*(-cs(n-2)) +END SUBROUTINE splinecof3_fast + +end module fastspline ! ------ third order spline: with testfunction, LSQ, smoothing ! @@ -61,7 +97,7 @@ !> INTEGER(I4B), PARAMETER :: VAR = 7 ... no of variables !> !> NEEDS: -!> solve_systems, calc_opt_lambda3 +!> calc_opt_lambda3 SUBROUTINE splinecof3_a(x, y, c1, cn, lambda1, indx, sw1, sw2, & a, b, c, d, m, f) !----------------------------------------------------------------------- @@ -75,6 +111,7 @@ SUBROUTINE splinecof3_a(x, y, c1, cn, lambda1, indx, sw1, sw2, & !(Bad performance for more than 1000 flux surfaces ~ (3*nsurf)^2) USE sparse_mod, ONLY : sparse_solve !! End Modifications by Andreas F. Martitsch (06.08.2014) + use fastspline, only: splinecof3_fast !--------------------------------------------------------------------- @@ -112,7 +149,23 @@ END FUNCTION f REAL(DP), DIMENSION(:), ALLOCATABLE :: inh, simqa, lambda, omega character(200) :: error_message - len_x = SIZE(x) + if (.not. m == 0) goto 100 ! skip if m is not 0 + if (.not. sw1 == 2 .and. sw2 == 4) goto 100 ! skip if not natural boundary condis + if (.not. abs(c1 - 0d0) < 1d-13) goto 100 ! skip if nonzero boundary condi + if (.not. abs(cn - 0d0) < 1d-13) goto 100 ! skip if nonzero boundary condi + if (.not. all(abs(lambda1 - 1d0) < 1d-13)) goto 100 ! skip if lambda1 is not 1 + + call splinecof3_fast(x, y, a, b, c, d) + a(size(x)) = 0d0 + b(size(x)) = 0d0 + c(size(x)) = 0d0 + d(size(x)) = 0d0 + + return + + + ! Cannot use fast splines, fall back to long spline routine +100 len_x = SIZE(x) len_indx = SIZE(indx) size_dimension = VAR * len_indx - 2 @@ -562,15 +615,15 @@ END FUNCTION f ! solve system CALL sparse_solve(MA, inh) - + ! take a(), b(), c(), d() DO i = 1, len_indx - a(i) = inh((i-1)*VAR+1) - b(i) = inh((i-1)*VAR+2) + a(i) = inh((i-1)*VAR+1) + b(i) = inh((i-1)*VAR+2) c(i) = inh((i-1)*VAR+3) - d(i) = inh((i-1)*VAR+4) + d(i) = inh((i-1)*VAR+4) END DO - + DEALLOCATE(MA, stat = i_alloc) IF(i_alloc /= 0) STOP 'splinecof3: Deallocation for arrays 1 failed!' @@ -608,7 +661,7 @@ SUBROUTINE reconstruction3_a(ai, bi, ci, di, h, a, b, c, d) REAL(DP), INTENT(IN) :: ai, bi, ci, di REAL(DP), INTENT(IN) :: h - REAL(DP), INTENT(OUT) :: a, b, c, d + REAL(DP), INTENT(OUT) :: a, b, c, d !--------------------------------------------------------------------- @@ -687,7 +740,7 @@ END FUNCTION f REAL(DP) :: h REAL(DP), DIMENSION(:), ALLOCATABLE :: xn, yn, lambda1 REAL(DP), DIMENSION(:), ALLOCATABLE :: ai, bi, ci, di - + no = SIZE(x) ns = SIZE(a) len_indx = SIZE(indx) @@ -869,11 +922,11 @@ END FUNCTION f ! weights: w(i)=0/1; if(w(i)==0) ... do not use this point w = 1 - + sw1 = 2 sw2 = 4 - c1 = 0.0D0 + c1 = 0.0D0 cn = 0.0D0 DO i = 1, no_cur From 5e92d4f937ce90d783b2d3262f5a24c6fea0af15 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 16 Jul 2025 17:32:05 +0200 Subject: [PATCH 2/5] Fix natural boundary conditions --- COMMON/spline_cof.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/COMMON/spline_cof.f90 b/COMMON/spline_cof.f90 index db4e9a27..7f01008c 100644 --- a/COMMON/spline_cof.f90 +++ b/COMMON/spline_cof.f90 @@ -150,7 +150,7 @@ END FUNCTION f character(200) :: error_message if (.not. m == 0) goto 100 ! skip if m is not 0 - if (.not. sw1 == 2 .and. sw2 == 4) goto 100 ! skip if not natural boundary condis + if (.not. (sw1 == 2 .and. sw2 == 4)) goto 100 ! skip if not natural boundary condis if (.not. abs(c1 - 0d0) < 1d-13) goto 100 ! skip if nonzero boundary condi if (.not. abs(cn - 0d0) < 1d-13) goto 100 ! skip if nonzero boundary condi if (.not. all(abs(lambda1 - 1d0) < 1d-13)) goto 100 ! skip if lambda1 is not 1 From 45f575cb54f0b5f4029b2513315754098418e389 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 16 Jul 2025 17:36:58 +0200 Subject: [PATCH 3/5] Fix index of spline coefficient --- COMMON/spline_cof.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/COMMON/spline_cof.f90 b/COMMON/spline_cof.f90 index 7f01008c..f880e9b5 100644 --- a/COMMON/spline_cof.f90 +++ b/COMMON/spline_cof.f90 @@ -46,8 +46,8 @@ SUBROUTINE splinecof3_fast(x, y, a, b, c, d) b(1) = r(1)/h(1) - h(1)/3d0*cs(1) b(2:n-2) = r(2:n-2)/h(2:n-2)-h(2:n-2)/3d0*(cs(2:n-2) + 2d0*cs(1:n-3)) b(n-1) = r(n-1)/h(n-1)-h(n-1)/3d0*(2d0*cs(n-2)) - c(4) = 0 - c(2:n-1) = cs + c(1) = 0d0 + c(2:n-1) = cs d(1) = 1d0/(3d0*h(1))*cs(1) d(2:n-2) = 1d0/(3d0*h(2:n-2))*(cs(2:n-2)-cs(1:n-3)) d(n-1) = 1d0/(3d0*h(n-1))*(-cs(n-2)) From a600c474d7e5c490c5b1299322ca54e210698a15 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 16 Jul 2025 17:40:20 +0200 Subject: [PATCH 4/5] Add error handling for LAPACK call --- COMMON/spline_cof.f90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/COMMON/spline_cof.f90 b/COMMON/spline_cof.f90 index f880e9b5..2dc83d07 100644 --- a/COMMON/spline_cof.f90 +++ b/COMMON/spline_cof.f90 @@ -42,6 +42,17 @@ SUBROUTINE splinecof3_fast(x, y, a, b, c, d) call dptsv(n-2, 1, ds, dl, cs, n-2, info) + if (info /= 0) then + if (info < 0) then + write(*,'(A,I0,A)') 'splinecof3_fast: LAPACK dptsv error - illegal value in argument ', -info, '.' + error stop 'splinecof3_fast: Invalid argument to dptsv' + else + write(*,'(A,I0,A)') 'splinecof3_fast: LAPACK dptsv error - diagonal element ', info, ' is zero.' + write(*,*) 'The tridiagonal system is singular and cannot be solved.' + error stop 'splinecof3_fast: Singular tridiagonal system in dptsv' + end if + end if + a(1:n-1) = y(1:n-1) b(1) = r(1)/h(1) - h(1)/3d0*cs(1) b(2:n-2) = r(2:n-2)/h(2:n-2)-h(2:n-2)/3d0*(cs(2:n-2) + 2d0*cs(1:n-3)) From a587ecc1dda0f829aa7964f9af9696c36dbbe3e4 Mon Sep 17 00:00:00 2001 From: Christopher Albert Date: Wed, 16 Jul 2025 17:45:46 +0200 Subject: [PATCH 5/5] Remove trailing whitespaces --- COMMON/spline_cof.f90 | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/COMMON/spline_cof.f90 b/COMMON/spline_cof.f90 index 2dc83d07..42e19a7d 100644 --- a/COMMON/spline_cof.f90 +++ b/COMMON/spline_cof.f90 @@ -1,6 +1,6 @@ !*********************************************************************** -! +! ! routines for calculating spline coefficients ! drivers ! @@ -13,7 +13,7 @@ !*********************************************************************** -! +! ! routines for third order spline ! !*********************************************************************** @@ -626,15 +626,15 @@ END FUNCTION f ! solve system CALL sparse_solve(MA, inh) - + ! take a(), b(), c(), d() DO i = 1, len_indx - a(i) = inh((i-1)*VAR+1) - b(i) = inh((i-1)*VAR+2) + a(i) = inh((i-1)*VAR+1) + b(i) = inh((i-1)*VAR+2) c(i) = inh((i-1)*VAR+3) - d(i) = inh((i-1)*VAR+4) + d(i) = inh((i-1)*VAR+4) END DO - + DEALLOCATE(MA, stat = i_alloc) IF(i_alloc /= 0) STOP 'splinecof3: Deallocation for arrays 1 failed!' @@ -672,7 +672,7 @@ SUBROUTINE reconstruction3_a(ai, bi, ci, di, h, a, b, c, d) REAL(DP), INTENT(IN) :: ai, bi, ci, di REAL(DP), INTENT(IN) :: h - REAL(DP), INTENT(OUT) :: a, b, c, d + REAL(DP), INTENT(OUT) :: a, b, c, d !--------------------------------------------------------------------- @@ -751,7 +751,7 @@ END FUNCTION f REAL(DP) :: h REAL(DP), DIMENSION(:), ALLOCATABLE :: xn, yn, lambda1 REAL(DP), DIMENSION(:), ALLOCATABLE :: ai, bi, ci, di - + no = SIZE(x) ns = SIZE(a) len_indx = SIZE(indx) @@ -933,11 +933,11 @@ END FUNCTION f ! weights: w(i)=0/1; if(w(i)==0) ... do not use this point w = 1 - + sw1 = 2 sw2 = 4 - c1 = 0.0D0 + c1 = 0.0D0 cn = 0.0D0 DO i = 1, no_cur