Fast spline implementation

COMMON/spline_cof.f90

@@ -17,7 +17,54 @@
 ! 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)
+
+  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))
+  b(n-1)   = r(n-1)/h(n-1)-h(n-1)/3d0*(2d0*cs(n-2))
+  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))
+END SUBROUTINE splinecof3_fast
+
+end module fastspline
 
 ! ------  third order spline: with testfunction, LSQ, smoothing
 !
@@ -61,7 +108,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 +122,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 +160,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