cubic interpolation subroutine 2

[starlite79]

I'd like to use IMPLICIT NONE for a subroutine called spline.f attached below. However, trying to declare all variables correctly is giving me erroneous results. The arrays are being passed by another routine that is calling spline.f. If I try to use IMPLICIT NONE as the only change, I get the obvious NO IMPLICIT type errors as well as a type/rank mismatch for argument U1.

Here is the code:

      subroutine spline(nn,x,u,s,a)
c
c$$$$ calls no other routines
c  finds array  s  for spline interpolator  eval.
c  nn  number of data points supplied (may be negative, see below)
c  x  array containing x-coordinates where function is sampled.  xx(1),xx(2),...
c     must be a strictly increasing sequence.
c  u  array containing sample values that are to be interpolated.
c  s  output array of 2nd derivative at sample points.
c  a  working space array of dimension at least  nn.
c  if the user wishes to force the derivatives at the ends of the series to
c  assume specified values, he should put du(1)/dx and du(n)/dx in s1,s2
c  and call the routine with nn=-number of terms in series.  normally a parabola
c  is fitted through the 1st and last 3 points to find the slopes.
c  if less than 4 points are supplied, straight lines are fitted.
c
      integer i,j,n,n1,nn
      real a,c,q,q1,qn,s,u,x
      dimension x(*),u(*),s(*),a(*)
c
      q(u1,x1,u2,x2)=(u1/x1**2-u2/x2**2)/(1.0/x1-1.0/x2)
c
      n=iabs(nn)
      if (n.le.3) then
c  series too short for cubic spline - use straight lines.
         do i=1,n
            s(i)=0.0
         enddo
         return
      endif
      q1=q(u(2)-u(1),x(2)-x(1),u(3)-u(1),x(3)-x(1))
      qn=q(u(n-1)-u(n),x(n-1)-x(n),u(n-2)-u(n),x(n-2)-x(n))
      if (nn.le.0) then
         q1=s(1)
         qn=s(2)
      endif
      s(1)=6.0*((u(2)-u(1))/(x(2)-x(1)) - q1)
      n1= n - 1
      do i=2,n1
         s(i)= (u(i-1)/(x(i)-x(i-1)) - u(i)*(1.0/(x(i)-x(i-1))+
     +    1.0/(x(i+1)-x(i))) + u(i+1)/(x(i+1)-x(i)))*6.0
      enddo
      s(n)=6.0*(qn + (u(n1)-u(n))/(x(n)-x(n1)))
      a(1)=2.0*(x(2)-x(1))
      a(2)=1.5*(x(2)-x(1)) + 2.0*(x(3)-x(2))
      s(2)=s(2) - 0.5*s(1)
      do i=3,n1
         c=(x(i)-x(i-1))/a(i-1)
         a(i)=2.0*(x(i+1)-x(i-1)) - c*(x(i)-x(i-1))
         s(i)=s(i) - c*s(i-1)
      enddo
      c=(x(n)-x(n1))/a(n1)
      a(n)=(2.0-c)*(x(n)-x(n1))
      s(n)=s(n) - c*s(n1)
c  back substitiute
      s(n)= s(n)/a(n)
      do j=1,n1
         i=n-j
         s(i) =(s(i) - (x(i+1)-x(i))*s(i+1))/a(i)
      enddo
      return
      end

I thought U1, U2, X1, and X2 are real and the array Q should be declared like this: Q(*,*,*,*), but this doesn't work (output from main program yields NaN). Could someone offer a hint about how to use IMPLICIT NONE correctly?

 

[xwb]

The problem with multi-dimensional arrays is that only the first dimension is variable. The compiler can't work out what to do. For instance

real x(8,*)
x(2, 7) is the ((2-1) * 8 + 7 - 1)th element

real y(4,8,*)
y(4,3,2) is the (((4 - 1) * 8 + 3 - 1) + 2 - 1)th element

If the second dimension is *, then it is not possible to work out which element it is.

 

[starlite79]

Hi,

I'm not sure I understand you completely, but I did the following to the code:

      IMPLICIT NONE
