spl, Fortran - rextester

Language: Layout:

program LabSpline
        integer :: n
        real :: res, arg
        real, dimension(7) :: x, y, b, c, d

n = 7
        h = 2/25.0
        x = (/-4, -2, 0, 2, 4, 6, 8/) 
        y = (/0, -2, -1, 4, 5, 2, 0/)
           
         call SPLINE(n, x, y, b, c, d)
         
         do i = 1, 151, 1
           arg = -4 + (i-1) * h
           res = SEVAL(n, arg, x, y, b, c, d)
           write(*,*) arg, res
         end do

end program LabSpline

SUBROUTINE SPLINE(N,X,Y,B,C,D)

INTEGER N

REAL X(N),Y(N),B(N),C(N),D(N)

INTEGER NM1,IB,I

REAL T

NM1=N-1

IF(N.LT.2) RETURN

IF(N.LT.3) GO TO 50

D(1)=X(2)-X(1)

C(2)=(Y(2)-Y(1))/D(1)

DO 10 I=2,NM1

D(I)=X(I+1)-X(I)

B(I)=2.*(D(I-1)+D(I))

C(I+1)=(Y(I+1)-Y(I))/D(I)

C(I)=C(I+1)-C(I)

10 CONTINUE

B(1)=-D(1)

B(N)=-D(N-1)

C(1)=0.

C(N)=0.

IF(N.EQ.3) GO TO 15

C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))

C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))

C(1)=C(1)*D(1)**2/(X(4)-X(1))

C(N)=-C(N)*D(N-1)**2/(X(N)-X(N-3))

15 DO 20 I=2,N

T=D(I-1)/B(I-1)

B(I)=B(I)-T*D(I-1)

C(I)=C(I)-T*C(I-1)

20 CONTINUE

C(N)=C(N)/B(N)

DO 30 IB=1,NM1

I=N-IB

C(I)=(C(I)-D(I)*C(I+1))/B(I)

30 CONTINUE

B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.*C(N))

DO 40 I=1,NM1

B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.*C(I))

D(I)=(C(I+1)-C(I))/D(I)

C(I)=3.*C(I)

40 CONTINUE

C(N)=3.*C(N)

D(N)=D(N-1)

RETURN

50 B(1)=(Y(2)-Y(1))/(X(2)-X(1))

C(1)=0.

D(1)=0.

B(2)=B(1)

C(2)=0.

D(2)=0.

RETURN

END

REAL FUNCTION SEVAL(N,U,X,Y,B,C,D)

INTEGER N

REAL U,X(N),Y(N),B(N),C(N),D(N)

INTEGER I,J,K

REAL DX

DATA I/1/

IF(I.GE.N) I=1

IF(U.LT.X(I)) GO TO 10

IF(U.LE.X(I+1)) GO TO 30

10 I=1

J=N+1

20 K=(I+J)/2

IF(U.LT.X(K))J=K

IF(U.GE.X(K))I=K

IF(J.GT.I+1)GO TO 20

30 DX=U-X(I)

SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))

RETURN

END

Show compiler warnings [ + ] Show input