Results 1 to 4 of 4

Thread: MergeSort based blocksorting algorithm

  1. #1
    Member
    Join Date
    Jun 2009
    Location
    Kraków, Poland
    Posts
    1,471
    Thanks
    26
    Thanked 120 Times in 94 Posts

    MergeSort based blocksorting algorithm

    Here's my newly invented (or not, probably someone did such thing before, as it's really simple) blocksorting algorithm, which I'll probably use as building block in my complex QRotSort algorithm.

    I've tried stuffing suffix.lcp'th byte from each suffix into suffix.lcp, so then sometimes only comparing those bytes could suffice to figure out the suffix order and lcps, but in the end the gain was really small (and complexity is higher) and probably some other technique will bring some improvement.

    Currently I'm thinking about multi-way MergeSort, but that will be tricky. Also the number of ways is yet to be decided.

    Code:
    /* 
     * File:   main.cpp
     * Author: Piotr Tarsa
     *
     */
    
    #include <cassert>
    #include <cstdio>
    #include <cstdlib>
    #include <ctime>
    
    #include <algorithm>
    
    #define COMPUTE_STATISTICS
    
    struct _index_t {
        int32_t index;
        int32_t lcp;
    };
    
    typedef _index_t index_t;
    
    size_t length;
    index_t *SAin;
    index_t *SAout;
    u_int8_t *sourceBuffer;
    u_int8_t *outputBuffer;
    
    int32_t min(int32_t former, int32_t latter) {
        return former < latter ? former : latter;
    }
    
    int32_t clampOnce(int32_t value, int32_t limit) {
        return (value >= limit) ? value - limit : value;
    }
    
    /*
     * 
     */
    int main(int argc, char** argv) {
    
        freopen("report.txt", "w", stderr);
    
        //FILE *sourceFile = fopen("/home/piotrek/Dokumenty/BWT/qsufsort.c", "rb");
        FILE *sourceFile = fopen("/home/piotrek/Pobrane/corpora/enwik/enwik8", "rb");
    
        fseek(sourceFile, 0, SEEK_END);
        length = ftell(sourceFile);
        const int32_t DepthLimit = min(8000000, length);
        printf("File length: %ld\n", length);
        fseek(sourceFile, 0, SEEK_SET);
    
        SAin = new index_t[length];
        SAout = new index_t[length];
    
        sourceBuffer = new u_int8_t[length];
        outputBuffer = new u_int8_t[length];
    
        size_t readBytes = fread(sourceBuffer, sizeof (u_int8_t), length, sourceFile);
        assert(readBytes == length);
    
        fclose(sourceFile);
    
        for (int32_t i = 0; i < length; i++) {
            SAin[i].index = i;
            SAin[i].lcp = 0;
        }
    
    
    #ifdef COMPUTE_STATISTICS
        int64_t comparisons = 0;
        int64_t quickChecks = 0;
    #endif
    
        int32_t sortedLength = 1;
        while (sortedLength < length) {
    
            int32_t leftStart = 0;
            int32_t leftEnd = sortedLength;
            int32_t rightStart = sortedLength;
            int32_t rightEnd = min(length, sortedLength * 2);
            int32_t outputIndex = 0;
    
            do {
                int32_t toLeftLcp = 0;
                int32_t toRightLcp = 0;
                int32_t lastIndex = 0;
    
                {
                    int32_t firstPtr = SAin[leftStart].index;
                    int32_t secondPtr = SAin[rightStart].index;
                    int32_t equal = 0;
                    for (; equal < DepthLimit; equal++) {
                        if (sourceBuffer[firstPtr] != sourceBuffer[secondPtr]) {
                            break;
                        }
                        firstPtr = ((firstPtr + 1) == length) ? 0 : firstPtr + 1;
                        secondPtr = ((secondPtr + 1) == length) ? 0 : secondPtr + 1;
                    }
    #ifdef COMPUTE_STATISTICS
                    comparisons++;
    #endif
                    if (sourceBuffer[firstPtr] <= sourceBuffer[secondPtr]) {
                        lastIndex = SAin[leftStart].index;
                        toLeftLcp = SAin[leftStart].lcp;
                        leftStart++;
                        toRightLcp = equal;
                    } else {
                        lastIndex = SAin[rightStart].index;
                        toRightLcp = SAin[rightStart].lcp;
                        rightStart++;
                        toLeftLcp = equal;
                    }
                }
    
                while ((leftStart < leftEnd) && (rightStart < rightEnd)) {
    
                    SAout[outputIndex].index = lastIndex;
    
                    if (toLeftLcp > toRightLcp) {
                        SAout[outputIndex++].lcp = toLeftLcp;
                        lastIndex = SAin[leftStart].index;
                        toLeftLcp = SAin[leftStart].lcp;
                        leftStart++;
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    } else if (toLeftLcp < toRightLcp) {
                        SAout[outputIndex++].lcp = toRightLcp;
                        lastIndex = SAin[rightStart].index;
                        toRightLcp = SAin[rightStart].lcp;
                        rightStart++;
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    } else { // toLeftLcp == toRightLcp
                        int32_t firstPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                        int32_t secondPtr = clampOnce(SAin[rightStart].index + toLeftLcp, length);
                        int32_t equal = toLeftLcp;
                        for (; equal < DepthLimit; equal++) {
                            if (sourceBuffer[firstPtr] != sourceBuffer[secondPtr]) {
                                break;
                            }
                            firstPtr = ((firstPtr + 1) == length) ? 0 : firstPtr + 1;
                            secondPtr = ((secondPtr + 1) == length) ? 0 : secondPtr + 1;
                        }
    #ifdef COMPUTE_STATISTICS
                        comparisons++;
    #endif
                        if (sourceBuffer[firstPtr] <= sourceBuffer[secondPtr]) {
                            SAout[outputIndex++].lcp = toLeftLcp;
                            lastIndex = SAin[leftStart].index;
                            toLeftLcp = SAin[leftStart].lcp;
                            leftStart++;
                            toRightLcp = equal;
                        } else {
                            SAout[outputIndex++].lcp = toRightLcp;
                            lastIndex = SAin[rightStart].index;
                            toRightLcp = SAin[rightStart].lcp;
                            rightStart++;
                            toLeftLcp = equal;
                        }
                    }
                }
    
                if (leftStart == leftEnd) {
                    SAout[outputIndex].index = lastIndex;
                    SAout[outputIndex].lcp = toRightLcp;
                    outputIndex++;
    #ifdef COMPUTE_STATISTICS
                    quickChecks++;
    #endif
                    while (rightStart < rightEnd) {
                        SAout[outputIndex++] = SAin[rightStart++];
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    }
                } else { // rightStart == rightEnd
                    SAout[outputIndex].index = lastIndex;
                    SAout[outputIndex].lcp = toLeftLcp;
                    outputIndex++;
    #ifdef COMPUTE_STATISTICS
                    quickChecks++;
    #endif
                    while (leftStart < leftEnd) {
                        SAout[outputIndex++] = SAin[leftStart++];
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    }
                }
    
                leftStart = rightEnd;
                leftEnd = min(leftStart + sortedLength, length);
                rightStart = leftEnd;
                rightEnd = min(rightStart + sortedLength, length);
            } while (rightStart < rightEnd);
    
            while (leftStart < leftEnd) {
                SAout[outputIndex++] = SAin[leftStart++];
    #ifdef COMPUTE_STATISTICS
                quickChecks++;
    #endif
            }
    
            {
                index_t *temp = SAin;
                SAin = SAout;
                SAout = temp;
            }
            sortedLength *= 2;
        }
    
    #ifdef COMPUTE_STATISTICS
        printf("Comparisons:           %11ld\n", comparisons);
        printf("Quick checks:          %11ld\n", quickChecks);
    #endif
    
    
        for (int32_t i = 0; i < length; i++) {
            outputBuffer[i] = sourceBuffer[clampOnce(SAin[i].index + length - 1, length)];
        }
    
        {
            clock_t time = clock();
            for (int32_t i = 1; i < length; i++) {
                int32_t firstPtr = SAin[i - 1].index;
                int32_t secondPtr = SAin[i].index;
                int32_t equal = 0;
                for (; equal < DepthLimit; equal++) {
                    if (sourceBuffer[firstPtr] > sourceBuffer[secondPtr]) {
                        printf("Wrong order on position: %d, equal prefix length: %d\n", i - 1, equal);
                        break;
                    } else if (sourceBuffer[firstPtr] < sourceBuffer[secondPtr]) {
                        break;
                    }
                    firstPtr = ((firstPtr + 1) == length) ? 0 : firstPtr + 1;
                    secondPtr = ((secondPtr + 1) == length) ? 0 : secondPtr + 1;
                }
            }
            printf("Rotation order checking time: %ld\n", clock() - time);
        }
    
        FILE *outputFile = fopen("/home/piotrek/Pulpit/bwtoutput.mergerotsort", "wb");
    
        size_t writtenBytes = fwrite(outputBuffer, sizeof (u_int8_t), length, outputFile);
        assert(writtenBytes == length);
    
        fclose(outputFile);
    
        return 0;
    }

  2. #2
    Member
    Join Date
    Jun 2009
    Location
    Kraków, Poland
    Posts
    1,471
    Thanks
    26
    Thanked 120 Times in 94 Posts

    ThreeWay MergeSort

    OK. I've managed to make three-way version. Unfortunately, it's only little faster, like 10 % on enwik8. But it's always something.


    Code:
    /* 
     * File:   main.cpp
     * Author: Piotr Tarsa
     *
     */
    
    #include <cassert>
    #include <cstdio>
    #include <cstdlib>
    #include <ctime>
    
    #include <algorithm>
    
    #define THREE_WAY_MERGE
    #define INITIAL_SORTED_LENGTH 5
    
    #define COMPUTE_STATISTICS
    
    struct _index_t {
        u_int32_t index;
        u_int32_t lcp;
    };
    
    typedef _index_t index_t;
    
    size_t length;
    index_t *SAin;
    index_t *SAout;
    u_int8_t *sourceBuffer;
    u_int8_t *outputBuffer;
    
    int32_t min(int32_t former, int32_t latter) {
        return former < latter ? former : latter;
    }
    
    int32_t clampOnce(int32_t value, int32_t limit) {
        return (value >= limit) ? value - limit : value;
    }
    
    /*
     * 
     */
    int main(int argc, char** argv) {
    
        freopen("report.txt", "w", stderr);
    
        //FILE *sourceFile = fopen("/home/piotrek/Dokumenty/BWT/qsufsort.c", "rb");
        FILE *sourceFile = fopen("/home/piotrek/Pobrane/corpora/enwik/enwik8", "rb");
    
        fseek(sourceFile, 0, SEEK_END);
        length = ftell(sourceFile);
        const int32_t DepthLimit = length;
        printf("File length: %ld\n", length);
        fseek(sourceFile, 0, SEEK_SET);
    
        SAin = new index_t[length];
        SAout = new index_t[length];
    
        sourceBuffer = new u_int8_t[length];
        outputBuffer = new u_int8_t[length];
    
        size_t readBytes = fread(sourceBuffer, sizeof (u_int8_t), length, sourceFile);
        assert(readBytes == length);
    
        fclose(sourceFile);
    
        for (u_int32_t i = 0; i < length; i++) {
            SAin[i].index = i;
            SAin[i].lcp = 0;
        }
        
    #ifdef COMPUTE_STATISTICS
        int64_t comparisons = 0;
        int64_t quickChecks = 0;
    #endif
    
    #if INITIAL_SORTED_LENGTH > 1
        for (u_int32_t currentStart = 0; currentStart < length; currentStart += INITIAL_SORTED_LENGTH) {
            u_int32_t currentEnd = currentStart + INITIAL_SORTED_LENGTH - 1;
            if (currentEnd >= length) {
                currentEnd = length - 1;
            }
            for (u_int32_t current = currentStart + 1; current <= currentEnd; current++) {
                int32_t currentIndex = SAin[current].index;
                u_int32_t currentLcp = 0;
                int32_t firstPtr = currentIndex;
                u_int8_t firstValue = sourceBuffer[firstPtr];
                u_int32_t j = current;
                do {
    #ifdef COMPUTE_STATISTICS
                    comparisons++;
    #endif
                    u_int32_t newLcp = currentLcp;
                    int32_t secondPtr = clampOnce(SAin[j - 1].index + currentLcp, length);
                    u_int8_t secondValue = sourceBuffer[secondPtr];
                    while ((firstValue == secondValue) && (newLcp < length)) {
                        newLcp++;
                        firstPtr = (firstPtr + 1) == length ? 0 : firstPtr + 1;
                        firstValue = sourceBuffer[firstPtr];
                        secondPtr = (secondPtr + 1) == length ? 0 : secondPtr + 1;
                        secondValue = sourceBuffer[secondPtr];
                    }
                    if (secondValue <= firstValue) {
                        SAin[j - 1].lcp = newLcp;
                        break;
                    }
                    currentLcp = newLcp;
                    do {
                        SAin[j] = SAin[j - 1];
                        j--;
                    } while ((j > currentStart) && (currentLcp < SAin[j - 1].lcp));
                } while ((j > currentStart) && (currentLcp <= SAin[j - 1].lcp));
                SAin[j].index = currentIndex;
                SAin[j].lcp = currentLcp;
            }
        }
    #endif
    
        u_int32_t sortedLength = INITIAL_SORTED_LENGTH;
        while (sortedLength < length) {
    
            clock_t time = clock();
    #ifdef COMPUTE_STATISTICS
            int64_t lastComparisons = comparisons;
            int64_t lastQuickChecks = quickChecks;
    #endif
    
            u_int32_t first = 0;
            int32_t outputIndex = 0;
    
            do {
                int32_t leftStart = first;
                int32_t leftEnd = min(length, leftStart + sortedLength);
                int32_t middleStart = leftEnd;
                int32_t middleEnd = min(length, middleStart + sortedLength);
                int32_t rightStart = middleEnd;
    #ifdef THREE_WAY_MERGE
                int32_t rightEnd = min(length, rightStart + sortedLength);
    #else
                int32_t rightEnd = rightStart;
    #endif
    
                int32_t toLeftLcp = 0;
                int32_t toMiddleLcp = 0;
                int32_t toRightLcp = 0;
                int32_t lastIndex = 0;
    
                if ((leftStart < leftEnd) && (middleStart < middleEnd) && (rightStart < rightEnd)) {
                    int32_t leftPtr = SAin[leftStart].index;
                    int32_t middlePtr = SAin[middleStart].index;
                    int32_t rightPtr = SAin[rightStart].index;
                    u_int8_t leftValue = sourceBuffer[leftPtr];
                    u_int8_t middleValue = sourceBuffer[middlePtr];
                    u_int8_t rightValue = sourceBuffer[rightPtr];
                    int32_t equal = 0;
                    while (true) {
                        if (leftValue != middleValue) {
                            if (leftValue < middleValue) {
                                if (leftValue < rightValue) {
                                    lastIndex = SAin[leftStart].index;
                                    toLeftLcp = SAin[leftStart].lcp;
                                    leftStart++;
                                    toMiddleLcp = equal;
                                    toRightLcp = equal;
                                } else if (leftValue > rightValue) {
                                    lastIndex = SAin[rightStart].index;
                                    toRightLcp = SAin[rightStart].lcp;
                                    rightStart++;
                                    toLeftLcp = equal;
                                    toMiddleLcp = equal;
                                } else {
                                    toMiddleLcp = equal;
                                    if (equal < DepthLimit) {
                                        do {
                                            equal++;
                                            leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                            leftValue = sourceBuffer[leftPtr];
                                            rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                            rightValue = sourceBuffer[rightPtr];
                                        } while ((leftValue == rightValue) && (equal < DepthLimit));
                                    }
                                    if (leftValue <= rightValue) {
                                        lastIndex = SAin[leftStart].index;
                                        toLeftLcp = SAin[leftStart].lcp;
                                        leftStart++;
                                        toRightLcp = equal;
                                    } else {
                                        lastIndex = SAin[rightStart].index;
                                        toRightLcp = SAin[rightStart].lcp;
                                        rightStart++;
                                        toLeftLcp = equal;
                                    }
                                }
                            } else if (middleValue > rightValue) {
                                lastIndex = SAin[rightStart].index;
                                toRightLcp = SAin[rightStart].lcp;
                                rightStart++;
                                toLeftLcp = equal;
                                toMiddleLcp = equal;
                            } else if (middleValue < rightValue) {
                                lastIndex = SAin[middleStart].index;
                                toMiddleLcp = SAin[middleStart].lcp;
                                middleStart++;
                                toLeftLcp = equal;
                                toRightLcp = equal;
                            } else {
                                toLeftLcp = equal;
                                if (equal < DepthLimit) {
                                    do {
                                        equal++;
                                        middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                        middleValue = sourceBuffer[middlePtr];
                                        rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                        rightValue = sourceBuffer[rightPtr];
                                    } while ((middleValue == rightValue) && (equal < DepthLimit));
                                }
                                if (middleValue <= rightValue) {
                                    lastIndex = SAin[middleStart].index;
                                    toMiddleLcp = SAin[middleStart].lcp;
                                    middleStart++;
                                    toRightLcp = equal;
                                } else {
                                    lastIndex = SAin[rightStart].index;
                                    toRightLcp = SAin[rightStart].lcp;
                                    rightStart++;
                                    toMiddleLcp = equal;
                                }
                            }
                            break;
                        } else if (leftValue != rightValue) {
                            if (leftValue > rightValue) {
                                lastIndex = SAin[rightStart].index;
                                toRightLcp = SAin[rightStart].lcp;
                                rightStart++;
                                toLeftLcp = equal;
                                toMiddleLcp = equal;
                            } else {
                                toRightLcp = equal;
                                if (equal < DepthLimit) {
                                    do {
                                        equal++;
                                        leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                        leftValue = sourceBuffer[leftPtr];
                                        middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                        middleValue = sourceBuffer[middlePtr];
                                    } while ((leftValue == middleValue) && (equal < DepthLimit));
                                }
                                if (leftValue <= middleValue) {
                                    lastIndex = SAin[leftStart].index;
                                    toLeftLcp = SAin[leftStart].lcp;
                                    leftStart++;
                                    toMiddleLcp = equal;
                                } else {
                                    lastIndex = SAin[middleStart].index;
                                    toMiddleLcp = SAin[middleStart].lcp;
                                    middleStart++;
                                    toLeftLcp = equal;
                                }
                            }
                            break;
                        }
                        if (equal < DepthLimit) {
                            leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                            leftValue = sourceBuffer[leftPtr];
                            middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                            middleValue = sourceBuffer[middlePtr];
                            rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                            rightValue = sourceBuffer[rightPtr];
                            equal++;
                        }
                        if (equal == DepthLimit) {
                            lastIndex = SAin[leftStart].index;
                            toLeftLcp = SAin[leftStart].lcp;
                            leftStart++;
                            toMiddleLcp = equal;
                            toRightLcp = equal;
                            break;
                        }
                    }
    #ifdef COMPUTE_STATISTICS
                    comparisons++;
    #endif
                } else if ((leftStart < leftEnd) && (middleStart < middleEnd)) {
                    int32_t leftPtr = SAin[leftStart].index;
                    int32_t middlePtr = SAin[middleStart].index;
                    int32_t equal = 0;
                    for (; equal < DepthLimit; equal++) {
                        if (sourceBuffer[leftPtr] != sourceBuffer[middlePtr]) {
                            break;
                        }
                        leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                        middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                    }
    #ifdef COMPUTE_STATISTICS
                    comparisons++;
    #endif
                    if (sourceBuffer[leftPtr] <= sourceBuffer[middlePtr]) {
                        lastIndex = SAin[leftStart].index;
                        toLeftLcp = SAin[leftStart].lcp;
                        leftStart++;
                        toMiddleLcp = equal;
                    } else {
                        lastIndex = SAin[middleStart].index;
                        toMiddleLcp = SAin[middleStart].lcp;
                        middleStart++;
                        toLeftLcp = equal;
                    }
                } else if (leftStart < leftEnd) {
                    lastIndex = SAin[leftStart].index;
                    toLeftLcp = SAin[leftStart].lcp;
                    leftStart++;
    #ifdef COMPUTE_STATISTICS
                    quickChecks++;
    #endif
                }
    
                while ((leftStart < leftEnd) && (middleStart < middleEnd) && (rightStart < rightEnd)) {
    
                    SAout[outputIndex].index = lastIndex;
    
                    if ((toLeftLcp > toMiddleLcp) && (toLeftLcp > toRightLcp)) {
                        SAout[outputIndex++].lcp = toLeftLcp;
                        lastIndex = SAin[leftStart].index;
                        toLeftLcp = SAin[leftStart].lcp;
                        leftStart++;
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    } else if ((toLeftLcp < toMiddleLcp) && (toMiddleLcp > toRightLcp)) {
                        SAout[outputIndex++].lcp = toMiddleLcp;
                        lastIndex = SAin[middleStart].index;
                        toMiddleLcp = SAin[middleStart].lcp;
                        middleStart++;
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    } else if ((toLeftLcp < toRightLcp) && (toMiddleLcp < toRightLcp)) {
                        SAout[outputIndex++].lcp = toRightLcp;
                        lastIndex = SAin[rightStart].index;
                        toRightLcp = SAin[rightStart].lcp;
                        rightStart++;
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    } else {
    #ifdef COMPUTE_STATISTICS
                        comparisons++;
    #endif
                        if (toLeftLcp < toMiddleLcp) { // && toMiddleLcp == toRightLcp
                            int32_t middlePtr = clampOnce(SAin[middleStart].index + toMiddleLcp, length);
                            int32_t rightPtr = clampOnce(SAin[rightStart].index + toMiddleLcp, length);
                            int32_t equal = toMiddleLcp;
                            for (; equal < DepthLimit; equal++) {
                                if (sourceBuffer[middlePtr] != sourceBuffer[rightPtr]) {
                                    break;
                                }
                                middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                            }
                            if (sourceBuffer[middlePtr] <= sourceBuffer[rightPtr]) {
                                SAout[outputIndex++].lcp = toMiddleLcp;
                                lastIndex = SAin[middleStart].index;
                                toMiddleLcp = SAin[middleStart].lcp;
                                middleStart++;
                                toRightLcp = equal;
                            } else {
                                SAout[outputIndex++].lcp = toRightLcp;
                                lastIndex = SAin[rightStart].index;
                                toRightLcp = SAin[rightStart].lcp;
                                rightStart++;
                                toMiddleLcp = equal;
                            }
                        } else if (toLeftLcp > toMiddleLcp) { // toLeftLcp == toRightLcp
                            int32_t leftPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                            int32_t rightPtr = clampOnce(SAin[rightStart].index + toLeftLcp, length);
                            int32_t equal = toLeftLcp;
                            for (; equal < DepthLimit; equal++) {
                                if (sourceBuffer[leftPtr] != sourceBuffer[rightPtr]) {
                                    break;
                                }
                                leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                            }
                            if (sourceBuffer[leftPtr] <= sourceBuffer[rightPtr]) {
                                SAout[outputIndex++].lcp = toLeftLcp;
                                lastIndex = SAin[leftStart].index;
                                toLeftLcp = SAin[leftStart].lcp;
                                leftStart++;
                                toRightLcp = equal;
                            } else {
                                SAout[outputIndex++].lcp = toRightLcp;
                                lastIndex = SAin[rightStart].index;
                                toRightLcp = SAin[rightStart].lcp;
                                rightStart++;
                                toLeftLcp = equal;
                            }
                        } else if (toLeftLcp > toRightLcp) { // && toLeftLcp == toMiddleLcp
                            int32_t leftPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                            int32_t middlePtr = clampOnce(SAin[middleStart].index + toLeftLcp, length);
                            int32_t equal = toLeftLcp;
                            for (; equal < DepthLimit; equal++) {
                                if (sourceBuffer[leftPtr] != sourceBuffer[middlePtr]) {
                                    break;
                                }
                                leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                            }
                            if (sourceBuffer[leftPtr] <= sourceBuffer[middlePtr]) {
                                SAout[outputIndex++].lcp = toLeftLcp;
                                lastIndex = SAin[leftStart].index;
                                toLeftLcp = SAin[leftStart].lcp;
                                leftStart++;
                                toMiddleLcp = equal;
                            } else {
                                SAout[outputIndex++].lcp = toMiddleLcp;
                                lastIndex = SAin[middleStart].index;
                                toMiddleLcp = SAin[middleStart].lcp;
                                middleStart++;
                                toLeftLcp = equal;
                            }
                        } else { // toLeftLcp == toMiddleLcp == toRightLcp
                            int32_t leftPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                            int32_t middlePtr = clampOnce(SAin[middleStart].index + toLeftLcp, length);
                            int32_t rightPtr = clampOnce(SAin[rightStart].index + toLeftLcp, length);
                            u_int8_t leftValue = sourceBuffer[leftPtr];
                            u_int8_t middleValue = sourceBuffer[middlePtr];
                            u_int8_t rightValue = sourceBuffer[rightPtr];
                            int32_t equal = toLeftLcp;
                            while (true) {
                                if (leftValue != middleValue) {
                                    if (leftValue < middleValue) {
                                        if (leftValue < rightValue) {
                                            SAout[outputIndex++].lcp = toLeftLcp;
                                            lastIndex = SAin[leftStart].index;
                                            toLeftLcp = SAin[leftStart].lcp;
                                            leftStart++;
                                            toMiddleLcp = equal;
                                            toRightLcp = equal;
                                        } else if (leftValue > rightValue) {
                                            SAout[outputIndex++].lcp = toRightLcp;
                                            lastIndex = SAin[rightStart].index;
                                            toRightLcp = SAin[rightStart].lcp;
                                            rightStart++;
                                            toLeftLcp = equal;
                                            toMiddleLcp = equal;
                                        } else {
                                            toMiddleLcp = equal;
                                            if (equal < DepthLimit) {
                                                do {
                                                    equal++;
                                                    leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                                    leftValue = sourceBuffer[leftPtr];
                                                    rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                                    rightValue = sourceBuffer[rightPtr];
                                                } while ((leftValue == rightValue) && (equal < DepthLimit));
                                            }
                                            if (leftValue <= rightValue) {
                                                SAout[outputIndex++].lcp = toLeftLcp;
                                                lastIndex = SAin[leftStart].index;
                                                toLeftLcp = SAin[leftStart].lcp;
                                                leftStart++;
                                                toRightLcp = equal;
                                            } else {
                                                SAout[outputIndex++].lcp = toRightLcp;
                                                lastIndex = SAin[rightStart].index;
                                                toRightLcp = SAin[rightStart].lcp;
                                                rightStart++;
                                                toLeftLcp = equal;
                                            }
                                        }
                                    } else if (middleValue > rightValue) {
                                        SAout[outputIndex++].lcp = toRightLcp;
                                        lastIndex = SAin[rightStart].index;
                                        toRightLcp = SAin[rightStart].lcp;
                                        rightStart++;
                                        toLeftLcp = equal;
                                        toMiddleLcp = equal;
                                    } else if (middleValue < rightValue) {
                                        SAout[outputIndex++].lcp = toMiddleLcp;
                                        lastIndex = SAin[middleStart].index;
                                        toMiddleLcp = SAin[middleStart].lcp;
                                        middleStart++;
                                        toLeftLcp = equal;
                                        toRightLcp = equal;
                                    } else {
                                        toLeftLcp = equal;
                                        if (equal < DepthLimit) {
                                            do {
                                                equal++;
                                                middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                                middleValue = sourceBuffer[middlePtr];
                                                rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                                rightValue = sourceBuffer[rightPtr];
                                            } while ((middleValue == rightValue) && (equal < DepthLimit));
                                        }
                                        if (middleValue <= rightValue) {
                                            SAout[outputIndex++].lcp = toMiddleLcp;
                                            lastIndex = SAin[middleStart].index;
                                            toMiddleLcp = SAin[middleStart].lcp;
                                            middleStart++;
                                            toRightLcp = equal;
                                        } else {
                                            SAout[outputIndex++].lcp = toRightLcp;
                                            lastIndex = SAin[rightStart].index;
                                            toRightLcp = SAin[rightStart].lcp;
                                            rightStart++;
                                            toMiddleLcp = equal;
                                        }
                                    }
                                    break;
                                } else if (leftValue != rightValue) {
                                    if (leftValue > rightValue) {
                                        SAout[outputIndex++].lcp = toRightLcp;
                                        lastIndex = SAin[rightStart].index;
                                        toRightLcp = SAin[rightStart].lcp;
                                        rightStart++;
                                        toLeftLcp = equal;
                                        toMiddleLcp = equal;
                                    } else {
                                        toRightLcp = equal;
                                        if (equal < DepthLimit) {
                                            do {
                                                equal++;
                                                leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                                leftValue = sourceBuffer[leftPtr];
                                                middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                                middleValue = sourceBuffer[middlePtr];
                                            } while ((leftValue == middleValue) && (equal < DepthLimit));
                                        }
                                        if (leftValue <= middleValue) {
                                            SAout[outputIndex++].lcp = toLeftLcp;
                                            lastIndex = SAin[leftStart].index;
                                            toLeftLcp = SAin[leftStart].lcp;
                                            leftStart++;
                                            toMiddleLcp = equal;
                                        } else {
                                            SAout[outputIndex++].lcp = toMiddleLcp;
                                            lastIndex = SAin[middleStart].index;
                                            toMiddleLcp = SAin[middleStart].lcp;
                                            middleStart++;
                                            toLeftLcp = equal;
                                        }
                                    }
                                    break;
                                }
                                if (equal < DepthLimit) {
                                    leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                                    leftValue = sourceBuffer[leftPtr];
                                    middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                                    middleValue = sourceBuffer[middlePtr];
                                    rightPtr = ((rightPtr + 1) == length) ? 0 : rightPtr + 1;
                                    rightValue = sourceBuffer[rightPtr];
                                    equal++;
                                }
                                if (equal == DepthLimit) {
                                    SAout[outputIndex++].lcp = toLeftLcp;
                                    lastIndex = SAin[leftStart].index;
                                    toLeftLcp = SAin[leftStart].lcp;
                                    leftStart++;
                                    toMiddleLcp = equal;
                                    toRightLcp = equal;
                                    break;
                                }
                            }
                        }
                    }
                }
    
                if ((leftStart == leftEnd) && (middleStart < middleEnd)) {
                    leftStart = middleStart;
                    leftEnd = middleEnd;
                    toLeftLcp = toMiddleLcp;
                } else if ((middleStart < middleEnd) && (rightStart == rightEnd)) {
                    rightStart = middleStart;
                    rightEnd = middleEnd;
                    toRightLcp = toMiddleLcp;
                }
    
                while ((leftStart < leftEnd) && (rightStart < rightEnd)) {
    
                    SAout[outputIndex].index = lastIndex;
    
                    if (toLeftLcp > toRightLcp) {
                        SAout[outputIndex++].lcp = toLeftLcp;
                        lastIndex = SAin[leftStart].index;
                        toLeftLcp = SAin[leftStart].lcp;
                        leftStart++;
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    } else if (toLeftLcp < toRightLcp) {
                        SAout[outputIndex++].lcp = toRightLcp;
                        lastIndex = SAin[rightStart].index;
                        toRightLcp = SAin[rightStart].lcp;
                        rightStart++;
    #ifdef COMPUTE_STATISTICS
                        quickChecks++;
    #endif
                    } else { // toLeftLcp == toRightLcp
                        int32_t firstPtr = clampOnce(SAin[leftStart].index + toLeftLcp, length);
                        int32_t secondPtr = clampOnce(SAin[rightStart].index + toLeftLcp, length);
                        int32_t equal = toLeftLcp;
                        for (; equal < DepthLimit; equal++) {
                            if (sourceBuffer[firstPtr] != sourceBuffer[secondPtr]) {
                                break;
                            }
                            firstPtr = ((firstPtr + 1) == length) ? 0 : firstPtr + 1;
                            secondPtr = ((secondPtr + 1) == length) ? 0 : secondPtr + 1;
                        }
    #ifdef COMPUTE_STATISTICS
                        comparisons++;
    #endif
                        if (sourceBuffer[firstPtr] <= sourceBuffer[secondPtr]) {
                            SAout[outputIndex++].lcp = toLeftLcp;
                            lastIndex = SAin[leftStart].index;
                            toLeftLcp = SAin[leftStart].lcp;
                            leftStart++;
                            toRightLcp = equal;
                        } else {
                            SAout[outputIndex++].lcp = toRightLcp;
                            lastIndex = SAin[rightStart].index;
                            toRightLcp = SAin[rightStart].lcp;
                            rightStart++;
                            toLeftLcp = equal;
                        }
                    }
                }
    
                if (rightStart < rightEnd) {
                    leftStart = rightStart;
                    leftEnd = rightEnd;
                    toLeftLcp = toRightLcp;
                }
    
                SAout[outputIndex].index = lastIndex;
                SAout[outputIndex].lcp = toLeftLcp;
                outputIndex++;
    
                while (leftStart < leftEnd) {
                    SAout[outputIndex++] = SAin[leftStart++];
    #ifdef COMPUTE_STATISTICS
                    quickChecks++;
    #endif
                }
    
    #ifdef THREE_WAY_MERGE
                first += sortedLength * 3;
    #else
                first += sortedLength * 2;
    #endif
            } while (first < length);
    
            if (outputIndex != length) {
                puts("dsflj");
            }
    
            {
                index_t *temp = SAin;
                SAin = SAout;
                SAout = temp;
            }
    
    #ifdef COMPUTE_STATISTICS
            printf("One pass merge time: %9ld, stride: %11u, comparisons: %11ld, quick checks: %11ld\n",
                    clock() - time, sortedLength, comparisons - lastComparisons,
                    quickChecks - lastQuickChecks);
    #else 
            printf("One pass merge time: %ld, stride: %u\n", clock() - time, sortedLength);
    #endif
    
    #ifdef THREE_WAY_MERGE
            sortedLength *= 3;
    #else
            sortedLength *= 2;
    #endif
    
        }
    
    #ifdef COMPUTE_STATISTICS
        printf("Comparisons:           %11ld\n", comparisons);
        printf("Quick checks:          %11ld\n", quickChecks);
    #endif
    
        //return 0;
    
    
        for (u_int32_t i = 0; i < length; i++) {
            outputBuffer[i] = sourceBuffer[clampOnce(SAin[i].index + length - 1, length)];
        }
    
        if (false) {
            int64_t lcpSum = 0;
            clock_t time = clock();
            for (u_int32_t i = 1; i < length; i++) {
                lcpSum += SAin[i - 1].lcp;
                int32_t leftPtr = SAin[i - 1].index;
                int32_t middlePtr = SAin[i].index;
                int32_t equal = 0;
                for (; equal < DepthLimit; equal++) {
                    if (sourceBuffer[leftPtr] > sourceBuffer[middlePtr]) {
                        printf("Wrong order on position: %d, equal prefix length: %d\n", i - 1, equal);
                        break;
                    } else if (sourceBuffer[leftPtr] < sourceBuffer[middlePtr]) {
                        break;
                    }
                    leftPtr = ((leftPtr + 1) == length) ? 0 : leftPtr + 1;
                    middlePtr = ((middlePtr + 1) == length) ? 0 : middlePtr + 1;
                }
                if (equal != SAin[i - 1].lcp) {
                    puts("błąd");
                }
            }
            printf("Rotation order checking time: %ld\n", clock() - time);
            printf("Average LCP: %lf\n", (double) lcpSum / length);
        }
    
        FILE *outputFile = fopen("/home/piotrek/Pulpit/bwtoutput.mergerotsort", "wb");
    
        size_t writtenBytes = fwrite(outputBuffer, sizeof (u_int8_t), length, outputFile);
        assert(writtenBytes == length);
    
        fclose(outputFile);
    
        return 0;
    }
    Last edited by Piotr Tarsa; 17th August 2011 at 13:36. Reason: Added optional InsertionSort

  3. #3
    Member
    Join Date
    Jul 2011
    Location
    boulder co usa
    Posts
    3
    Thanks
    0
    Thanked 0 Times in 0 Posts

    Thinking multi-core?

    I'm not sure what is the complete context of your code, but would a multi-core parallelism implementation be more helpful than low-level code optimization? Given that almost everyone has two cores now, and servers have 8-64, this would seem to be an obvious thing. Have you looked at using e.g. CILK++ or Intel TBB or OpenMP?

  4. #4
    Member
    Join Date
    Jun 2009
    Location
    Kraków, Poland
    Posts
    1,471
    Thanks
    26
    Thanked 120 Times in 94 Posts
    It's supposed to be a part of a multithreaded concurrent block-sorter. It's pointless to make a threaded algorithm which will be equally performant on four cores than best non-threaded algorithm on one core. I also care about the performance on low number of cores. I am focusing on eg performance/ watt, so both scaling and single-core performance are equally importatnt.

Similar Threads

  1. BWT based ZPAQ model
    By Matt Mahoney in forum Data Compression
    Replies: 0
    Last Post: 17th March 2011, 04:14
  2. New layer 0 - compression algorithm
    By abocut in forum Data Compression
    Replies: 5
    Last Post: 28th May 2010, 01:32
  3. The best algorithm for high compression
    By Wladmir in forum Data Compression
    Replies: 8
    Last Post: 18th April 2010, 14:54
  4. BCM v0.01 - New BWT+CM-based compressor
    By encode in forum Data Compression
    Replies: 81
    Last Post: 9th February 2009, 16:47
  5. PREDICTOR algorithm
    By encode in forum Forum Archive
    Replies: 30
    Last Post: 16th February 2008, 19:28

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •