gnunesjr
03-02-2007, 02:22 PM
I have a problem that is caused by NumRec's aggressive use of the efficiencies offered by Fortran 90. Maybe someone can help me unravel this...
I need to calculate some function g(z) at a large number of different z values.
At each value of z, evaluating g(z) involves a numerical integral over the function f(x). To evaluate f(x) at a particular x also requires a numerical integral.
Rather than perform the same integral thousands of times, my strategy is to calculate f(x) only once over the entire allowed range of x, and store the results in an array. Then, when qsimp needs to evaluate f(x), my function would use polynomial interpolation to quickly return the correct value.
Now we get to the problem: when qsimp calls my function f(x), x is not a simple scalar. It is an array of some size (which I guess was set by the routine trapzd in the interests of efficiency).
But the routine polint wants a scalar argument. How do I correctly call polint on all the x's and return the right kind of object to qsimp? Some F90 judo is required here, and I've only got a green belt.
I need to calculate some function g(z) at a large number of different z values.
At each value of z, evaluating g(z) involves a numerical integral over the function f(x). To evaluate f(x) at a particular x also requires a numerical integral.
Rather than perform the same integral thousands of times, my strategy is to calculate f(x) only once over the entire allowed range of x, and store the results in an array. Then, when qsimp needs to evaluate f(x), my function would use polynomial interpolation to quickly return the correct value.
Now we get to the problem: when qsimp calls my function f(x), x is not a simple scalar. It is an array of some size (which I guess was set by the routine trapzd in the interests of efficiency).
But the routine polint wants a scalar argument. How do I correctly call polint on all the x's and return the right kind of object to qsimp? Some F90 judo is required here, and I've only got a green belt.