@mecej4 Thank you! Your code is correct.
Yes we indeed prevented the use of cumsum function/subroutine explicitly. Uhm, however the calculation of cumsum is still there in the loop, and the time cost is the same if not more. But nonetheless your code/suggestion is insightful!
The reason I am curious about cumsum is from Thomas Schon’s Matlab code,
http://user.it.uu.se/~thosc112/research/rao-blackwellized-particle.html
In the pf.m, he uses his sysresample function which uses cumsum, Matlab code is below,
function i=sysresample(q)
qc=cumsum(q); M=length(q);
u=([0:M-1]+rand(1))/M;
i=zeros(1,M); k=1;
for j=1:M
while (qc(k)<u(j))
k=k+1;
end
i(j)=k;
end
end
I support putting cumsum into Fortran’s intrinsic function like the Matlab one, Cumulative sum - MATLAB cumsum
I found a link, leonfoks seems already did something about this,
PS.
Just to be clear, I do not have any intention to waste your guy’s time to solve my problem for me.
In this post I was simply just asking about cumsum function because it is simple enough.
I know cumsum not something really slow down my code because it is O(N) which is not too bad after all.
Just in case someone want to see a minimal working example, I put my not even baked MWE in the link, Files · master · Chen Rong / SDE-PF · GitLab
If you have visual studio + intel OneAPI installed, just click the .sln file and it will open the project, just build and run, it should take 0.15 sec or less.
main.f90 is main program. The cumsum is in the file pf.f90, that is the module pf (particle filter), in the sysresample subroutine. around line 153. Above it, which is commented now is my translated version with cumsum, below line 153 is the one @mecej4 suggested without explicit cumsum.
The very naive SDE solver (which is not used in this example in this thread) I slightly wrote is based on Kasdin’s paper and John Burkardt’s code
Note that John’s stochastic_RK is a scalar version which only solve 1 equation. While Kasdin’s paper provides a general solution which can solve a system of SDEs, so it is actually a vector version.
I plan to learn more from FLINT ODE code by @prince_mahajan
and make my SDE solver much better than it is now.
About SDE and Kasdin’s method, a relevant S.O post is,
The comments there by Lutz Lehmann about SDE is very useful.
I also asked Kasdin’s method at
Chris Rack told me that Kasdin’s Runge-Kutta method is like order 0.5 than order 4. But nonetheless here I will begin with Kasdin’s method.
Thing is, in Fortran, we does not seem to have some good SDE solver now it seems.