ftnlf

Check-in [e62b3cf44f]
Login

Many hyperlinks are disabled.
Use anonymous login to enable hyperlinks.

Overview
Comment:Add polynomial fit/evaluation module.
Downloads: Tarball | ZIP archive | SQL archive
Timelines: family | ancestors | trunk
Files: files | file ages | folders
SHA1:e62b3cf44f4746e4217d656d9ca99f52ed0c66f8
User & Date: vadim 2019-06-16 18:22:12
Context
2019-06-16
18:22
Add polynomial fit/evaluation module. Leaf check-in: e62b3cf44f user: vadim tags: trunk
17:16
Move interpolation objects to separate module. check-in: 8f85de2af5 user: vadim tags: trunk
Changes
Hide Diffs Unified Diffs Ignore Whitespace Patch

Changes to CMakeLists.txt.

54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81






82
83
84
85

86





87
88
89
90
91
92
93
...
121
122
123
124
125
126
127










128
129
130
131
132
133
134
135
136
137
138
add_custom_command(
  OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/l2f_testmodule.f90
  COMMAND lua l2f.lua testmod.lua ${CMAKE_CURRENT_BINARY_DIR}/l2f_testmodule.f90 testmodule
  DEPENDS l2f.lua testmod.lua
  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
  )

set (src_lib
  ftnlf.f90
  ftnlf_fxcore.f90
  ftnlf_fxinterp.f90
  luaf/luaf.f90
  luaf/luafe.f90
  ${CMAKE_CURRENT_BINARY_DIR}/luaf_conf.fi
  )

set (src_lib_core
  ftnlf_fxcore.f90
  luaf/luaf.f90
  luaf/luafe.f90
  ${CMAKE_CURRENT_BINARY_DIR}/luaf_conf.fi
  # export everything
  # ftnlf_fxcore.def
  )

set (src_lib_interp
  ftnlf_fxinterp.f90
  )







set (src_lib_fx
  ${src_lib_core}
  ${src_lib_interp}

)






set (src_test
  test.f90
  ${CMAKE_CURRENT_BINARY_DIR}/l2f_testmodule.f90
  )

