Replacing calls to E01BAF

The SLATEC equivalent of this routine is PDA_DBINTK with order K = 4. The NAG code would look like

      INTEGER M, IFAIL
      DOUBLE PRECISION X(M), Y(M), T(M+4), C(M+4), WRK(6*M+16)
      IFAIL = 1
      CALL E01BAF( M, X, Y, T, C, M+4, WRK, 6*M+16, IFAIL )
      IF ( IFAIL .NE. 0 ) THEN
         An error has occurred
      END IF

While the size of the coefficient vector can be reduced to M for PDA_DBINTK, you now need two work spaces. The major difference is that PDA_DBINTK needs to be given the knots. So you have to calculate them in the same way as E01BAF would have done. The handling of the status is different.

      INTEGER I, K, M, IFAIL
      PARAMETER ( K = 4 )
      DOUBLE PRECISION X(M), Y(M), T(M+K), C(M)
      DOUBLE PRECISION WRK1( (2*K-1)*M ), WRK2( 2*K )
      DO 1 I = 1, K
         T(I)   = X(1)
         T(M+I) = X(M)
    1 CONTINUE
      DO 2 I = K+1, M
         T(I) = X(I-K/2)      ! Note: K is even
    2 CONTINUE
      IFAIL = 0
      CALL PDA_DBINTK( X, Y, T, M, K, C, WRK1, WRK2, IFAIL )
      IF ( IFAIL .NE. 0 ) THEN
         An error has occurred
      END IF