symmetric.c

// gcc -march=native -g -std=gnu99 -Wall -Wextra -O3 symmetric.c -o symmetric -DUSE_ALG
//   (where USE_ALG is one of USE_NATE, USE_KARIM, USE_BASIC, or USE_CONDITIONAL

// Or if using https://code.google.com/p/likwid/ with -m markers:
// gcc -march=native -g -std=gnu99 -Wall -Wextra -O3 symmetric.c -o symmetric -DLIKWID -llikwid -lpthread -lm -DUSE_ALG

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

#define DEFAULT_NUM_ROWS 10000
#define DEFAULT_NUM_REPEAT 10
#define DEFAULT_DENSITY 0.5

void dieUsage(char *programName) {
    printf("Usage: %s [numRows] [numRepeat] [density]\n", programName);
    printf("  Convert an upper right covariance matrix to be symmetric.\n");
    printf("\n");
    printf("  numRows:   size of square matrix (default %d)\n", DEFAULT_NUM_ROWS);
    printf("  numRepeat: iterations to repeat (default %d\n", DEFAULT_NUM_REPEAT);
    printf("  density:   fraction of set cell values (default %f)\n", DEFAULT_DENSITY);
    exit(1);
}

#ifdef LIKWID
#include <likwid.h>
#else
#define likwid_markerInit()
#define likwid_markerThreadInit()
#define likwid_markerStartRegion(name)
#define likwid_markerStopRegion(name)
#define likwid_markerClose()
#endif // LIKWID

#if defined(USE_NATE)
// proceed row-by-row through lower left of squareMatrix
// for example, for numRows=4 the order is
//    ,0 ,1 ,2 ,3
//    ___________
// 0,| *  *  *  *
// 1,| 1  *  *  *
// 2,| 2  3  *  *
// 3,| 4  5  6  *
// [1,0] [2,0] [2,1] [3,0] [3,1] [3,2]
//
// Since C is "row major" the order through memory is
// 4 : 8, 9 : 12, 13, 14
//
// The corresponding inputs are
// [1,0] [0,2] [1,2] [0,3] [1,3] [2,3]
// with memory order
// 1 : 2, 6 : 3, 7, 11

void convert(uint32_t *squareMatrix, uint32_t numRows) {
    for (uint32_t rowNum = 1; rowNum < numRows; rowNum++) {
        uint32_t *outPtr = squareMatrix + rowNum * numRows; // [rowNum, 0]
        uint32_t *inPtr = squareMatrix + rowNum;            // [0, rowNum]
        for (uint32_t colNum = 0; colNum < rowNum; colNum++) {
            *outPtr = *inPtr;         // copy value from inPtr to outPtr
            outPtr = outPtr + 1;      // advance outPtr across the row
            inPtr = inPtr + numRows;  // advance inPtr down the column
        }
    }
}
#elif defined(USE_KARIM)
/*
 * Convert 2D Adjacency Matrix “G” from directed to undirected graph
 *
 * dim: array dimension
 */
