# What is the fastest way to do cumulative sum in Fortran to mimic Matlab cumsum? **URL:** \ **Category:** Uncategorized\ **Created:** [October 1, 2021, 2:27pm UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976 "2021-10-01T14:27:54Z")\ **Posts on this page:** 20\ **Page:** 1
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 1, 2021, 2:27pm UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/1 "2021-10-01T14:27:54Z")
Deal all, I wanted to create a cumsum subroutine which mimic the the following cumsum function in Matlab, [B = cumsum(A)](https://www.mathworks.com/help/matlab/ref/cumsum.html#d123e264509) Simple. Eg, B and A are both 1d array with the same size and shape. I wanted to create B as cumulative sum of A such that for each element i, ```auto B(i) = sum(A(:i)) ``` Below is my cumsum sloppy subroutine for illustration, ``` subroutine cumsum(A,B) real(kind=r8) :: A(:),B(size(A)) integer(kind=i8) :: i B(1) = A(1) do i = 2, size(A) B(i) = B(i-1) + A(i) enddo return end subroutine cumsum ``` I wonder, is there faster way to do cumsum what perhaps fully utilize vectorization or something? Thank you very much in advance! In response to a question at Stack Overflow, ‘King’ has provided a useful link, > ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [October 1, 2021, 5:36pm UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/4 "2021-10-01T17:36:52Z")
Your code is faulty: in the DO loop, when i = 1 the statement in the loop references qout(0), which is a bounds violation. You can fix this bug easily. Please define what you mean by “best way” and what you mean by “mimic Matlab”. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 1, 2021, 6:23pm UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/5 "2021-10-01T18:23:19Z")
Thank you. I updated my description. Thank you for pointing out the potential issue, my code is more for illustration, I am sorry it is a little sloppy. I would require A and B in the input start from index 1. I would actually wish A can be any 1D array starting from any index, does not have to start from 1. In the i loop, i should start from the lbound(A), end at lbound(A) + size(A) -1. But that is not the most relevant for now. I wish to find the ‘fastest’ way to do cumsum. ---
**Author:** ![everythingfunctional](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/everythingfunctional/32/176_2.png) [@everythingfunctional](https://fortran-lang.discourse.group/u/everythingfunctional)\ **Post date:** [October 1, 2021, 6:26pm UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/6 "2021-10-01T18:26:27Z")
```auto cumsum = [(sum(a(1:i)), i = 1, size(a))] ``` The compiler should be allowed to vectorize, parallelize, do whatever optimizations it likes with that. ---
**Author:** ![stavros](https://avatars.discourse-cdn.com/v4/letter/s/51bf81/32.png) [@stavros](https://fortran-lang.discourse.group/u/stavros)\ **Post date:** [October 2, 2021, 5:56am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/11 "2021-10-02T05:56:16Z")
Correction inside loop: ```auto do i = 2, size(A) B(i) = A(i-1) + A(i) enddo ``` Also, the most elegant way is how @everythingfunctional suggested, but with `cumsum` corrected as an array not scalar. ```auto pure function cumsum(a) result (r) real(dp), intent(in) :: a(:) real(dp) :: r(size(a)) integer :: i r(:) = [(sum(a(1:i)),i=1,size(a))] end function cumsum ``` then you simply call it as you initially intended: `b=cumsum(a)` here is a working example [link](https://godbolt.org/z/KTv84caK9) ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 2, 2021, 6:19am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/12 "2021-10-02T06:19:33Z")
Thank you so much @stavros Yes, this is truly an elegant one line solution. Indeed, let me say that I do not mean to offend anyone, nothing personal. Just from business point of view, this elegant code may not be efficient enough. The cost is O( (size(a))^2 ). The code I show is O( size(a) ) which is a factor of size(a) faster (however I wish to further improve it and make it faster.) ``` subroutine cumsum(a,b) real(kind=r8) :: a(:),b(size(a)) integer(kind=i8) :: i b(1)=a(1) do i = 2, size(a) b(i) = b(i-1) + a(i) enddo return end subroutine cumsum ``` You could check the speed difference with a array whose size(a) = 1000000. Note that I did not set the size = 1000000 just to show that the elegant code is slow. I set it as 1000000 is because I need this number of particles in my particle filter (sequential Monte Carlo). I can use a + random array to propose the move of all the particles randomly and fully using vectorization. This may be faster than do loop which loop over 1000000 particles. Slightly more detailed, in the below code, np = 1000000, ``` noise = matmul(reshape(gaussian(gbl%nx*np),[np,gbl%nx]),sq) x0tmp = spread(gbl%x0,1,np) call mf(x0tmp,gbl%kel) x = x0tmp + noise ! this is the a + random array subroutine mf(x,kel) real(kind=r8) :: x(:,:),kel x(:,1) = exp(-0.2_r8*kel)*x(:,1) return end subroutine mf subroutine sysresample(q,i) use random integer(kind=i8) :: m,j,k,l real(kind=r8) :: q(:),qc(size(q)),rn(1),u(size(q)) integer(kind=i8) :: i(size(q)) call cumsum(q,qc) !qc=cumsum(q) m = size(q) rn=randn(1) u = ( [(l, l = 0, m-1)] + rn(1))/m i=0; k=1 do j=1,m do while (qc(k) **Author:** ![stavros](https://avatars.discourse-cdn.com/v4/letter/s/51bf81/32.png) [@stavros](https://fortran-lang.discourse.group/u/stavros)\ **Post date:** [October 2, 2021, 7:22am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/13 "2021-10-02T07:22:03Z") Oh, no offense at all!!! I am actually agreeing with you. My comment was only about _elegant_ 😀 In my actual code, which by the way is in a similar field to yours (Monte Carlo and particles) I use nothing functional within the main loops, perhaps only during the preprocessor part where speed doesn’t really matter, and even then not for large arrays. Of course, there is nothing inherently wrong with the functional style, the problem is on the compiler’s side. Such expressions should be automatically vectorized/unrolled. What really happens is that every time such expression is called a temporary array is created. The difference in speed is in orders of magnitudes. By the way, to further optimize use only subroutines, for \*some reason code generated from subroutines is better optimized resulting in minor speed improvement, it is noticeable in large sizes. \*Extra code is generated to reallocate the allocatable array. ---
**Author:** ![themos](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/themos/32/521_2.png) [@themos](https://fortran-lang.discourse.group/u/themos)\ **Post date:** [October 2, 2021, 9:40am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/14 "2021-10-02T09:40:52Z")
> **[Prefix sum](https://en.wikipedia.org/wiki/Prefix_sum)** > > In computer science, the prefix sum, cumulative sum, inclusive scan, or simply scan of a sequence of numbers x0, x1, x2, ... is a second sequence of numbers y0, y1, y2, ..., the sums of prefixes (running totals) of the input sequence: > For instance, the prefix sums of the natural numbers are the triangular numbers: > Prefix sums are trivial to compute in sequential models of computation, by using the formula yi = yi − 1 + xi to compute each output value in sequence order. However, despite their ea... ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [October 2, 2021, 11:17am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/15 "2021-10-02T11:17:18Z")
CRquantum, I am afraid that this thread illustrates why we should heed [Knuth’s warning](http://wiki.c2.com/?PrematureOptimization), “Premature optimization is the root of all evil”. You focused on a small part (i.e., the formation of **cumsum** ) of the overall algorithm, and until you posted #7, we did not know what you planned to do with the computed cumsum array. But now, it turns out that even if you could compute that array in _zero time,_ your overall calculation is going to suffer from the slowness of **linear search**. The cumsum array is already in increasing order (strictly so, if none of the elements of q are zero), so you could do a **binary search** for u(j) in the array qc. Further study of the values that you may expect for the arrays u and qc could enable the binary search to be conducted with a segment of the array qc rather than with the entire array, … ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 2, 2021, 4:55pm UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/16 "2021-10-02T16:55:23Z")
@mecej4 Thanks. I got you. The piece I show is a small piece in my code. True, I totally agree, one should not spend 90% of the time optimizing a part of the code which at most give like 1% speedup. Uhm, but if there is a chance to speedup each piece and learn something from you guys, why not? 🙂 The thing I mainly focused, actually is something else. But that may be a little off topic now (I may put these issues into smaller specific question and ask you guys one by one, as @VladimirF always reminds me). If I may, The true issue, is that perhaps you may ask, why do you need 10^6 particles to represent a distribution? In the code, bringing down the number of particle from 10^6 to 10^3 or so, is where true speedup occur. Because in my code that reduces the number of stochastic differential equations (SDE) to be solved after all. The subroutine sysresample is subroutine for resampling. ``` sysresample(q,i) ``` q is an array size 10^6, each of its value means the normalized weight of a particle. I need to resample the new 10^6 particles according to weights. The i is the index to help me choose particles to form the new resampled array x from the old x. like below, ``` q = q/sum(q) call sysresample(q,idx) x = x(idx,:) ! the x on the LHS is the resampled x, using the index from the x on the RHS. ``` I believe the sysresample subroutine I took from Thomas Schon, [Software and data sets | Thomas Schön](http://user.it.uu.se/~thosc112/research/software.html) is at most N\*log(N) algorithm (which is not too bad) where N mean the number of particles which is the size of the array. The Carpenter etal paper I showed in thread 7 is a perhaps O(N) algorithm. Resampling, or sometimes people call bootstrap sampling or importance-sampling-resampling (SIR). Eg, see a paper, ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/7ce75e1b038735758a1e5397da3d2147113a4ad3.png) The binary search you mentioned is like perhaps bisection method or something like that, is to find some particular values in an array. I could be wrong, but perhaps it is not exactly what I needed here now for resampling. But I do thank you very much indeed for your thoughts which is indeed very helpful, not only for me, but also for everyone here. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [October 3, 2021, 12:15am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/17 "2021-10-03T00:15:10Z")
> [@CRquantum](#): > > Just from business point of view, this elegant code may not be efficient enough. The cost is O( (size(a))^2 ). I think you missed the point here. The [semantics](https://en.wikipedia.org/wiki/Semantics_(computer_science)) of the one-line expression using _implied do_ and the _sum_ intrinsic contain sufficient information, such that an optimizing compiler could recognize this expression as the _cumulative sum_ operation, and choose the fastest set of machine instructions for the target architecture. At least for the hardware (x86-64) and compilers (`gfortran` and `ifort`) available to me, this does not appear to be the case today. A blog post from Daniel Lemire on [“Is C++ worth it”](https://lemire.me/blog/2012/07/23/is-cc-worth-it/) provides an anecdote about cumulative sum performance in C++ and Java with a somewhat-related takeaway message. Now to receive the answer you are looking for, you probably need to ask the more accurate question > How to obtain good performance for the cumulative sum for double precision reals on hardware A with compiler B? ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 3, 2021, 12:58am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/18 "2021-10-03T00:58:33Z")
@ivanpribec One line method is elegant that is obvious. I learned from this line a lot. However, please try the one line method, and the regular do loop cumsum, for size=1000000 array. On the one hand, on my laptop one-line method for size=1000000 array, after I hit enter it does not seem to give me answer. Modern laptop’s CPU is Ghz level, so roughly 10^9 operations per second. For the one line method, 10^6\*10^6=10^12 operations will likely take 10^3 seconds. Considering some vectorization, it may took 10 to 10^2 seconds or so. That is my rough estimation. I am afraid the Intel OneAPI compiler I use may not be that smart to magically reduce a N^2 method to N. On the other hand, the do loop one stupid version gives me answer in less than 0.1 second. I think my question and description is not too unclear for most people, and I am afraid I am not very capable of understanding your accurate question. Perhaps ‘King’ at S.O gives me answer already, > Again, I can be absolutely wrong and I may indeed miss the point, but if you find the one-line method is indeed faster than the stupid do loop subroutine, please let me know with your concrete timing. I feel that may help more. I accept all of your criticism, and I admit most people here have much higher ability in Fortran than me. That is why I come here and learn. But I also think that, if we can just focus on this small topic, how to write cumsum as fast as possible from technical point of view, it may be more positive and useful. But anyway, since there may not be too much room to improve for this cumsum, we may just move on to new topics in other threads. ---
**Author:** ![Beliavsky](https://avatars.discourse-cdn.com/v4/letter/b/ba8739/32.png) [@Beliavsky](https://fortran-lang.discourse.group/u/Beliavsky)\ **Post date:** [October 3, 2021, 1:23am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/19 "2021-10-03T01:23:21Z")
> [@CRquantum](#): > > Again, I can be absolutely wrong and I may indeed miss the point, but if you find the one-line method is indeed faster than the stupid do loop, please let me know with your concrete timing. I think @ivanpribec said that the one-liner is indeed slower (note especially the last sentence below) when he wrote > The [semantics](https://en.wikipedia.org/wiki/Semantics_(computer_science)) of the one-line expression using _implied do_ and the _sum_ intrinsic contain sufficient information, such that an optimizing compiler could recognize this expression as the _cumulative sum_ operation, and choose the fastest set of machine instructions for the target architecture. > > At least for the hardware (x86-64) and compilers ( `gfortran` and `ifort` ) available to me, this does not appear to be the case today ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [October 3, 2021, 2:05am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/20 "2021-10-03T02:05:28Z")
> [@CRquantum](#): > > Uhm, but if there is a chance to speedup each piece and learn something from you guys, why not? That statement opens up a teachable moment. The answer is, “[Because Worse is Better](https://cs.stanford.edu/people/eroberts/cs201/projects/2010-11/WorseIsBetter/index.php/Worse-is-better.html)”. In other words, if I have to spend considerable amounts of time describing and learning good techniques for each sub-task, by the time I learn all the techniques and put together a program to solve the complete problem, I may be too old to care or someone else may have solved it much sooner. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 3, 2021, 2:12am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/21 "2021-10-03T02:12:23Z")
Thank you @Beliavsky . If @ivanpribec could simply add However, and delete At least, in the beginning of the 2nd paragraph, such that " … However, for the hardware (x86-64), … … " Then the logic seems more clear and I can understand where my point is missing more easily. So he is saying that he believe the one-liner gives the compiler enough information such that the compiler ‘intelligently’ know this is ‘cumulative sum’. However, in reality compilers are not that smart. Actually I am not missing the point entirely. But feel free to say anything @ivanpribec , I am not really bothered by any criticism, as long as the logic is perhaps a little bit more clear. (9th floor the comment by @mecej4 is more clear to me.) I know one-liner concise code is elegant and usually gives better performance, particularly when deal with array operations which can be vectorized. See here, > This part by [veryreverie](https://stackoverflow.com/users/8559300/veryreverie) ```auto do i=1,4 ts(i) = t + sum(as(:i-1,i)) * h xs(i) = x + dot_product(as(:i-1,i), ks(:i-1)) ks(i) = h * (fi(xs(i)) + gi(xs(i))*warray(i)) enddo ``` I understand dot\_product(as(:i-1,i), ks(:i-1)) is efficient. In that code, nothing is repeatedly calcuated more than twice. But in the cumsum case, because you know there is some dependence there, like in order to go to step k+1, need to calculate step k first. Therefore although I wish, I do not know how the compiler is intelligent enough to realize it is a cumsum. Unless the engineer perhaps particularly hard-code something in the compiler. But that can be another topic for compiler engineer how to more intelligently optimized user’s code, particularly realize the logic behind the code. However, just for cumsum problem, even if the compiler realize it is cumsum, it will still calculate the cumsum use the stupid loop method perhaps. So it may come back the question again, how to write the fastest cumsum. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [October 3, 2021, 3:19am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/23 "2021-10-03T03:19:19Z")
I have spent two evenings with the cumulative sum in hopes of providing some closure to this thread. To this effect I’ve ignored the warning from Knuth and the many earlier threads on Discourse which have discussed how micro-benchmarking rarely provides good evidence. The way the question from the OP is framed, > is there faster way to do cumsum what perhaps fully utilize vectorization or something? tells us that he’s interested in architecture-specific solutions. The answer will obviously depend upon the hardware, compiler, compiler settings, array size, memory layout, and other factors making it difficult to give any kind of general answer. Best I can do is try with experimentation. In post #10, @mecej4 has brought up Knuth’s warning, which also appears as the epigraph at the beginning of Chapter 1 in the book [“Scientific Software Design”](https://www.cambridge.org/de/academic/subjects/engineering/engineering-general-interest/scientific-software-design-object-oriented-way?format=PB&utm_source=shortlink&utm_medium=shortlink&utm_campaign=rouson) by Rouson, Xia & Xu. Reading the chapter further, Rouson & coauthors introduce the central thesis of the book: > **Your Time is Worth More Than Your Computer’s Time** > > Total solution time and cost can be reduced in greater proportion by reducing \> development time and costs than by reducing computing time and costs. If you are developing an application where performance is critical, I can only suggest you first profile your code using a tool like [gprof](https://ftp.gnu.org/old-gnu/Manuals/gprof-2.9.1/html_mono/gprof.html) or the [Intel® VTune™ Profiler](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/vtune-profiler.html#gs.c9fc7e). This will help you to identify the sections of your code which occupy most of the runtime, and where you should focus your efforts to gain speed-up. * * * Before you reject the canonical implemention with the do loop, it’s fair to do a comparison with other programming languages. I’ve chosen the following three: - C++: [`std::partial_sum`](https://en.cppreference.com/w/cpp/algorithm/partial_sum) - Python: [`numpy.cumsum`](https://numpy.org/doc/stable/reference/generated/numpy.cumsum.html) - MATLAB: [`cumsum`](https://uk.mathworks.com/help/matlab/ref/cumsum.html) I’ve placed all the benchmark codes I used here: > [GitHub - ivan-pi/cumsum\_benchmark: Some benchmarks for the cumulative sum of an array](https://github.com/ivan-pi/cumsum_benchmark) Results I got for C++, Python, and MATLAB are the following: ![ref](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/2b174d90bac0b153f59c3935b8a556799885fc54.png) As you can see the C++ `std::partial_sum` implementation leads with just over 1000 million doubles per second. MATLAB is not trailing very far behind. For Fortran `cumsum`, I’ve prepared some variations following the [C++ codes of Prof. Daniel Lemire](https://github.com/lemire/Code-used-on-Daniel-Lemire-s-blog/tree/master/2012/07/23/). Hopefully, there are no mistakes left. With gfortran I get the following results (columns correspond to array sizes N = 1000, 10000, 100000, 1000000, rows are distinct trials): ```auto # Compiler version: GCC version 10.3.0 # Compiler options: -mtune=generic -march=x86-64 -O3 -Wall # do_loop 391.389 828.912 1024.947 749.708 994.036 1066.553 1028.987 748.018 979.432 1006.239 967.090 753.746 1037.344 1070.320 1039.696 757.796 1018.330 1073.768 1041.558 753.090 1057.082 1074.691 1025.178 759.203 1009.082 1071.582 1033.506 747.828 # step 2 1592.357 1639.882 1334.971 746.350 1515.152 1682.369 1351.461 748.210 1605.136 1554.243 1352.174 665.873 1420.455 1618.123 1245.051 750.325 1517.451 1632.120 1347.073 744.820 1305.483 1500.600 1364.647 751.274 1515.152 1670.844 1322.576 735.441 # step 2 unrolled 1027.749 1040.908 1029.569 639.030 1001.001 1024.695 1014.641 701.242 1025.641 1073.537 1016.085 627.270 1018.330 996.711 933.158 623.454 1003.009 1032.418 864.655 686.328 1061.571 1066.780 1027.359 749.817 1060.445 1041.775 1036.506 753.150 ``` before experiencing a crash with the 3-step method. As you can see the simple do loop provides equal performance as the C++ standard library. The `step2` version is able to achieve a 50 % performance increase for small array sizes, but with large array sizes, the advantage seems to vanish. The second two step version `step2_unrolled` provides no benefit compared to the canonical do loop. Moral of the story is - don’t try to be smarter than the compiler, unless you are a low-level performance tuning engineer with intimate knowledge of the cache organization and CPU instructions sets. I also made a run with the Intel Fortran compiler: ```auto # Compiler version: Intel(R) Fortran Intel(R) 64 Compiler Classic for applicati ons running on Intel(R) 64, Version 2021.3.0 Build 20210609_000000 # Compiler options: -warn all -O2 -o f_cumsum # do_loop 280.269 443.027 524.172 458.147 532.481 532.992 514.099 490.898 513.611 517.813 510.485 496.566 534.474 522.166 524.929 475.125 518.941 517.518 506.827 497.843 533.333 519.292 516.054 484.909 533.333 522.384 516.316 482.826 # step 2 1216.545 1280.082 1071.432 719.369 1485.884 1548.947 1321.528 707.293 1479.290 1604.879 1374.627 652.008 1519.757 1534.684 1351.936 683.041 1519.757 1546.551 1353.858 742.416 1524.390 1540.357 1306.626 708.752 1483.680 1586.798 1366.139 651.351 # step 2 unrolled 326.584 333.400 344.710 334.904 350.877 336.655 345.354 334.240 340.252 351.136 352.399 338.091 356.379 344.187 345.642 331.063 344.353 344.982 344.578 328.011 339.789 335.435 328.517 325.161 344.234 333.901 332.197 323.531 # step 3 993.049 987.069 900.163 623.206 1199.041 1139.601 935.121 632.438 1447.178 1176.471 972.328 697.442 1552.795 1265.342 1035.540 691.362 1468.429 1104.850 1039.685 653.648 1108.647 1159.420 940.230 632.437 1335.113 1144.427 1054.863 708.983 # step four 585.823 659.979 614.881 286.593 683.527 708.617 635.300 294.837 796.813 726.111 624.836 288.141 772.798 738.007 595.745 275.712 816.327 745.268 649.199 255.521 653.595 667.824 539.110 264.546 788.644 745.935 642.422 282.847 ``` Here the `step2_unrolled` version gives performance roughly equal to numpy. The basic `do_loop` version is 50 % slower than gfortran. The `step2` and `step3` versions appear to perform reasonably well. The `step_four` version is taken from the Stack Exchange thread linked in the original post, but it also doesn’t provide much benefit against the simple do loop. In conclusion based upon my limited and imperfect evidence, I agree with @mecej4’s wisdom, “Because Worse is Better”. Using the simple `do` loop is a good return on the invested coding effort. Given the 6 hours I just invested in preparing the code for this post, I don’t think the 50 % increase in two-step unrolled version is really worth it. I would probably prefer to wait a few years, and buy a new CPU. @CRquantum, please feel free to expand the code I put on GitHub with your own `cumsum` variations or borrow it for other performance-related questions. * * * Concerning the implied do loop suggested by Brad (@everythingfunctional), the performance is in fact terribly slow, that I could not afford to wait for the results of the N \> 10^4 cases. As indicated by @CRquantum, the work appears to scale as O(N^2) introducing a severe performance penalty. For N = 1000 the implied do version is roughly 500 times slower for the smallest array size, and drops another order of magnitude for N = 10000: ```auto # implied_do 2.250 0.211 2.245 0.208 2.166 0.206 2.148 0.208 2.086 0.210 2.204 0.212 2.238 0.210 ``` Edit: the results were measured on a ThinkPad T530 laptop with an Intel(R) Core™ i5-3320M (first released 2012). It would be interesting to learn what speed new processors achieve, and whether the gains already outperform the tweaks like loop unrolling. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 3, 2021, 5:11am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/24 "2021-10-03T05:11:34Z")
Thank you @ivanpribec that is truly remarkable work! 👍 My motivation is just that I saw matlab have cumsum, but cumsum in matlab is just like an API, I do not know its source code. So I wonder in Fortran can we have the fastest cunsum algorithm. I just feel that a fastest algorithm overall, should not depend on hardware too much. I have downloaded your code, and I will read your code carefully. Right now, you may use ``` procedure() :: ``` so that there is no need to write the interface blocks. Eg, in your fortran file line 198, ``` subroutine run_test(ntrials,nreps,n,cumsum) integer, intent(in) :: ntrials, nreps integer, intent(in) :: n(:) interface subroutine cumsum(a,b) import dp real(dp), intent(in) :: a(:) real(dp), intent(out) :: b(:) end subroutine end interface ``` I believe you may do ``` subroutine run_test(ntrials,nreps,n,cumsum) use cumsum_mod integer, intent(in) :: ntrials, nreps integer, intent(in) :: n(:) procedure(cumsum_do_loop) :: cumsum ``` The one line ``` procedure(cumsum_do_loop) :: cumsum ``` will automatically take care of the interface for the cumsum. Actually I think you may pick any name of your cumsum subroutines inside the () in procedure(), because they have the same interface. I just ran your code, here is the results. I am using Thinkpad P72 with Xeon-2186M CPU, perhaps it may of some help for your to see how the performance of CPU increased over the years (I do not really think single core performance increased much since Intel’s Haswell architecture) ![image](https://global.discourse-cdn.com/free1/uploads/fortran_lang/original/1X/aeff3b7ee791d6760cd1b30d0b218fa58b094da9.png) Your subroutine cumsum\_step2(a, b) looks like the best cumsum. The block trick is new to me, and it looks very interesting. Perhaps it fully takes advantage of AVX2, perhaps for CPU with AVX512, step2 may be even better. ``` subroutine cumsum_step2(a, b) real(dp), intent(in) :: a(:) real(dp), intent(out) :: b(:) real(dp) :: s0, s1 integer :: i, n s0 = 0.0_dp s1 = 0.0_dp n = size(a) do i = 1, n, 2 block real(dp) :: x0, x1 x0 = a(i) x1 = a(i+1) s0 = s1 + x0 s1 = s1 + (x1 + x0) b(i) = s0 b(i+1) = s1 end block end do if (mod(n,2) /= 0) then b(n) = b(n-1) + a(n) end if end subroutine ``` ---
**Author:** ![mecej4](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/mecej4/32/855_2.png) [@mecej4](https://fortran-lang.discourse.group/u/mecej4)\ **Post date:** [October 3, 2021, 9:33am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/25 "2021-10-03T09:33:50Z")
Here is a modified version of the main loop in CRQuantum’s sysresample() subroutine in his post #7, ```auto k=1 qcum = q(k) do j=1,m do while (qcum < u(j)) k = k + 1 qcum = qcum+q(k) enddo i(j) = k enddo ``` _I have not tested the code, so it may need correction_. Note that it makes the **qc** array and, therefore, the function (or subroutine) **cumsum** superfluous. If this version works well, we have in this thread discussed some great answers to the wrong question as far as the task of doing CRQuantum’s Monte Carlo simulation is concerned. On the other hand, if there is sufficient demand for **cumsum** , it could be added to the list of Fortran intrinsic functions. ---
**Author:** ![CRquantum](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/crquantum/32/730_2.png) [@CRquantum](https://fortran-lang.discourse.group/u/CRquantum)\ **Post date:** [October 3, 2021, 2:43pm UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/26 "2021-10-03T14:43:36Z")
@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](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) **[cumsum – Fortran Program](https://leonfoks.github.io/coretran/interface/cumsum.html)** > > Modern well documented core Fortran routines for cleaner, smarter code. > > > //github.com/leonfoks/coretran/tree/master/src/maths * * * 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](https://gitlab.com/CRquantum/sde-pf/-/tree/master) 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 > [@A small question in John Burkardt's stochastic\_rk.f90 code](https://fortran-lang.discourse.group/t/a-small-question-in-john-burkardts-stochastic-rk-f90-code/1827): > > Dear all, I have a quick question wrt the famous guy John Burkardt’s stochastic\_rk.f90 code. The main description is here, [https://people.math.sc.edu/Burkardt/f\_src/stochastic\_rk/stochastic\_rk.html](https://people.math.sc.edu/Burkardt/f_src/stochastic_rk/stochastic_rk.html) The stochastic\_rk.f90 code is below, [https://people.math.sc.edu/Burkardt/f\_src/stochastic\_rk/stochastic\_rk.f90](https://people.math.sc.edu/Burkardt/f_src/stochastic_rk/stochastic_rk.f90) Now, In the subroutine for RK4 rk4\_ti\_step ( x, t, h, q, fi, gi, seed, xstar ) in the lines, t4 = t1 + a41 \* h + a42 \* h + a43 \* h x4 = x1 + a41 \* k1 + a42 \* k2 I wonder, should th… 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 > **[GitHub - princemahajan/FLINT: Fortran Library for numerical INTegration of...](https://github.com/princemahajan/FLINT)** > > Fortran Library for numerical INTegration of differential equations - GitHub - princemahajan/FLINT: Fortran Library for numerical INTegration of differential equations 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](https://stackoverflow.com/users/3088138/lutz-lehmann) about SDE is very useful. I also asked Kasdin’s method at > **[Has anyone used Kasdin's Runge-Kutta SDE algorithm?](https://discourse.julialang.org/t/has-anyone-used-kasdins-runge-kutta-sde-algorithm/68198)** > > Dear all, I am doing SDE now. I know Julia has fancy differentialequations.jl. It has several SDE solver like SRIW1 or something like that, which uses adaptive Runge-Kutta (RK). I also found there is a fortran version of fix step RK method for... 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. ---
**Author:** ![ivanpribec](https://yyz2.discourse-cdn.com/free1/user_avatar/fortran-lang.discourse.group/ivanpribec/32/3290_2.png) [@ivanpribec](https://fortran-lang.discourse.group/u/ivanpribec)\ **Post date:** [October 4, 2021, 10:44am UTC](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976/27 "2021-10-04T10:44:14Z")
> [@CRquantum](#): > > Uhm, however the calculation of cumsum is still there in the loop, and the time cost is the same if not more. I don’t see why it would cost more. By using a scalar variable `qcum` to store the partial sum, the code from mecej4 got rid of a reference to the `qc` array. Instead of needing to fetch both `qc` and `u` values from the cache, only `u` will be fetched and `qcum` will likely remain in one the registers. I would expect this to provide a small but perhaps significant speed-up to the `sysresample` routine. I am happy to be proven wrong. Of course we can only be certain after performing an experiment and/or analyzing the output assembly code. This reminds me of some projects I came across recently: - [Evidence-Based Performance Tuning](https://stair.center/archives/projects/ebt) - [Evidence-Based Software Engineering](http://community.dur.ac.uk/ebse/) [Next page](https://fortran-lang.discourse.group/t/what-is-the-fastest-way-to-do-cumulative-sum-in-fortran-to-mimic-matlab-cumsum/1976.md?page=2)