add_library(ftnlf STATIC ${src_lib})
................................................................................
target_link_libraries(interp Core)
set_target_properties(interp PROPERTIES COMPILE_FLAGS "${FCFLAGS90}")
set_target_properties(interp PROPERTIES LINK_FLAGS "${LFLAGS_FX}")
set_target_properties(interp PROPERTIES PREFIX "")
if(CMAKE_BUILD_TYPE STREQUAL "Release" AND CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
  set_target_properties(interp PROPERTIES INTERPROCEDURAL_OPTIMIZATION 1)
endif()











add_executable(ftnlf_test ${src_test})
# Link static
target_link_libraries(ftnlf_test ftnlf ${LUA_STATIC_LIB} ${CMAKE_DL_LIBS})
# Link dynamic
# target_link_libraries(ftnlf ${LUA_LIBRARIES})
set_target_properties(ftnlf_test PROPERTIES COMPILE_FLAGS "${FCFLAGS90}")
set_target_properties(ftnlf_test PROPERTIES LINK_FLAGS "${LFLAGS}")

# Ad-hoc
# set(CMAKE_SHARED_LIBRARY_LINK_Fortran_FLAGS "")







<
<
<
<
<
<
<
<
<












>
>
>
>
>
>




>

>
>
>
>
>







 







>
>
>
>
>
>
>
>
>
>











54
55
56
57
58
59
60









61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
...
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
add_custom_command(
  OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/l2f_testmodule.f90
  COMMAND lua l2f.lua testmod.lua ${CMAKE_CURRENT_BINARY_DIR}/l2f_testmodule.f90 testmodule
  DEPENDS l2f.lua testmod.lua
  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
  )










set (src_lib_core
  ftnlf_fxcore.f90
  luaf/luaf.f90
  luaf/luafe.f90
  ${CMAKE_CURRENT_BINARY_DIR}/luaf_conf.fi
  # export everything
  # ftnlf_fxcore.def
  )

set (src_lib_interp
  ftnlf_fxinterp.f90
  )

set (src_lib_polyfit
  xpolft.f90
  dpolft.f
  dp1vlu.f
  )

set (src_lib_fx
  ${src_lib_core}
  ${src_lib_interp}
  ${src_lib_polyfit}
)

set (src_lib
  ftnlf.f90
  ${src_lib_fx}
  )

set (src_test
  test.f90
  ${CMAKE_CURRENT_BINARY_DIR}/l2f_testmodule.f90
  )

add_library(ftnlf STATIC ${src_lib})
................................................................................
target_link_libraries(interp Core)
set_target_properties(interp PROPERTIES COMPILE_FLAGS "${FCFLAGS90}")
set_target_properties(interp PROPERTIES LINK_FLAGS "${LFLAGS_FX}")
set_target_properties(interp PROPERTIES PREFIX "")
if(CMAKE_BUILD_TYPE STREQUAL "Release" AND CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
  set_target_properties(interp PROPERTIES INTERPROCEDURAL_OPTIMIZATION 1)
endif()

# FX.polyfit library
add_library(polyfit MODULE ${src_lib_polyfit})
target_link_libraries(polyfit Core)
set_target_properties(polyfit PROPERTIES COMPILE_FLAGS "${FCFLAGS90}")
set_target_properties(polyfit PROPERTIES LINK_FLAGS "${LFLAGS_FX}")
set_target_properties(polyfit PROPERTIES PREFIX "")
if(CMAKE_BUILD_TYPE STREQUAL "Release" AND CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
  set_target_properties(polyfit PROPERTIES INTERPROCEDURAL_OPTIMIZATION 1)
endif()

add_executable(ftnlf_test ${src_test})
# Link static
target_link_libraries(ftnlf_test ftnlf ${LUA_STATIC_LIB} ${CMAKE_DL_LIBS})
# Link dynamic
# target_link_libraries(ftnlf ${LUA_LIBRARIES})
set_target_properties(ftnlf_test PROPERTIES COMPILE_FLAGS "${FCFLAGS90}")
set_target_properties(ftnlf_test PROPERTIES LINK_FLAGS "${LFLAGS}")

# Ad-hoc
# set(CMAKE_SHARED_LIBRARY_LINK_Fortran_FLAGS "")

Added dp1vlu.f.















































































































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
*DECK DP1VLU
      SUBROUTINE DP1VLU (L, NDER, X, YFIT, YP, A)
C***BEGIN PROLOGUE  DP1VLU
C***PURPOSE  Use the coefficients generated by DPOLFT to evaluate the
C            polynomial fit of degree L, along with the first NDER of
C            its derivatives, at a specified point.
C***LIBRARY   SLATEC
C***CATEGORY  K6
C***TYPE      DOUBLE PRECISION (PVALUE-S, DP1VLU-D)
C***KEYWORDS  CURVE FITTING, LEAST SQUARES, POLYNOMIAL APPROXIMATION
C***AUTHOR  Shampine, L. F., (SNLA)
C           Davenport, S. M., (SNLA)
C***DESCRIPTION
C
C     Abstract
C
C     The subroutine  DP1VLU  uses the coefficients generated by  DPOLFT
C     to evaluate the polynomial fit of degree  L , along with the first
C     NDER  of its derivatives, at a specified point.  Computationally
C     stable recurrence relations are used to perform this task.
C
C     The parameters for  DP1VLU  are
C
C     Input -- ALL TYPE REAL variables are DOUBLE PRECISION
C         L -      the degree of polynomial to be evaluated.  L  may be
C                  any non-negative integer which is less than or equal
C                  to  NDEG , the highest degree polynomial provided
C                  by  DPOLFT .
C         NDER -   the number of derivatives to be evaluated.  NDER
C                  may be 0 or any positive value.  If NDER is less
C                  than 0, it will be treated as 0.
C         X -      the argument at which the polynomial and its
C                  derivatives are to be evaluated.
C         A -      work and output array containing values from last
C                  call to  DPOLFT .
C
C     Output -- ALL TYPE REAL variables are DOUBLE PRECISION
C         YFIT -   value of the fitting polynomial of degree  L  at  X
C         YP -     array containing the first through  NDER  derivatives
C                  of the polynomial of degree  L .  YP  must be
C                  dimensioned at least  NDER  in the calling program.
C
C***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston,
C                 Curve fitting by polynomials in one variable, Report
C                 SLA-74-0270, Sandia Laboratories, June 1974.
C***ROUTINES CALLED  
C***REVISION HISTORY  (YYMMDD)
C   740601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890911  Removed unnecessary intrinsics.  (WRB)
C   891006  Cosmetic changes to prologue.  (WRB)
C   891006  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900510  Convert XERRWV calls to XERMSG calls.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DP1VLU
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER I,IC,ILO,IN,INP1,IUP,K1,K1I,K2,K3,K3P1,K3PN,K4,K4P1,K4PN,
     * KC,L,LM1,LP1,MAXORD,N,NDER,NDO,NDP1,NORD
      DOUBLE PRECISION A(*),CC,DIF,VAL,X,YFIT,YP(*)
      CHARACTER*8 XERN1, XERN2
C***FIRST EXECUTABLE STATEMENT  DP1VLU
      IF (L .LT. 0) GO TO 12
      NDO = MAX(NDER,0)
      NDO = MIN(NDO,L)
      MAXORD = A(1) + 0.5D0
      K1 = MAXORD + 1
      K2 = K1 + MAXORD
      K3 = K2 + MAXORD + 2
      NORD = A(K3) + 0.5D0
      IF (L .GT. NORD) GO TO 11
      K4 = K3 + L + 1
      IF (NDER .LT. 1) GO TO 2
      DO 1 I = 1,NDER
 1      YP(I) = 0.0D0
 2    IF (L .GE. 2) GO TO 4
      IF (L .EQ. 1) GO TO 3
C
C L IS 0
C
      VAL = A(K2+1)
      GO TO 10
C
C L IS 1
C
 3    CC = A(K2+2)
      VAL = A(K2+1) + (X-A(2))*CC
      IF (NDER .GE. 1) YP(1) = CC
      GO TO 10
C
C L IS GREATER THAN 1
C
 4    NDP1 = NDO + 1
      K3P1 = K3 + 1
      K4P1 = K4 + 1
      LP1 = L + 1
      LM1 = L - 1
      ILO = K3 + 3
      IUP = K4 + NDP1
      DO 5 I = ILO,IUP
 5      A(I) = 0.0D0
      DIF = X - A(LP1)
      KC = K2 + LP1
      A(K4P1) = A(KC)
      A(K3P1) = A(KC-1) + DIF*A(K4P1)
      A(K3+2) = A(K4P1)
C
C EVALUATE RECURRENCE RELATIONS FOR FUNCTION VALUE AND DERIVATIVES
C
      DO 9 I = 1,LM1
        IN = L - I
        INP1 = IN + 1
        K1I = K1 + INP1
        IC = K2 + IN
        DIF = X - A(INP1)
        VAL = A(IC) + DIF*A(K3P1) - A(K1I)*A(K4P1)
        IF (NDO .LE. 0) GO TO 8
        DO 6 N = 1,NDO
          K3PN = K3P1 + N
          K4PN = K4P1 + N
 6        YP(N) = DIF*A(K3PN) + N*A(K3PN-1) - A(K1I)*A(K4PN)
C
C SAVE VALUES NEEDED FOR NEXT EVALUATION OF RECURRENCE RELATIONS
C
        DO 7 N = 1,NDO
          K3PN = K3P1 + N
          K4PN = K4P1 + N
          A(K4PN) = A(K3PN)
 7        A(K3PN) = YP(N)
 8      A(K4P1) = A(K3P1)
 9      A(K3P1) = VAL
C
C NORMAL RETURN OR ABORT DUE TO ERROR
C
 10   YFIT = VAL
      RETURN
C
   11 WRITE (XERN1, '(I8)') L
      WRITE (XERN2, '(I8)') NORD
      write (0,*) 'SLATEC.DP1VLU:' //
     *   'THE ORDER OF POLYNOMIAL EVALUATION, L = '// XERN1 //
     *   ' REQUESTED EXCEEDS THE HIGHEST ORDER FIT, NORD = ' // XERN2 //
     *   ', COMPUTED BY DPOLFT -- EXECUTION TERMINATED.'
      RETURN
C
   12 write (0,*) 'SLATEC.DP1VLU:' //
     +   'INVALID INPUT PARAMETER.  ORDER OF POLYNOMIAL EVALUATION ' //
     +   'REQUESTED IS NEGATIVE.'
      RETURN
      END

Added dpolft.f.









































































































































































































































































































































































































































































































































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
*DECK DPOLFT
      SUBROUTINE DPOLFT (N, X, Y, W, MAXDEG, NDEG, EPS, R, IERR, A)
C***BEGIN PROLOGUE  DPOLFT
C***PURPOSE  Fit discrete data in a least squares sense by polynomials
C            in one variable.
C***LIBRARY   SLATEC
C***CATEGORY  K1A1A2
C***TYPE      DOUBLE PRECISION (POLFIT-S, DPOLFT-D)
C***KEYWORDS  CURVE FITTING, DATA FITTING, LEAST SQUARES, POLYNOMIAL FIT
C***AUTHOR  Shampine, L. F., (SNLA)
C           Davenport, S. M., (SNLA)
C           Huddleston, R. E., (SNLL)
C***DESCRIPTION
C
C     Abstract
C
C     Given a collection of points X(I) and a set of values Y(I) which
C     correspond to some function or measurement at each of the X(I),
C     subroutine  DPOLFT  computes the weighted least-squares polynomial
C     fits of all degrees up to some degree either specified by the user
C     or determined by the routine.  The fits thus obtained are in
C     orthogonal polynomial form.  Subroutine  DP1VLU  may then be
C     called to evaluate the fitted polynomials and any of their
C     derivatives at any point.  The subroutine  DPCOEF  may be used to
C     express the polynomial fits as powers of (X-C) for any specified
C     point C.
C
C     The parameters for  DPOLFT  are
C
C     Input -- All TYPE REAL variables are DOUBLE PRECISION
C         N -      the number of data points.  The arrays X, Y and W
C                  must be dimensioned at least  N  (N .GE. 1).
C         X -      array of values of the independent variable.  These
C                  values may appear in any order and need not all be
C                  distinct.
C         Y -      array of corresponding function values.
C         W -      array of positive values to be used as weights.  If
C                  W(1) is negative,  DPOLFT  will set all the weights
C                  to 1.0, which means unweighted least squares error
C                  will be minimized.  To minimize relative error, the
C                  user should set the weights to:  W(I) = 1.0/Y(I)**2,
C                  I = 1,...,N .
C         MAXDEG - maximum degree to be allowed for polynomial fit.
C                  MAXDEG  may be any non-negative integer less than  N.
C                  Note -- MAXDEG  cannot be equal to  N-1  when a
C                  statistical test is to be used for degree selection,
C                  i.e., when input value of  EPS  is negative.
C         EPS -    specifies the criterion to be used in determining
C                  the degree of fit to be computed.
C                  (1)  If  EPS  is input negative,  DPOLFT  chooses the
C                       degree based on a statistical F test of
C                       significance.  One of three possible
C                       significance levels will be used:  .01, .05 or
C                       .10.  If  EPS=-1.0 , the routine will
C                       automatically select one of these levels based
C                       on the number of data points and the maximum
C                       degree to be considered.  If  EPS  is input as
C                       -.01, -.05, or -.10, a significance level of
C                       .01, .05, or .10, respectively, will be used.
C                  (2)  If  EPS  is set to 0.,  DPOLFT  computes the
C                       polynomials of degrees 0 through  MAXDEG .
C                  (3)  If  EPS  is input positive,  EPS  is the RMS
C                       error tolerance which must be satisfied by the
C                       fitted polynomial.  DPOLFT  will increase the
C                       degree of fit until this criterion is met or
C                       until the maximum degree is reached.
C
C     Output -- All TYPE REAL variables are DOUBLE PRECISION
C         NDEG -   degree of the highest degree fit computed.
C         EPS -    RMS error of the polynomial of degree  NDEG .
C         R -      vector of dimension at least NDEG containing values
C                  of the fit of degree  NDEG  at each of the  X(I) .
C                  Except when the statistical test is used, these
C                  values are more accurate than results from subroutine
C                  DP1VLU  normally are.
C         IERR -   error flag with the following possible values.
C             1 -- indicates normal execution, i.e., either
C                  (1)  the input value of  EPS  was negative, and the
C                       computed polynomial fit of degree  NDEG
C                       satisfies the specified F test, or
C                  (2)  the input value of  EPS  was 0., and the fits of
C                       all degrees up to  MAXDEG  are complete, or
C                  (3)  the input value of  EPS  was positive, and the
C                       polynomial of degree  NDEG  satisfies the RMS
C                       error requirement.
C             2 -- invalid input parameter.  At least one of the input
C                  parameters has an illegal value and must be corrected
C                  before  DPOLFT  can proceed.  Valid input results
C                  when the following restrictions are observed
C                       N .GE. 1
C                       0 .LE. MAXDEG .LE. N-1  for  EPS .GE. 0.
C                       0 .LE. MAXDEG .LE. N-2  for  EPS .LT. 0.
C                       W(1)=-1.0  or  W(I) .GT. 0., I=1,...,N .
C             3 -- cannot satisfy the RMS error requirement with a
C                  polynomial of degree no greater than  MAXDEG .  Best
C                  fit found is of degree  MAXDEG .
C             4 -- cannot satisfy the test for significance using
C                  current value of  MAXDEG .  Statistically, the
C                  best fit found is of order  NORD .  (In this case,
C                  NDEG will have one of the values:  MAXDEG-2,
C                  MAXDEG-1, or MAXDEG).  Using a higher value of
C                  MAXDEG  may result in passing the test.
C         A -      work and output array having at least 3N+3MAXDEG+3
C                  locations
C
C     Note - DPOLFT  calculates all fits of degrees up to and including
C            NDEG .  Any or all of these fits can be evaluated or
C            expressed as powers of (X-C) using  DP1VLU  and  DPCOEF
C            after just one call to  DPOLFT .
C
C***REFERENCES  L. F. Shampine, S. M. Davenport and R. E. Huddleston,
C                 Curve fitting by polynomials in one variable, Report
C                 SLA-74-0270, Sandia Laboratories, June 1974.
C***ROUTINES CALLED  DP1VLU
C***REVISION HISTORY  (YYMMDD)
C   740601  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891006  Cosmetic changes to prologue.  (WRB)
C   891006  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900911  Added variable YP to DOUBLE PRECISION declaration.  (WRB)
C   920501  Reformatted the REFERENCES section.  (WRB)
C   920527  Corrected erroneous statements in DESCRIPTION.  (WRB)
C***END PROLOGUE  DPOLFT
      INTEGER I,IDEGF,IERR,J,JP1,JPAS,K1,K1PJ,K2,K2PJ,K3,K3PI,K4,
     * K4PI,K5,K5PI,KSIG,M,MAXDEG,MOP1,N,NDEG,NDER,NFAIL
      DOUBLE PRECISION TEMD1,TEMD2
      DOUBLE PRECISION A(*),DEGF,DEN,EPS,ETST,F,FCRIT,R(*),SIG,SIGJ,
     * SIGJM1,SIGPAS,TEMP,X(*),XM,Y(*),YP,W(*),W1,W11
      DOUBLE PRECISION CO(4,3)
      SAVE CO
      DATA  CO(1,1), CO(2,1), CO(3,1), CO(4,1), CO(1,2), CO(2,2),
     1      CO(3,2), CO(4,2), CO(1,3), CO(2,3), CO(3,3),
     2  CO(4,3)/-13.086850D0,-2.4648165D0,-3.3846535D0,-1.2973162D0,
     3          -3.3381146D0,-1.7812271D0,-3.2578406D0,-1.6589279D0,
     4          -1.6282703D0,-1.3152745D0,-3.2640179D0,-1.9829776D0/
C***FIRST EXECUTABLE STATEMENT  DPOLFT
      M = ABS(N)
      IF (M .EQ. 0) GO TO 30
      IF (MAXDEG .LT. 0) GO TO 30
      A(1) = MAXDEG
      MOP1 = MAXDEG + 1
      IF (M .LT. MOP1) GO TO 30
      IF (EPS .LT. 0.0D0 .AND.  M .EQ. MOP1) GO TO 30
      XM = M
      ETST = EPS*EPS*XM
      IF (W(1) .LT. 0.0D0) GO TO 2
      DO 1 I = 1,M
        IF (W(I) .LE. 0.0D0) GO TO 30
 1      CONTINUE
      GO TO 4
 2    DO 3 I = 1,M
 3      W(I) = 1.0D0
 4    IF (EPS .GE. 0.0D0) GO TO 8
C
C DETERMINE SIGNIFICANCE LEVEL INDEX TO BE USED IN STATISTICAL TEST FOR
C CHOOSING DEGREE OF POLYNOMIAL FIT
C
      IF (EPS .GT. (-.55D0)) GO TO 5
      IDEGF = M - MAXDEG - 1
      KSIG = 1
      IF (IDEGF .LT. 10) KSIG = 2
      IF (IDEGF .LT. 5) KSIG = 3
      GO TO 8
 5    KSIG = 1
      IF (EPS .LT. (-.03D0)) KSIG = 2
      IF (EPS .LT. (-.07D0)) KSIG = 3
C
C INITIALIZE INDEXES AND COEFFICIENTS FOR FITTING
C
 8    K1 = MAXDEG + 1
      K2 = K1 + MAXDEG
      K3 = K2 + MAXDEG + 2
      K4 = K3 + M
      K5 = K4 + M
      DO 9 I = 2,K4
 9      A(I) = 0.0D0
      W11 = 0.0D0
      IF (N .LT. 0) GO TO 11
C
C UNCONSTRAINED CASE
C
      DO 10 I = 1,M
        K4PI = K4 + I
        A(K4PI) = 1.0D0
 10     W11 = W11 + W(I)
      GO TO 13
C
C CONSTRAINED CASE
C
 11   DO 12 I = 1,M
        K4PI = K4 + I
 12     W11 = W11 + W(I)*A(K4PI)**2
C
C COMPUTE FIT OF DEGREE ZERO
C
 13   TEMD1 = 0.0D0
      DO 14 I = 1,M
        K4PI = K4 + I
        TEMD1 = TEMD1 + W(I)*Y(I)*A(K4PI)
 14     CONTINUE
      TEMD1 = TEMD1/W11
      A(K2+1) = TEMD1
      SIGJ = 0.0D0
      DO 15 I = 1,M
        K4PI = K4 + I
        K5PI = K5 + I
        TEMD2 = TEMD1*A(K4PI)
        R(I) = TEMD2
        A(K5PI) = TEMD2 - R(I)
 15     SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
      J = 0
C
C SEE IF POLYNOMIAL OF DEGREE 0 SATISFIES THE DEGREE SELECTION CRITERION
C
      IF (EPS) 24,26,27
C
C INCREMENT DEGREE
C
 16   J = J + 1
      JP1 = J + 1
      K1PJ = K1 + J
      K2PJ = K2 + J
      SIGJM1 = SIGJ
C
C COMPUTE NEW B COEFFICIENT EXCEPT WHEN J = 1
C
      IF (J .GT. 1) A(K1PJ) = W11/W1
C
C COMPUTE NEW A COEFFICIENT
C
      TEMD1 = 0.0D0
      DO 18 I = 1,M
        K4PI = K4 + I
        TEMD2 = A(K4PI)
        TEMD1 = TEMD1 + X(I)*W(I)*TEMD2*TEMD2
 18     CONTINUE
      A(JP1) = TEMD1/W11
C
C EVALUATE ORTHOGONAL POLYNOMIAL AT DATA POINTS
C
      W1 = W11
      W11 = 0.0D0
      DO 19 I = 1,M
        K3PI = K3 + I
        K4PI = K4 + I
        TEMP = A(K3PI)
        A(K3PI) = A(K4PI)
        A(K4PI) = (X(I)-A(JP1))*A(K3PI) - A(K1PJ)*TEMP
 19     W11 = W11 + W(I)*A(K4PI)**2
C
C GET NEW ORTHOGONAL POLYNOMIAL COEFFICIENT USING PARTIAL DOUBLE
C PRECISION
C
      TEMD1 = 0.0D0
      DO 20 I = 1,M
        K4PI = K4 + I
        K5PI = K5 + I
        TEMD2 = W(I)*((Y(I)-R(I))-A(K5PI))*A(K4PI)
 20     TEMD1 = TEMD1 + TEMD2
      TEMD1 = TEMD1/W11
      A(K2PJ+1) = TEMD1
C
C UPDATE POLYNOMIAL EVALUATIONS AT EACH OF THE DATA POINTS, AND
C ACCUMULATE SUM OF SQUARES OF ERRORS.  THE POLYNOMIAL EVALUATIONS ARE
C COMPUTED AND STORED IN EXTENDED PRECISION.  FOR THE I-TH DATA POINT,
C THE MOST SIGNIFICANT BITS ARE STORED IN  R(I) , AND THE LEAST
C SIGNIFICANT BITS ARE IN  A(K5PI) .
C
      SIGJ = 0.0D0
      DO 21 I = 1,M
        K4PI = K4 + I
        K5PI = K5 + I
        TEMD2 = R(I) + A(K5PI) + TEMD1*A(K4PI)
        R(I) = TEMD2
        A(K5PI) = TEMD2 - R(I)
 21     SIGJ = SIGJ + W(I)*((Y(I)-R(I)) - A(K5PI))**2
C
C SEE IF DEGREE SELECTION CRITERION HAS BEEN SATISFIED OR IF DEGREE
C MAXDEG  HAS BEEN REACHED
C
      IF (EPS) 23,26,27
C
C COMPUTE F STATISTICS  (INPUT EPS .LT. 0.)
C
 23   IF (SIGJ .EQ. 0.0D0) GO TO 29
      DEGF = M - J - 1
      DEN = (CO(4,KSIG)*DEGF + 1.0D0)*DEGF
      FCRIT = (((CO(3,KSIG)*DEGF) + CO(2,KSIG))*DEGF + CO(1,KSIG))/DEN
      FCRIT = FCRIT*FCRIT
      F = (SIGJM1 - SIGJ)*DEGF/SIGJ
      IF (F .LT. FCRIT) GO TO 25
C
C POLYNOMIAL OF DEGREE J SATISFIES F TEST
C
 24   SIGPAS = SIGJ
      JPAS = J
      NFAIL = 0
      IF (MAXDEG .EQ. J) GO TO 32
      GO TO 16
C
C POLYNOMIAL OF DEGREE J FAILS F TEST.  IF THERE HAVE BEEN THREE
C SUCCESSIVE FAILURES, A STATISTICALLY BEST DEGREE HAS BEEN FOUND.
C
 25   NFAIL = NFAIL + 1
      IF (NFAIL .GE. 3) GO TO 29
      IF (MAXDEG .EQ. J) GO TO 32
      GO TO 16
C
C RAISE THE DEGREE IF DEGREE  MAXDEG  HAS NOT YET BEEN REACHED  (INPUT
C EPS = 0.)
C
 26   IF (MAXDEG .EQ. J) GO TO 28
      GO TO 16
C
C SEE IF RMS ERROR CRITERION IS SATISFIED  (INPUT EPS .GT. 0.)
C
 27   IF (SIGJ .LE. ETST) GO TO 28
      IF (MAXDEG .EQ. J) GO TO 31
      GO TO 16
C
C RETURNS
C
 28   IERR = 1
      NDEG = J
      SIG = SIGJ
      GO TO 33
 29   IERR = 1
      NDEG = JPAS
      SIG = SIGPAS
      GO TO 33
 30   IERR = 2
      write (0,*) 'SLATEC.DPOLFT: INVALID INPUT PARAMETER.'
      GO TO 37
 31   IERR = 3
      NDEG = MAXDEG
      SIG = SIGJ
      GO TO 33
 32   IERR = 4
      NDEG = JPAS
      SIG = SIGPAS
C
 33   A(K3) = NDEG
C
C WHEN STATISTICAL TEST HAS BEEN USED, EVALUATE THE BEST POLYNOMIAL AT
C ALL THE DATA POINTS IF  R  DOES NOT ALREADY CONTAIN THESE VALUES
C
      IF(EPS .GE. 0.0  .OR.  NDEG .EQ. MAXDEG) GO TO 36
      NDER = 0
      DO 35 I = 1,M
        CALL DP1VLU (NDEG,NDER,X(I),R(I),YP,A)
 35     CONTINUE
 36   EPS = SQRT(SIG/XM)
 37   RETURN
      END

Changes to ftnlf.f90.

12
13
14
15
16
17
18


19
20
21
22
23
24
25
26
27

28
29
30
31
32
33
34
..
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56


57
58
59
60
61
62
63

    type(c_ptr) :: L_st

    type(luaFE_FunctionEntry), allocatable :: fx_loaders(:)
    
    character(*), parameter :: FX_Core_name = 'FX.Core'
    character(*), parameter :: FX_interp_name = 'FX.interp'


    
    integer(c_int) :: ref_array ! reference to 'array' function

contains
    ! Initialization EP
    function ftnlf_init(fname, loadlist) result(r)
        use, intrinsic :: iso_c_binding, only: c_associated, c_funloc, c_null_ptr
        use ftnlf_fxcore
        use ftnlf_fxinterp

        use luaFE
        character(*), intent(in) :: fname
        type(luaFE_FunctionEntry), optional :: loadlist(:)
        logical :: r
        integer(4) :: res, nload
    
        r = .false. ! Default
................................................................................
    103     format('Failed to initialize Lua.')
            return
        end if

        ! initialize loaders
        if (present(loadlist)) then
            nload = size(loadlist)
            allocate(fx_loaders(2+nload))
            fx_loaders(3:2+nload) = loadlist
        else
            allocate(fx_loaders(2))
        end if
        ! Core
        fx_loaders(1)%name = FX_Core_name
        fx_loaders(1)%f => ldr_fx_core
        fx_loaders(2)%name = FX_interp_name
        fx_loaders(2)%f => ldr_fx_interp



        ! The 1st stage of initialization
        res = lua_cpcall(L_st, c_funloc(linit1), c_null_ptr)
        r = lcheck(res)
        if (.not. r) return
        
        ! The 2nd stage -- load database file







>
>
|








>







 







|
|

|






>
>







12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
..
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

    type(c_ptr) :: L_st

    type(luaFE_FunctionEntry), allocatable :: fx_loaders(:)
    
    character(*), parameter :: FX_Core_name = 'FX.Core'
    character(*), parameter :: FX_interp_name = 'FX.interp'
    character(*), parameter :: FX_polyfit_name = 'FX.polyfit'
    integer, parameter :: n_modules = 3
 
    integer(c_int) :: ref_array ! reference to 'array' function

contains
    ! Initialization EP
    function ftnlf_init(fname, loadlist) result(r)
        use, intrinsic :: iso_c_binding, only: c_associated, c_funloc, c_null_ptr
        use ftnlf_fxcore
        use ftnlf_fxinterp
        use xpolft
        use luaFE
        character(*), intent(in) :: fname
        type(luaFE_FunctionEntry), optional :: loadlist(:)
        logical :: r
        integer(4) :: res, nload
    
        r = .false. ! Default
................................................................................
    103     format('Failed to initialize Lua.')
            return
        end if

        ! initialize loaders
        if (present(loadlist)) then
            nload = size(loadlist)
            allocate(fx_loaders(n_modules+nload))
            fx_loaders(n_modules+1:n_modules+nload) = loadlist
        else
            allocate(fx_loaders(n_modules))
        end if
        ! Core
        fx_loaders(1)%name = FX_Core_name
        fx_loaders(1)%f => ldr_fx_core
        fx_loaders(2)%name = FX_interp_name
        fx_loaders(2)%f => ldr_fx_interp
        fx_loaders(3)%name = FX_polyfit_name
        fx_loaders(3)%f => ldr_fx_polyfit

        ! The 1st stage of initialization
        res = lua_cpcall(L_st, c_funloc(linit1), c_null_ptr)
        r = lcheck(res)
        if (.not. r) return
        
        ! The 2nd stage -- load database file

Changes to ftnlf_fxinterp.f90.

3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
...
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
module ftnlf_fxinterp
    use LuaF
    use LuaFE
    use, intrinsic :: iso_c_binding, only: c_int, c_ptr
    use ftnlf_fxcore, only: fx_fa_new
    implicit none
    private
    ! Public API: loader function, userdata ids,
    ! Fortran functions to create and copy arrays
    public ldr_fx_interp
    public mt_interp

    character(*), parameter :: mt_interp = 'ftnlf.fx.interp'

contains

................................................................................
        
        type(c_ptr) :: ud
        integer(c_size_t) :: lentbl
        real(8), pointer :: arr(:)
        real(8) :: x, y
        
        ! Check 1st arg, get length
        ! @@@ name of userdata?
        ud = luaFE_checkudSUV(L, 1, 'interp userdata')
        ! bytes to double
        lentbl = lua_objlen(L, 1)/8

        ! Check 2nd arg
        x = luaL_checknumber(L, 2)








|
<







 







<







3
4
5
6
7
8
9
10

11
12
13
14
15
16
17
...
175
176
177
178
179
180
181

182
183
184
185
186
187
188
module ftnlf_fxinterp
    use LuaF
    use LuaFE
    use, intrinsic :: iso_c_binding, only: c_int, c_ptr
    use ftnlf_fxcore, only: fx_fa_new
    implicit none
    private
    ! Public API: loader function, userdata ids

    public ldr_fx_interp
    public mt_interp

    character(*), parameter :: mt_interp = 'ftnlf.fx.interp'

contains

................................................................................
        
        type(c_ptr) :: ud
        integer(c_size_t) :: lentbl
        real(8), pointer :: arr(:)
        real(8) :: x, y
        
        ! Check 1st arg, get length

        ud = luaFE_checkudSUV(L, 1, 'interp userdata')
        ! bytes to double
        lentbl = lua_objlen(L, 1)/8

        ! Check 2nd arg
        x = luaL_checknumber(L, 2)

Changes to test.f90.

105
106
107
108
109
110
111



112
113
114
115
116
117
118
    write (*,'(L2,1X,1PE11.4)') r, fval

    r = luafun_s('v6', [-3.5d0], fval, 'dir1', 'tbl2')
    write (*,'(L2,1X,1PE11.4)') r, fval

    r = luafun_s('v7', [1.d0,-3.5d0], fval, 'dir1', 'tbl2')
    write (*,'(L2,1X,1PE11.4)') r, fval




    r = luacache_x('tblq', 'dir1', 'tbl2')
    write (*,'(L2)') r

    r = luafun_s('v6', [1.d0], fval, cache=.true., tblname='tblq')
    write (*,'(L2,1X,1PE11.4)') r, fval








>
>
>







105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
    write (*,'(L2,1X,1PE11.4)') r, fval

    r = luafun_s('v6', [-3.5d0], fval, 'dir1', 'tbl2')
    write (*,'(L2,1X,1PE11.4)') r, fval

    r = luafun_s('v7', [1.d0,-3.5d0], fval, 'dir1', 'tbl2')
    write (*,'(L2,1X,1PE11.4)') r, fval

    r = luafun_s('v8', [4.0d0], fval, 'dir1', 'tbl2')
    write (*,'(L2,1X,1PE11.4)') r, fval

    r = luacache_x('tblq', 'dir1', 'tbl2')
    write (*,'(L2)') r

    r = luafun_s('v6', [1.d0], fval, cache=.true., tblname='tblq')
    write (*,'(L2,1X,1PE11.4)') r, fval

Changes to testdb.lua.

1
2

3

4
5
6
7
8
9
10
..
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
local FXC = require('FX.Core')
local FXi = require('FX.interp')

iii = FXi.interp({0,1,10,0,15,-3})


tbl1 = { val1 = 17.4e-1,
v2 = function(x,y) return (x/y) end,
v3 = iii,
valN = function(x,y,ai,am)
   local ao = FXC.array(10)
   ao[1] = (x+y)*ai[1]
................................................................................
end,
}

dir1 = {
   {},
   { v5b = function(x,y) return x-y+1 end },
   tbl1 = { v5 = function(x,y) return x-y end },
   tbl2 = { v6 = FXi.interp({-10,1,10,70}), v7 = -29.3e-4},

}

if true then

local T = require('FX.Test')
print(T.testfun(-1.4))


>

>







 







|







1
2
3
4
5
6
7
8
9
10
11
12
..
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
local FXC = require('FX.Core')
local FXi = require('FX.interp')
local FXp = require('FX.polyfit')
iii = FXi.interp({0,1,10,0,15,-3})
local pfit = FXp.polyfit(FXC.apackt{1., 2., 3.}, FXC.apackt{0., 1., 4.}, 2, 0.1)

tbl1 = { val1 = 17.4e-1,
v2 = function(x,y) return (x/y) end,
v3 = iii,
valN = function(x,y,ai,am)
   local ao = FXC.array(10)
   ao[1] = (x+y)*ai[1]
................................................................................
end,
}

dir1 = {
   {},
   { v5b = function(x,y) return x-y+1 end },
   tbl1 = { v5 = function(x,y) return x-y end },
   tbl2 = { v6 = FXi.interp({-10,1,10,70}), v7 = -29.3e-4, v8 = pfit },

}

if true then

local T = require('FX.Test')
print(T.testfun(-1.4))

Changes to testfx.lua.

1
2

3
4
5
6
7
8
9
..
88
89
90
91
92
93
94

















local FXC = require('FX.Core')
local FXi = require('FX.interp')


local A = FXC.array(4)
print(#A)
print(A[1])
A[1] = 96.7
print(A[1])
A[1] = 0.1
................................................................................

local X1 = FXC.apack(nil,nil)
print('#[]=',#X1)
print(X1,getmetatable(X1))

local Z = FXC.array(0)
print(Z,getmetatable(Z))



















>







 







>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
..
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
local FXC = require('FX.Core')
local FXi = require('FX.interp')
local FXp = require('FX.polyfit')

local A = FXC.array(4)
print(#A)
print(A[1])
A[1] = 96.7
print(A[1])
A[1] = 0.1
................................................................................

local X1 = FXC.apack(nil,nil)
print('#[]=',#X1)
print(X1,getmetatable(X1))

local Z = FXC.array(0)
print(Z,getmetatable(Z))

print('polyfit test')
local fit, st, eps, yev = FXp.polyfit(FXC.apackt{1., 2., 3.}, FXC.apackt{3., 4.9, 7.2}, 1, 0.1)
print(st, eps, fit[#fit])

fit, st, eps, yev = FXp.polyfit(FXC.apackt{1., 2., 3.}, FXC.apackt{0., 1., 4.}, 2, 0.1)
print(st, eps, fit[#fit])

print(FXp.polyval(4, fit))
local y, yder = FXp.polyval(4, fit,0)
print(y, #yder)
y, yder = FXp.polyval(4, fit, 1)
print(y, yder[1])
print(FXp.polyval(0, fit, nil, 1))
y, yder = FXp.polyval(4, fit, 1, 1)
print(y, yder[1])
print(fit(4), fit(0))

Added xpolft.f90.

























































































































































































































































































































































































































































































































































>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
! Fortran extensions to use in Lua modules
! (polynomial fit)
module xpolft
    use LuaF
    use LuaFE
    use, intrinsic :: iso_c_binding, only: c_int, c_ptr
    use ftnlf_fxcore, only: fx_fa_new
    implicit none
    private
    public ldr_fx_polyfit
    ! Public API: loader function, userdata ids
    public mt_polyfit

    character(*), parameter :: mt_polyfit = 'ftnlf.fx.polyfit'

contains

    ! FX.polyfit module loader
    function ldr_fx_polyfit(L) bind(C, name='luaopen_FX_polyfit') result(r)
        use, intrinsic :: iso_c_binding, only: c_funloc
        type(c_ptr), value, intent(in) :: L
        integer(c_int) :: r

        ! Pop package name
        call lua_pop(L, 1)

        call lua_newtable(L)
        call luaFE_register(L, 'polyfit', l_polyfit)
        call luaFE_register(L, 'polyval', l_polyval)

        r = 1
    end function ldr_fx_polyfit

    ! calculate polynomial fit
    ! syntax
    ! fit, st, eps, yev = polyfit(x, y, maxdeg, eps)
    ! where:
    ! (input params)
    ! x and y are Fortran arrays of points with the same length
    ! maxdeg and eps are scalar numeric parameters of DPOLFT
    ! (output params)
    ! fit is output fit array (A and ndeg) or false if failed
    ! st is status (IERR)
    ! eps is output RMS error
    ! yev are fit evaluated at x points (R array)

    function l_polyfit(L) bind(C) result(r)
        use ftnlf_fxcore

        use, intrinsic :: iso_c_binding, only: c_f_pointer, c_intptr_t
        type(c_ptr), value, intent(in) :: L
        integer(c_int) :: r

        real(8) :: eps
        integer :: maxdeg, ierr, ndeg

        ! Arrays for processing
        real(8), pointer, dimension(:) :: a_x, a_y, a_a, a_r
        real(8), allocatable, dimension(:) :: a_w

        integer :: n_x, n_y, len_a
        type(c_ptr) :: ud

        ! Extract scalar vars
        maxdeg = luaL_checkint(L, 3)
        eps = luaL_checknumber(L, 4)

        ! Extract Fortran arrays
        ! Check arg x, get length
        ud = luaL_checkudata(L, 1, F_C_STR(mt_FA))
        ! bytes to double
        n_x = INT(lua_objlen(L, 1)/8, 4)
        ! Associate
        call c_f_pointer(ud, a_x, [n_x])
        ! Check arg y, get length
        ud = luaL_checkudata(L, 2, F_C_STR(mt_FA))
        ! bytes to double
        n_y = INT(lua_objlen(L, 2)/8, 4)
        call luaFE_check(L, n_x == n_y, &
            & 'x and y arrays are of different length')
        ! Associate
        call c_f_pointer(ud, a_y, [n_y])

        ! Discard extra parameters
        call lua_settop(L, 4)

        ! create A array and push it as the 5th element
        len_a = 3*n_x + 3*maxdeg + 3 + 1
        call fx_fa_new(L, len_a, mt_polyfit)
        ud = lua_touserdata(L, 5)
        ! Associate
        call c_f_pointer(ud, a_a, [len_a])
        
        ! create R array and push it as the 6th element
        call fx_fa_new(L, n_x)
        ud = lua_touserdata(L, 6)
        ! Associate
        call c_f_pointer(ud, a_r, [n_x])

        ! Allocate weight arrays
        allocate(a_w(n_x))
        ! Set weights
        a_w(1:n_x) = 1.d0

        ! calculate fit
        call DPOLFT(n_x, a_x, a_y, a_w, maxdeg, ndeg, eps, a_r, ierr, a_a)

        ! store output parameters

        if (ierr == 2) then
            ! Invalid parameters, no calculated values
            call lua_pushnil(L)

            ! store st
            call lua_pushinteger(L, INT(ierr, C_INTPTR_T))

            r = 2
        else
            ! store ndeg
            a_a(len_a) = dble(ndeg)
            ! store fit (a_a)
            call lua_pushvalue(L, 5)
            ! Get fit metatable
            r = lua_getmetatable(L, -1)
            ! Add call metamethod
            call luaFE_registerSUV(L, '__call', l_polyfit_call)
            ! Pop metatable
            call lua_pop(L, 1)

            ! store st
            call lua_pushinteger(L, INT(ierr, C_INTPTR_T))
            ! store eps
            call lua_pushnumber(L, eps)
            ! store yev (a_r)
            call lua_pushvalue(L, 6)
            r = 4
        end if
        
        ! disassociate pointers
        a_a => null()
        a_r => null()
        a_x => null()
        a_y => null()
        deallocate(a_w)

    end function l_polyfit

    ! evaluate polynomial fit
    ! syntax
    ! y, yder = polyval(x, fit, nder, ndeg)
    ! where:
    ! (input params)
    ! x is scalar input point
    ! fit is array returned by polyfit
    ! nder is the number of derivatives requested
    !   nil/none if no derivatives are requested
    ! ndeg is the degree of polynomial evaluated
    !   nil/none if degree is stored in the fit array
    ! (output params)
    ! y is evaluated value
    ! yder is derivatives array (if requested)

    function l_polyval(L) bind(C) result(r)
        use ftnlf_fxcore

        use, intrinsic :: iso_c_binding, only: c_f_pointer
        type(c_ptr), value, intent(in) :: L
        integer(c_int) :: r

        integer :: nder, ndeg
        real(8) :: x, y
        logical :: der

        ! Arrays for processing
        real(8), pointer, dimension(:) :: a_a, yder

        integer :: len_a
        type(c_ptr) :: ud

        real(8), target :: yder0(0)

        ! Extract vars
        x = luaL_checknumber(L, 1)
        ! Fit array
        ! Check arg, get length
        ud = luaL_checkudata(L, 2, F_C_STR(mt_polyfit))
        ! bytes to double
        len_a = INT(lua_objlen(L, 2)/8, 4)
        ! Associate
        call c_f_pointer(ud, a_a, [len_a])

        ! Fill absent parameters and discard extra parameters
        call lua_settop(L, 4)

        if (lua_isnil(L, 3)) then
            nder = 0
            der = .false.
        else
            nder = luaL_checkint(L, 3)
            der = .true.
        end if

        if (lua_isnil(L, 4)) then
            ndeg = int(a_a(len_a) + 0.5d0)
        else
            ndeg = luaL_checkint(L, 4)
        end if

        if (der) then
            ! create yder array and push it as the 5th element
            call fx_fa_new(L, nder)
            ud = lua_touserdata(L, 5)
            ! Associate
            call c_f_pointer(ud, yder, [nder])
        else
            yder => yder0
        end if
        
        ! calculate value
        call DP1VLU(ndeg, nder, x, y, yder, a_a)

        ! store output parameters
        call lua_pushnumber(L, y)
        r = 1
        if (der) then
            ! store yder
            call lua_pushvalue(L, 5)
            r  = 2
        end if
        
        ! disassociate pointers
        a_a => null()
        yder => null()
    
    end function l_polyval

    ! polyfit userdata metatable __call method
    function l_polyfit_call(L) bind(C) result(r)
        use, intrinsic :: iso_c_binding, only: c_funloc
        type(c_ptr), value, intent(in) :: L
        integer(c_int) :: r

        type(c_ptr) :: ud
        
        ! Check 1st arg
        ud = luaFE_checkudSUV(L, 1, 'polyfit userdata')

        ! Adjust to one input variable
        call lua_settop(L, 2)

        ! Push polyval function
        call lua_pushcfunction(L, c_funloc(l_polyval))

        ! stack is now: fit arg1 <polyval>

        call lua_insert(L, 1)
        call lua_insert(L, 2)
        ! stack is now: <polyval> arg1 fit

        ! call polyval with 2 input, 1 output parameters
        call lua_call(L, 2, 1)
        
        ! return its output

        r = 1
    end function l_polyfit_call

end module xpolft