void convert(uint32_t *G, long dim)
{
    unsigned long rowMajorArrSize;
    unsigned long i,h,t,r;
    unsigned long currentHeadRow, currentTailRow;
    unsigned long symmetricCell;
    unsigned long column;
    unsigned long dimMinusOne = dim-1;
    unsigned long headPointer, tailPointer;
    const unsigned long blockSize = 8;
    /* assign HEAD and TAIL pointers */
    headPointer = 0;
    rowMajorArrSize = dim*dim;
    // note: dim has to be long or do casting
    tailPointer = (rowMajorArrSize)-blockSize;

    unsigned long blocksLength = (rowMajorArrSize/blockSize)/2;
    for (i = 0; i <blocksLength; i++)
        {
            // HEAD block
            for (h = 0; h < blockSize ; h++)
                {
                    currentHeadRow = headPointer/dim;
                    column = headPointer-currentHeadRow*dim;
                    symmetricCell = headPointer -(( (currentHeadRow-column)*(dimMinusOne)));;

                    if ( G[headPointer] ==1 ) G[symmetricCell] = 1;
                    headPointer++;
                }
            // TAIL block
            for (t = 0; t < blockSize; t++)
                {
                    currentTailRow = (tailPointer/dim);
                    column = tailPointer-currentTailRow*dim;
                    symmetricCell = tailPointer -(( (currentTailRow-column)*(dimMinusOne)));;
                    if ( G[tailPointer] ==1 ) G[symmetricCell] = 1;
                    tailPointer++;
                }
            // move TAIL pinter to the next symmetrical block with HEAD
            tailPointer = tailPointer - (blockSize*2);
        }
    // reversing the last operation since TAIL became smaller than HEAD
    tailPointer = tailPointer+blockSize;

    // REMAINING CELLS
    for (r = headPointer; r < tailPointer; r++)
        {
            currentHeadRow = headPointer/dim;
            column = headPointer-currentHeadRow*dim;
            symmetricCell = headPointer -(( (currentHeadRow-column)*(dimMinusOne)));;

            if ( G[headPointer] ==1 ) G[symmetricCell] = 1;
            headPointer++;
        }

}
#elif defined(USE_BASIC)
void convert(uint32_t *G, uint32_t dim) {
    uint64_t row, col;
    for (row = 0; row < dim; row++) {
        for (col = row + 1; col < dim; col++) {
            G[col*dim + row] = G[row*dim + col];
        }
    }
}
#elif defined(USE_CONDITIONAL)
void convert(uint32_t *G, uint32_t dim) {
    uint64_t row, col;
    for (row = 0; row < dim; row++) {
        for (col = row + 1; col < dim; col++) {
            uint32_t val = G[row*dim + col];
            if (val) G[col*dim + row] = 1;
        }
    }
}
#else
#error Define USE_NATE, USE_KARIM, USE_BASIC, or USE_CONDITIONAL to specify algorithm
#endif

int main(int argc, char **argv) {
    char *programName = argv[0];
    uint64_t numRows = DEFAULT_NUM_ROWS;
    if (argc > 1) numRows = atoi(argv[1]);
    if (numRows == 0) dieUsage(programName);

    uint64_t numRepeat = DEFAULT_NUM_REPEAT;
    if (argc > 2) numRepeat = atoi(argv[2]);
    if (numRepeat == 0) dieUsage(programName);

    double density = DEFAULT_DENSITY;
    if (argc > 3) density = atof(argv[3]);
    if (density == 0.0) dieUsage(programName);
    if (density > 1.0) dieUsage(programName);

    uint64_t numCols = numRows; // matrix is always square
    uint64_t numCells = numRows * numCols;
    uint32_t *squareMatrix = calloc(numCells, sizeof(uint32_t));
    if (squareMatrix == NULL) {
        double sizeGB = numRows * numCols * sizeof(uint32_t) /
            (double) (1000 * 1000 * 1000);
        printf("Unable to allocate a %ld x %ld matrix (%.1f GB)\n",
               numRows, numCols, sizeGB);
        exit(1);
    }

    // initialize upper right of matrix according to density
    for (uint64_t rowNum = 0; rowNum < numRows; rowNum++) {
        // start each row at first entry after the diagonal
        uint32_t *outPtr = squareMatrix + rowNum * numCols + rowNum + 1;
        for (uint64_t colNum = rowNum + 1; colNum < numRows; colNum++) {
            // and write a 1 or 0 until the end of the row
            float rand = drand48();
            uint32_t cellVal = 0;
            if (rand < density) cellVal = 1;
            *outPtr++ = cellVal;
        }
    }

    // convert the matrix to symmetric (safe to repeat)
    likwid_markerInit();
    likwid_markerThreadInit();
    likwid_markerStartRegion("convert");
    for (uint64_t i = 0; i < numRepeat; i++) {
        convert(squareMatrix, numRows);
    }
    likwid_markerStopRegion("convert");
    likwid_markerClose();

#if 0
    for (uint64_t rowNum = 0; rowNum < numRows; rowNum++) {
        for (uint64_t colNum = rowNum + 1; colNum < numRows; colNum++) {
            uint32_t in =  squareMatrix[colNum + rowNum * numCols];
            uint32_t out = squareMatrix[rowNum + colNum * numRows];
            if (in != out) {
                printf("Bad value at [%ld, %ld] / [%ld, %ld]\n",
                       rowNum, colNum, colNum, rowNum);
                exit(1);
            }
        }
    }
#endif

    free(squareMatrix);
    exit(0);
}