1

I tried to implement an algorithm which sorts the array r of length DIM*n using an array of length n. I don't see where my code is wrong. I don't get the expected result. The result should look like a space filling Morton curve. But as you can see, the result consists of a lot of zeros. I don't know where they come from? Can you please help me to find the error? Here is my executable code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define DIM 2

// Sort list "r" using list "mInt"
void sort(unsigned int *mInt, double *r, int n){ 
   unsigned int i, j, ph0;
   double ph1, ph2;

   for(i = 1; i <= n-1; i++)
      for(j = 1; j <= n-i; j++)
         if(mInt[j-1] >= mInt[j])
         {
            // 1
            ph1 = r[DIM*(j-1)+0];
            ph2 = r[DIM*(j-1)+1];
            ph0 = mInt[j-1];

            // 2
            mInt[j-1] = mInt[j];
            r[DIM*(j-1)+0] = r[DIM*j+0];
            r[DIM*(j-1)+1] = r[DIM*j+1];

            // 3
            mInt[j] = ph0;
            r[DIM*j+0] = ph1;
            r[DIM*j+1] = ph2;
         }
}

// Create morton key
inline unsigned int mortoncode(unsigned int x, unsigned int y){
    int answer = 0;
    for (unsigned int i = 0; i < (sizeof(unsigned int)* 8)/2; ++i) {
        answer |= (((x & ((unsigned int)1 << i)) << 2*i) | ((y & ((unsigned int)1 << i)) << (2*i + 1)));
    }
    return answer;
}

// Find max / min values
double maxValue(double *r, int n, int d){
    double max = r[d];
    for(int k=0; k<n; k++){
        if(max < r[DIM*k+d]){
            max = r[DIM*k+d];
        }
    }
    return max;
}
double minValue(double *r, int n, int d){
    double min = r[d];
    for(int k=0; k<n; k++){
        if(min > r[DIM*k+d]){
            min = r[DIM*k+d];
        }
    }
    return min;
}

int main(int argc, char **argv){  
    FILE *f = fopen("data.dat", "w"); 
    int n = 100;
    double r[n*DIM];

    // Initialize data
    double x1 = 0;
    double y1 = 0;
    for(int k=0; k<n; k++){
        r[DIM*k+0] = x1;
        r[DIM*k+1] = y1;
        x1 += 0.1;
        if(k % 10 == 0){
            y1 += 0.1;
            x1 = 0;
        }
        printf("%lf %lf\n", r[DIM*k+0], r[DIM*k+1]);
    }

    // Get max/min values
    double rMin[DIM];
    double rMax[DIM];
    for(int d=0; d<DIM; d++){
        rMin[d] = minValue(r, n, d);
        rMax[d] = maxValue(r, n, d);
    }

    // Convert double data to integers
    printf("\n");
    unsigned int rInt[n];
    for(int k=0; k<n; k++){
        for(int d=0; d<DIM; d++){
            int idx=DIM*k+d;
            double map = floor(((2097151)/(rMax[d]-rMin[d]))*r[idx]-rMin[d]);
            rInt[idx] = (int)map;
        }
        printf("%d %d\n", rInt[DIM*k+0], rInt[DIM*k+1]);
    }

    // Convert rInt[x_1,y_1,x_2,y_2,...] to Morton key
    printf("\n");
    unsigned int rMor[n];
    for(int k=0; k<n; k++){
        int idx = DIM*k;
        rMor[k] = mortoncode(rInt[idx+0], rInt[idx+1]);
    }

    // Sort data 
    sort(rMor, r, n);

    for(int k=0; k<n; k++){
        printf("%lf %lf\n", r[DIM*k+0], r[DIM*k+1]);
        fprintf(f, "%lf, %lf\n", r[DIM*k+0], r[DIM*k+1]);

    }

    return 0;
}
11
  • unsigned int rInt[n] is not large enough. You probably get error before it reaches sort Commented Nov 26, 2016 at 1:15
  • @cdlane data.dat is an output file, where I save my data. You see the output also in the terminal. See last part of code printf("%lf %lf\n", r[DIM*k+0], r[DIM*k+1]); and fprintf(f, "%lf, %lf\n", r[DIM*k+0], r[DIM*k+1]);. Commented Nov 26, 2016 at 1:21
  • @cdlane That would not help very much. It is just the defective version of array r. Something with the sort function must be wrong. But I don't find the error. Commented Nov 26, 2016 at 1:36
  • BTW, what is rand01pm()? A special function or a global replace gone bad? I checked @BarmakShemirani's rInt[n] but I don't see map values much greater than 2.0e+06 so they should fit into an unsigned int. Commented Nov 26, 2016 at 1:39
  • It should be rInt[n*DIM] and rMor[n*DIM] or more, otherwise you get overflow. Also this first part of your code is not relevant to the question. You could just fill the arrays with random values and focus on the sort. Commented Nov 26, 2016 at 1:45

1 Answer 1

3

