I ran into a limitation with the FSL module for performing N-point
Gauss-Legendre integration, gauleg.4th (FSL Algorithm #27). While the
existing algorithm is useful for numerically integrating one-dimensional
functions, it cannot be easily used for multi-dimensional integrals,
such as 2D integrals. A simple example is integrating the 2D function,

f(x, y) = x^2 + y^2

over the rectangle defined by x: 0->2 and y: 0->1. This integral may
easily be evaluated analytically to obtain the answer, 3 1/3.

The basic problem with the implementation of the integrator,

)GL-INTEGRATE

is that it is not reentrant. Thus, in a 2D integral, the instance of
)GL-INTEGRATE used to perform outer integral, say over x, cannot call
)GL-INTEGRATE to perform the inner integral over y. This is evident
because )GL-INTEGRATE stores information about the integrand function
and the arrays containing abscissas and weights for N-point
Gauss-Legendre integration in external variables:

---
: )gl-integrate ( func &x &w n -- z )

\ validate the parameter N
DUP 1 < ABORT" bad value of N (must be > 0) for )gl-integrate "

>R
& w{ &! & x{ &! is f(x)

0.0E0
R> 0 DO
x{ I } F@ f(x)
w{ I } F@ F*
F+
LOOP
; ---

This word may be rewritten to be re-entrant as follows:
---
: )gl-integrate ( func &x &w n -- z )
\ validate the parameter N
DUP 1 < ABORT" bad value of N (must be > 0) for )gl-integrate "
0.0E0
0 DO \ func &x &w ; F: rsum
over I } F@ 2 pick execute
dup I } F@ F*
F+ \ func &x &w ; F: rsum2
LOOP
2drop drop ;
---

The simple example may be evaluated using the revised )GL-INTEGRATE --
the code may be found at the link,

https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/demo/2D_integral_01.4th

The revised version of gauleg.4th, allowing for the integrator to be
reentrant, may be found at

https://github.com/mynenik/kForth-64/blob/master/forth-src/fsl/gauleg.4th

---
$ kforth64

include 2D_integral_01.4th

Type 'integral_xy fs.' to compute the 2D integral of
the function, f[x, y] = x^2 + y^2, over the rectangular
region: x: 0 -> 2, y: 0 -> 1

ok
integral_xy fs.
3.33333333333332e+00 ok
---

--
Krishna Myneni

On Saturday, March 26, 2022 at 3:46:46 AM UTC+1, Krishna Myneni wrote:
> I ran into a limitation with the FSL module for performing N-point
> Gauss-Legendre integration, gauleg.4th (FSL Algorithm #27).

I had somewhat of the reverse problem when making the NRC library
available to iForth. Some of the NRC functions, e.g. the
Davidon-Fletcher-Powell algorithm, work with function pointers.
This means C will call Forth code.

As the data transfer in the NRC is mostly done with matrices,
a Forth word in the callback must keep its state (easy: use the
stack or locals) and dynamically allocate memory for the needed
matrices. The word must be careful when calling NRC
functions when it executes, but luckily these C functions are
re-entrant already (iForth does not support the Fortran NRC :-)

There is of course a real problem when a Forth word being
callbacked calls non-re-entrant other Forth words, but this is
relatively straight-forward and something the programmer will
have to take care of himself.

-marcel

On 3/26/22 03:23, Marcel Hendrix wrote:
> On Saturday, March 26, 2022 at 3:46:46 AM UTC+1, Krishna Myneni wrote:
>> I ran into a limitation with the FSL module for performing N-point
>> Gauss-Legendre integration, gauleg.4th (FSL Algorithm #27).
>
> I had somewhat of the reverse problem when making the NRC library
> available to iForth. Some of the NRC functions, e.g. the
> Davidon-Fletcher-Powell algorithm, work with function pointers.
> This means C will call Forth code.
>

Do you provide a compiled C library of the NRC routines?

I've spent a lot of time in the past interfacing Forth (kForth-32) to
external C and Fortran libraries. For most of my needs these days, it's
better off to work with Forth source libraries, and, when efficiency is
required, a mix of Forth source and assembly language. If it's too
daunting a task to write a Forth source library, I will work in another
language. But, having a Forth source library makes the code immensely
more valuable for study (interactively from the Forth environment) and
to modify to suit my particular computing need.

> As the data transfer in the NRC is mostly done with matrices,
> a Forth word in the callback must keep its state (easy: use the
> stack or locals) and dynamically allocate memory for the needed
> matrices. The word must be careful when calling NRC
> functions when it executes, but luckily these C functions are
> re-entrant already (iForth does not support the Fortran NRC :-)
>
> There is of course a real problem when a Forth word being
> callbacked calls non-re-entrant other Forth words, but this is
> relatively straight-forward and something the programmer will
> have to take care of himself.
>

I don't think the FSL algorithms have undergone much scrutiny for
reentrancy problems, in cases where there is a possibility of occurring.
I expect that such cases are unusual. It is not likely, for example,
that a matrix inversion algorithm needs to be reentrant, but there may
be cases for block diagonal matrices.

--
Krishna