Loop-Level Prefix-Sum: Messed up partial sums

Question

#include <omp.h>
#include <stdio.h>


void prefix(int n, int A[n]){
    
    //end of recursion
    if(n==1) return;
    
    int n2 = n/2;
    
    //Array split
    int L[n2];
    int R[n-n2];
    
    #pragma omp parallel for
    for(int i=0; i<n2; i++){
        L[i]=A[i];
    }
    

    #pragma omp parallel for
    for(int i=n2; i<n-n2; i++){
        R[i-n2]=A[i];
    }
         
    
    //recursion
    #pragma omp task
    prefix(n2,L);
    
    #pragma omp task
    prefix(n-n2,R);
    
    #pragma omp taskwait
    #pragma omp parallel for
    for(int i = 0; i<n2; i++){
        A[n2+i]=A[n2]+A[n2+i];
    }
    for(int i = 0; i<n; i++){
        printf("%d\n",A[i]);
    }
    printf("END\n");
}

int main(){
    
    int A[4] = {1,2,3,4};
    prefix(4,A);
    
    for(int i=0;i<4;i++){
        printf("%d",A[i]);
    }
}

Outputs: 1 END 32 END 1 4 END 16 END 32 END 16 0 END 1 2 6 7 END 1267

So I know that the partial recursion sums are messed up due to 32 in the first right half. But I am currently struggling to find the reason. The initial algorithm uses arrays from 1 to n. Is it something with the ranges?

Thanks for any help.

Edit:

Now I get:

1 2 END 6684773 7864425 END 1 2 3 7 END 1237

What is it causing to fail to calculate the prefix sum?

Solved it now:

#include <omp.h>
#include <stdio.h>


void prefix(int n, int *A, int start, int end){
    
    //end of recursion
    if(n==1) return;
    
    int n2 = n/2;
    
    //recursion
    #pragma omp task
    prefix(n2,A, start, end-n2);
    
    #pragma omp task
    prefix(n-n2,A,end-n2+1, end);   
    #pragma omp taskwait
   
    #pragma omp parallel for
    for(int i = 0; i<n2; i++){
        A[end-n2+1+i]=A[end-n2]+A[end-n2+1+i];
    }
    for(int i = 0; i<n; i++){
        printf("%d\n",A[i]);
    }
    printf("END\n");
}

int main(){
    
    int A[4] = {1,1,1,1};
    prefix(4,A,0,3);
    
    for(int i=0;i<4;i++){
        printf("%d",A[i]);
    }
}

If `n` is zero, you should return zero, not `A[0]` (which is an out-of-bounds access in that case, assuming zero-sized arrays are meaningful). If `n` is one, then you can return `A[0]`. — Tom Karzes, Aug 21 '21 at 02:26
The `prefix` function shouldn't be returning anything, since all three calls to `prefix` ignore the return value. Declare the function as `void prefix(int n, int A[n]);` and change `return A[0];` to `return;`. — user3386109, Aug 21 '21 at 02:52
The code doesn't properly handle odd values of `n`. Keep in mind that integer division truncates towards 0. For example, if `n=3` then `L[n/2]` and `R[n/2]` are in fact `L[1]` and `R[1]`. That's a problem since the starting value of `n` is 4, so the recursive calls are `n=2` and then `n=1`, and 1 is odd. — user3386109, Aug 21 '21 at 02:58
Compute `n2 = n/2` once, at the top of the function. The two arrays will then have lengths `n2` and `n - n2`. Note that `n - n2` will equal `n2` if `n` is even, and `n2+1` if `n` is odd. — Tom Karzes, Aug 21 '21 at 03:15
What exactly `prefix` function should do? Please clarify it. Besides indexing problems, recursive calls using local arrays (`L`,`R`) does not have any effect on array `A` in `main`. — Laci, Aug 21 '21 at 05:10
So my question: Why does my code result in outputting 16 and 32? — Rapiz, Aug 21 '21 at 20:37

score 0 · Answer 1 · answered Aug 22 '21 at 01:05

#include <omp.h>
#include <stdio.h>


void prefix(int n, int *A, int start, int end){
    
    //end of recursion
    if(n==1) return;
    
    int n2 = n/2;
    
    //recursion
    #pragma omp task
    prefix(n2,A, start, end-n2);
    
    #pragma omp task
    prefix(n-n2,A,end-n2+1, end);   
    #pragma omp taskwait
   
    #pragma omp parallel for
    for(int i = 0; i<n2; i++){
        A[end-n2+1+i]=A[end-n2]+A[end-n2+1+i];
    }
    for(int i = 0; i<n; i++){
        printf("%d\n",A[i]);
    }
    printf("END\n");
}

int main(){
    
    int A[4] = {1,1,1,1};
    prefix(4,A,0,3);
    
    for(int i=0;i<4;i++){
        printf("%d",A[i]);
    }
}

Laci · Answer 2 · 2021-08-22T08:57:37.500

0

Note that the program in your answer does not give correct result if n is odd. For this simple task you do not need a recursive algorithm unless it is your homework and you have to use a recursive algorithm and OpenMP tasks. You should use a much simpler algorithm such as

void prefix(int n, int* A){    
    if (n>1) {
        int sum=A[0];
        for(int i=1; i<n;i++){
            sum+=A[i];
            A[i]=sum;
        }
    }
}

edited Aug 22 '21 at 08:57

answered Aug 22 '21 at 04:57

Laci

2,738
1
13
22

Wouldn't it work for odd numbers by just enlarging the initial array with a zero and remove it at the end? – Rapiz Aug 22 '21 at 15:22
No. For `int A[] = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1};` your code gives: `1,1,1,3,4,4,6,7,7,9,10,10,12,13,13,15,16,16,18,19,19,21,25,26,27`. – Laci Aug 22 '21 at 17:55
1

If you need a parallel algorithm to calculate prefix-sum you can find many on the net, e.g.: http://www2.hawaii.edu/~nodari/teaching/f17/slides/L06-segmented-prefix-sums.pdf – Laci Aug 22 '21 at 18:35
Would work to enlarge it with 0 until the length is a potency of 2. – Rapiz Aug 23 '21 at 19:00

Loop-Level Prefix-Sum: Messed up partial sums

2 Answers2