I believe that @BarmakShemirani had the answer half a dozen comments ago, you claim:

rInt is n*DIM big.

But you wrote:

unsigned int rInt[n];

Fixing that, I tested your sort() routine by not using it but rather by making r and rMor into a single array of structs and calling qsort() on it. They basically produce the same result except for duplicate indexes where one puts them out reversed relative to the other:

           qsort                          your sort
index      r0       r1         index      r0       r1
2456659099 0.400000 0.500000   2456659099 0.400000 1.000000
2456659099 0.400000 1.000000   2456659099 0.400000 0.500000

The modified code:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define DIM 2

typedef struct element {
    unsigned int m;
    double r[DIM];
} ELEMENT;

// Sort element structure on member 'm'

int comparator(const void *p, const void *q) {

    ELEMENT *a = (ELEMENT *) p;
    ELEMENT *b = (ELEMENT *) q;

    return (a->m > b->m) - (a->m < b->m); // compare idiom
}

// Create morton key
unsigned int mortoncode(unsigned int x, unsigned int y) {
    int answer = 0;

    for (unsigned int i = 0; i < (sizeof(unsigned int) * 8) / 2; i++) {
        answer |= (((x & (1u << i)) << 2 * i) | ((y & (1u << i)) << (2 * i + 1)));
    }

    return answer;
}

// Find max / min values
double maxValue(ELEMENT data[], int n, int d) {
    double max = data[0].r[d];

    for (int k = 0; k < n; k++) {
        if (max < data[k].r[d]) {
            max = data[k].r[d];
        }
    }

    return max;
}

double minValue(ELEMENT data[], int n, int d) {
    double min = data[0].r[d];

    for (int k = 0; k < n; k++) {
        if (min > data[k].r[d]) {
            min = data[k].r[d];
        }
    }

    return min;
}

int main(int argc, char **argv) {  
    FILE *f = fopen("data.dat", "w"); 
    int n = 100;
    ELEMENT data[n];

    // Initialize data
    double x1 = 0;
    double y1 = 0;

    for (int k = 0; k < n; k++) {
        data[k].r[0] = x1;
        data[k].r[1] = y1;
        x1 += 0.1;
        if (k % 10 == 0) {
            y1 += 0.1;
            x1 = 0;
        }
        printf("%lf %lf\n", data[k].r[0], data[k].r[1]);
    }
    printf("\n");

    // Get max/min values
    double rMin[DIM];
    double rMax[DIM];

    for (int d = 0; d < DIM; d++) {
        rMin[d] = minValue(data, n, d);
        rMax[d] = maxValue(data, n, d);
    }

    // Convert double data to integers

    unsigned int rInt[DIM * n];

    for (int k = 0; k < n; k++) {
        for (int d = 0; d < DIM; d++) {
            int idx = DIM * k + d;
            double map = floor((2097151 / (rMax[d] - rMin[d])) * data[k].r[d] - rMin[d]);
            rInt[idx] = (int) map;
        }

        printf("%d %d\n", rInt[DIM * k + 0], rInt[DIM * k + 1]);
    }
    printf("\n");

    // Convert rInt[x_1, y_1, x_2, y_2, ...] to Morton key

    for (int k = 0; k < n; k++) {
        int idx = DIM * k;
        data[k].m = mortoncode(rInt[idx + 0], rInt[idx + 1]);
    }

    // Sort data 
    qsort(data, n, sizeof(ELEMENT), comparator);

    for (int k = 0; k < n; k++) {
        printf("%u %lf %lf\n", data[k].m, data[k].r[0], data[k].r[1]);
        fprintf(f, "%lf, %lf\n", data[k].r[0], data[k].r[1]);
    }

    return 0;
}

Any time you find yourself creating parallel arrays like r and rMor, it's usually a sign that your missing a real data structure.

Sign up to request clarification or add additional context in comments.

4 Comments

In fairness to the OP, he says that he has not yet learned about structures in C. That doesn't stop it being the better way to process the data. Your use of double r[DIM]; in the structure instead of the double x; double y; I suggested in a comment is an improvement; it keeps the minValue and maxValue functions simpler.
Thank you @JonathanLeffler for your comment. I worked with the assumption that DIM could increase (even though not all the code is currently set up for that). I switched to a structure approach to use qsort() to test the OP's own sort() upon which he was pinning the blame. The differences I noted when two index elements compared the same was a missed opportunity to discuss the concept of a stable sort.
@cdlane Thank you for your comment. Unfortunately the result of your code does not fit to my expectations. The resulting elements of r after the sorting do not create a Z-Curve when I plot them with lines. I tried the following for nine points: For the highest coordinates in r, we get the largest Morton-Keys and vice versa. When I compile your code the highest Morton-Key is not longer associated with the highest coordinates. Do you see why this is?
@cdlane Ok, I checked your code again and I think it does the right thing. But why do I not get an Z-Curve if I plot the sorted r. All I get is a really mixed up plot. Can someone help me to understand what I am doing wrong?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.