*      INTEGER I,J,N,N1,NN
      INTEGER I,J,N,N1,NN,NMAX
*      REAL A,C,Q,Q1,QN,S,U,X
      PARAMETER (NMAX = 20)
      REAL A,C,Q,Q1,QN,S,U,X,U1,U2,X1,X2
*      DIMENSION X(*),U(*),S(*),A(*)
      DIMENSION X(NMAX),U(NMAX),S(NMAX),A(NMAX)

      U1 = U(1)
      X1 = X(1)
      U2 = U(2)
      X2 = X(2)

      Q(U1,X1,U2,X2)=(U1/X1**2-U2/X2**2)/(1.0/X1-1.0/X2)
      

      N = IABS(NN)

      IF (N.LE.3) THEN

*  series too short for cubic spline - use straight lines.
         DO I=1,N
            S(I)=0.0
         ENDDO
         RETURN
      ENDIF

      Q1=Q(U(2)-U(1),X(2)-X(1),U(3)-U(1),X(3)-X(1))
      QN=Q(U(N-1)-U(N),X(N-1)-X(N),U(N-2)-U(N),X(N-2)-X(N))

      IF (NN.LE.0) THEN
         Q1=S(1)
         QN=S(2)
      ENDIF

      S(1)=6.0*((U(2)-U(1))/(X(2)-X(1)) - Q1)
      N1= N - 1

      DO I=2,N1
         S(I)= (U(I-1)/(X(I)-X(I-1)) - U(I)*(1.0/(X(I)-X(I-1))
     .   + 1.0/(X(I+1)-X(I))) + U(I+1)/(X(I+1)-X(I)))*6.0
      ENDDO

      S(N)=6.0*(QN + (U(N1)-U(N))/(X(N)-X(N1)))
      A(1)=2.0*(X(2)-X(1))
      A(2)=1.5*(X(2)-X(1)) + 2.0*(X(3)-X(2))
      S(2)=S(2) - 0.5*S(1)

      DO I=3,N1
         C=(X(I)-X(I-1))/A(I-1)
         A(I)=2.0*(X(I+1)-X(I-1)) - C*(X(I)-X(I-1))
         S(I)=S(I) - C*S(I-1)
      ENDDO

      C=(X(N)-X(N1))/A(N1)
      A(N)=(2.0-C)*(X(N)-X(N1))
      S(N)=S(N) - C*S(N1)

*  Back substitute
      S(N)= S(N)/A(N)

      DO J=1,N1
         I=N-J
         S(I) =(S(I) - (X(I+1)-X(I))*S(I+1))/A(I)
      ENDDO
      RETURN

I used the parameter NMAX=20 because in the subroutine that calls spline, all of the arrays (which appear to be only one dimensional arrays) have at most 20 constituents. Now I get the following error when I try to compile:

In file SPLINE.F:84

Q(U1,X1,U2,X2)=(U1/X1**2-U2/X2**2)/(1.0/X1-1.0/X2)
1
Error: Unexpected STATEMENT FUNCTION statement at (1)
make: *** [SPLINE.o] Error 1

I don't understand how Q can be a four dimensional array. If it is, does it need to be declared as Q(NMAX, NMAX, NMAX, NMAX)? What puzzles me the most is how, when the routine is not using IMPLICIT NONE, it compiles and gives what are possibly legitimate results...

Thanks again in advance for assistance.

 

[mikrom]

I don't understand how Q can be a four dimensional array.

Q is not an 4D array. It's STATEMENT FUNCTION.

 

[mikrom]

Statement function should appear after variables declaration and before all executable statements in the program unit - see:

http://www.obliquity.com/computer/fortran/function.html

 

[starlite79]

OK, thanks to mikrom for clarifying what Q really is. The compiler thinks this statement function is unexpected. Have my dummy arguments been declared correctly? I don't want anything to be implicitly defined.

 

[starlite79]

SOLVED:

I removed the statements,
* U1 = U(1)
* X1 = X(1)
* U2 = U(2)
* X2 = X(2)
and now it compiles and produces reasonable results. Thanks for explaining that what I was looking at was not a multi-dimensional array but rather a statement function.