Re: GCC for Mac OS X

["Timothy J. Wood" <tjw at omnigroup dot com>]

>From: Stan Shebs <shebs@apple.com>
>David Young wrote:
>> 
>> As I understand it: the altivec patches to gcc do NOT perform automatic
>> vectorization (loop unrolling). The code above is NOT altivec ready.
>> 
>> Altivec defines additional assembly instructions that operate on an 
extended
>> data types called "vector" and "matrix", much like "float" is defined in
>> regular c. The operators are overloaded to perform the arithmetic, much 
like
>> a C++ library that defines vector types and overloads the obvious 
operators.
>
>That's right.  The code I'm looking at doesn't actually define a literal 
"matrix"
>type, instead it has "vector" and "pixel" types, and seems to allow things
>to be defined as "vector pixel".
>
>As you might imagine, the changes to support these new types and connect 
them
>to code generation are threaded all through GCC.


  Here is an example (that I was compiling with MrC a while back).  The 
support in MrC and gcc look just about identical as far as I can tell 
(having both been done by Motorola).

  The code below is definitely suboptimial in a lot of ways, but what it 
does show is that I was able to write some fairly straight-forward C code 
that controlled the instruction selection mechanism.  While I wouldn't 
discourage anyone from writing some experimental patches that could 
automatically vectorize code, it would be nice to still be able to control 
the instruction selection from C and have the compiler do mundane tasks like 
register assignment and scheduling.

  Since the support for Altivec in gcc is mostly broken in the version on 
Apple's ftp site, the code below only compiles with MrC.  I ended up 
rewriting the core functions in hand coded assembly.  It sure would have 
been nicer to use a higher level language (at least to get the general form 
of the assembly and then maybe perform some hand tweaking).

-tim




---- julia.c ----

//
// This computes a quaternion Julia set using the Altivec instructions on 
the G4
//

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

#define RADIUS (4.0)
#define COUNT  (100)

typedef struct _RGBColor {
    unsigned char r, g, b;
} RGBColor;


// There is a bug in the current compiler that requires we reverse these 
args
#define VEC_MADD(a, b, c) vec_madd(a, c, b)

static int debug = 0;

static void generateJulia(vector float c, unsigned int pixelSize, float 
unitSize, unsigned long *pixelIter
ations);
static vector float floatVector(float a, float b, float c, float d);
static float vectorFloatElement(vector float f, unsigned int i);
static unsigned long vectorLongElement(vector unsigned long v, unsigned int 
i);
static void floatVectorGet(vector float f, float *a, float *b, float *c, 
float *d);

#define COUNT_OPS
#ifdef COUNT_OPS
unsigned long long altivecOps = 0;
#endif

int main(int argc, char *argv[])
{
    vector float   c = (vector float)(0,0,0,0);
    unsigned int   size;
        RGBColor      *pixels;
        unsigned long *pixelIterations, pixelCount;
        FILE          *ppmFile;
        
        if (argc == 1) {
                size = 512;
                //c = (vector float)(-0.15652, -1.03225, 0.0, 0.0);
                c = (vector float)(-0.12, 0.74, 0.0, 0.0);
        } else if (argc != 6) {
                fprintf(stderr, "usage: %s size r i j k\n", argv[0]);
                exit(1);
        } else {
                float        *rijk;

                size = atoi(argv[1]);
        
                rijk = (float *)malloc(sizeof(float) * 4 * 2);
                rijk = (float *)((unsigned long)rijk & ~(16 - 1));
        
                rijk[0] = atof(argv[2]);
                rijk[1] = atof(argv[3]);
                rijk[2] = atof(argv[4]);
                rijk[3] = atof(argv[5]);
        
                c = vec_ld(0, rijk);
        }
        
        if (size % 4) {
            fprintf(stderr, "size must be a multiple of four.\n");
                exit(1);
        }
        
        pixelIterations = (unsigned long *)calloc(size * size, 
sizeof(unsigned long));
        generateJulia(c, size, 4.0, pixelIterations);

        pixels = (RGBColor *)malloc(size * size * sizeof(RGBColor));
        for (pixelCount = 0; pixelCount < size * size; pixelCount++) {
                // Plain gray-scale for now
                pixels[pixelCount].r = (unsigned char)((255 * 
pixelIterations[pixelCount]) / COUNT);
                pixels[pixelCount].g = (unsigned char)((255 * 
pixelIterations[pixelCount]) / COUNT);
                pixels[pixelCount].b = (unsigned char)((255 * 
pixelIterations[pixelCount]) / COUNT);
        }
        
        // Write out a ppm file
        ppmFile = fopen("julia.ppm", "w");
        if (!ppmFile) {
            fprintf(stderr, "Cannot open julia.ppm\n");
                exit(1);
        }
        fprintf(ppmFile, "P6\012%d %d\012255\012", size, size);
        fwrite(pixels, sizeof(RGBColor), size * size, ppmFile);
        
#ifdef COUNT_OPS
    fprintf(stderr, "num ops = %d:%d\n", (unsigned int)(altivecOps >> 32), 
(unsigned int) altivecOps);
#endif

        return 0;
}

static vector bool long selectMask[4] = {
        (vector bool long)(0xffffffff, 0x00000000, 0x00000000, 0x00000000),
        (vector bool long)(0x00000000, 0xffffffff, 0x00000000, 0x00000000),
        (vector bool long)(0x00000000, 0x00000000, 0xffffffff, 0x00000000),
        (vector bool long)(0x00000000, 0x00000000, 0x00000000, 0xffffffff)
};
        
static vector unsigned long zero    = (vector unsigned long)(0x00000000, 
0x00000000, 0x00000000, 0x0000000
0);

// Replace the current element at the indicated position with the next item 
to compute
static void insertElement(vector float *zrNext, vector float *ziNext, vector 
float *zjNext, vector float *
zkNext,
                          vector unsigned long *iteration,
                          vector float *zr, vector float *zi, vector float 
*zj, vector float *zk,
                          unsigned int position)
{
    vector bool long sel;

    if (debug)
        fprintf(stderr, "  Loading element %d\n", position);

    // Shift the bitmask over the correct amount so that we'll
    // put the new element in the correct position
    sel = selectMask[position];

    // Fill in the z coordinates
    *zr = vec_sel(*zr, *zrNext, sel);
    *zi = vec_sel(*zi, *ziNext, sel);
    *zj = vec_sel(*zj, *zjNext, sel);
    *zk = vec_sel(*zk, *zkNext, sel);

    // Zero out the iteration counter for this field
    *iteration = vec_sel(*iteration, zero, sel);

#ifdef COUNT_OPS
    altivecOps += 5;
#endif
}
                                                  
                                                  
static void generateJulia(vector float c, unsigned int pixelSize, float 
unitSize, unsigned long *pixelIter
ations)
{
    vector float         zrNext, ziNext, zjNext, zkNext;
    vector float         zr, zi, zj, zk;
    vector float         cr, ci, cj, ck;
    vector float         xStep, yStep;
    vector unsigned long iteration = (vector unsigned long)(0, 0, 0, 0);
    vector unsigned long oneI = (vector unsigned long)(1, 1, 1, 1);
    vector float         radius = (vector float)(RADIUS, RADIUS, RADIUS, 
RADIUS);
    vector unsigned long count  = (vector unsigned long)(COUNT, COUNT, 
COUNT, COUNT);
    unsigned int         pixelsUnfinished, xNext, yNext, pointX[4], 
pointY[4];

    yStep = floatVector(unitSize/pixelSize, unitSize/pixelSize, 
unitSize/pixelSize, unitSize/pixelSize);
    xStep = floatVector(unitSize/pixelSize, unitSize/pixelSize, 
unitSize/pixelSize, unitSize/pixelSize);

    zrNext = floatVector(-unitSize/2.0, -unitSize/2.0, -unitSize/2.0, 
-unitSize/2.0);
    ziNext = floatVector(unitSize/2.0, unitSize/2.0, unitSize/2.0, 
unitSize/2.0);
    zjNext = (vector float)(0.0, 0.0, 0.0, 0.0);
    zkNext = (vector float)(0.0, 0.0, 0.0, 0.0);

    cr = floatVector(vectorFloatElement(c, 0), vectorFloatElement(c, 0), 
vectorFloatElement(c, 0), vectorF
loatElement(c, 0));
    ci = floatVector(vectorFloatElement(c, 1), vectorFloatElement(c, 1), 
vectorFloatElement(c, 1), vectorF
loatElement(c, 1));
    cj = floatVector(vectorFloatElement(c, 2), vectorFloatElement(c, 2), 
vectorFloatElement(c, 2), vectorF
loatElement(c, 2));
    ck = floatVector(vectorFloatElement(c, 3), vectorFloatElement(c, 3), 
vectorFloatElement(c, 3), vectorF
loatElement(c, 3));

    // Load the first four positions
    insertElement(&zrNext, &ziNext, &zjNext, &zkNext, &iteration, &zr, &zi, 
&zj, &zk, 0);
    zrNext = vec_add(zrNext, xStep);
    pointX[0] = 0; pointY[0] = 0;

    insertElement(&zrNext, &ziNext, &zjNext, &zkNext, &iteration, &zr, &zi, 
&zj, &zk, 1);
    zrNext = vec_add(zrNext, xStep);
    pointX[1] = 1; pointY[1] = 0;

    insertElement(&zrNext, &ziNext, &zjNext, &zkNext, &iteration, &zr, &zi, 
&zj, &zk, 2);
    zrNext = vec_add(zrNext, xStep);
    pointX[2] = 2; pointY[2] = 0;

    insertElement(&zrNext, &ziNext, &zjNext, &zkNext, &iteration, &zr, &zi, 
&zj, &zk, 3);
    zrNext = vec_add(zrNext, xStep);
    pointX[3] = 3; pointY[3] = 0;

    pixelsUnfinished = pixelSize * pixelSize;
    xNext = 4; yNext = 0;
    while (pixelsUnfinished) {
        vector float ijk2Sum, zr2, zi2, zj2, zk2;
        vector float zero = (vector float)(0.0, 0.0, 0.0, 0.0);
        vector float two  = (vector float)(2.0, 2.0, 2.0, 2.0);
        vector float mags = zero;
        vector bool long doneCountMask, doneRadiusMask, doneMask, checkMask;
        unsigned int position;

        // Process this set of points until one of them is done.  Use a 
do-while
        // loop to avoid having to clear 'mags' when points are added to the
        // active set.
        do {
            iteration = vec_add(iteration, oneI);
#ifdef COUNT_OPS
            altivecOps += 1;
#endif

            if (debug) {        
                for (position = 0; position < 4; position++)
                    fprintf(stderr, "  z[%d] = %f %f %f %f, iteration %d, 
mags = %f\n",
                            position,
                            vectorFloatElement(zr, position),
                            vectorFloatElement(zi, position),
                            vectorFloatElement(zj, position),
                            vectorFloatElement(zk, position),
                            vectorLongElement(iteration, position),
                            vectorFloatElement(mags, position));

                fprintf(stderr, "zr = %f %f %f %f\n",
                        vectorFloatElement(zr, 0),
                        vectorFloatElement(zr, 1),
                        vectorFloatElement(zr, 2),
                        vectorFloatElement(zr, 3));
            }
            
            zr2 = VEC_MADD(zr, zr, zero);
            zi2 = VEC_MADD(zi, zi, zero);
            zj2 = VEC_MADD(zj, zj, zero);
            zk2 = VEC_MADD(zk, zk, zero);
#ifdef COUNT_OPS
            altivecOps += 4;
#endif

            if (debug) {
                fprintf(stderr, "zr2 = %f %f %f %f\n",
                        vectorFloatElement(zr2, 0),
                        vectorFloatElement(zr2, 1),
                        vectorFloatElement(zr2, 2),
                        vectorFloatElement(zr2, 3));

                for (position = 0; position < 4; position++)
                    fprintf(stderr, "  z2[%d] = %f %f %f %f,\n",
                            position,
                            vectorFloatElement(zr2, position),
                            vectorFloatElement(zi2, position),
                            vectorFloatElement(zj2, position),
                            vectorFloatElement(zk2, position));
                }

            ijk2Sum = vec_add(vec_add(zi2, zj2), zk2);
#ifdef COUNT_OPS
            altivecOps += 2;
#endif

            zk = VEC_MADD(VEC_MADD(two, zr, zero), zk, ck);
            zj = VEC_MADD(VEC_MADD(two, zr, zero), zj, cj);
            zi = VEC_MADD(VEC_MADD(two, zr, zero), zi, ci);
            zr = vec_add(vec_sub(zr2, ijk2Sum), cr);
#ifdef COUNT_OPS
            altivecOps += 8;
#endif

            mags = vec_add(zr2, ijk2Sum);
#ifdef COUNT_OPS
            altivecOps += 1;
#endif



#ifdef COUNT_OPS
            // they won't always both be used, but usually...
            altivecOps += 2;
#endif
            } while (vec_all_le(iteration, count) && vec_all_lt(mags, 
radius));

        if (debug) {
            for (position = 0; position < 4; position++)
                fprintf(stderr, "  z[%d] = %f %f %f %f, iteration %d, mags = 
%f\n",
                        position,
                        vectorFloatElement(zr, position),
                        vectorFloatElement(zi, position),
                        vectorFloatElement(zj, position),
                        vectorFloatElement(zk, position),
                        vectorLongElement(iteration, position),
                        vectorFloatElement(mags, position));
            }

        // Build a mask of which points are done -- there could be more than 
one!
        doneCountMask  = vec_cmpgt(iteration, count);
        doneRadiusMask = vec_cmpge(mags, radius);
        doneMask       = vec_or(doneCountMask, doneRadiusMask);
#ifdef COUNT_OPS
        altivecOps += 3;
#endif

        if (debug) {
            fprintf(stderr, "  doneMask = 0x%08x 0x%08x 0x%08x 0x%08x\n",
                    vectorLongElement((vector unsigned long)doneMask, 0),
                    vectorLongElement((vector unsigned long)doneMask, 1),
                    vectorLongElement((vector unsigned long)doneMask, 2),
                    vectorLongElement((vector unsigned long)doneMask, 3));
            }

        // Determine which fields are done
        for (position = 0; position < 4; position++) {
            checkMask = vec_and(doneMask, selectMask[position]);

            // This should either be zero or equal to selectMask[position]
#ifdef COUNT_OPS
            altivecOps += 2;
#endif
            if (vec_all_eq((vector unsigned long)checkMask, (vector unsigned 
long)selectMask[position])) {
                if (debug)
                    fprintf(stderr, "  Position %d done\n", position);

                // We are done with a pixel
                pixelsUnfinished--;

                // Set the pixel
                pixelIterations[pointX[position] + pointY[position] * 
pixelSize] = vectorLongElement(itera
tion, position);

                // Load the next value
                insertElement(&zrNext, &ziNext, &zjNext, &zkNext, 
&iteration, &zr, &zi, &zj, &zk, position
);
                pointX[position] = xNext;
                pointY[position] = yNext;
                xNext++;
                if (xNext >= pixelSize) {
                    xNext = 0;
                    yNext++;
                    if (yNext >= pixelSize) {
                        // all done!
                        return;
                        }
                    zrNext = floatVector(-unitSize/2.0, -unitSize/2.0, 
-unitSize/2.0, -unitSize/2.0);
                    ziNext = vec_sub(ziNext, yStep); // We are stepping down 
on the y axis
#ifdef COUNT_OPS
                    altivecOps += 2;
#endif
                    fprintf(stderr, "Starting row %d\n", yNext);
                    } else {
                        zrNext = vec_add(zrNext, xStep);
#ifdef COUNT_OPS
                        altivecOps += 1;
#endif
                        }
                }
            }

        if (debug) {
            if (getchar() == 'q')
                return;
            }
        }
}

static vector float floatVector(float a, float b, float c, float d)
{
    float *f;

    f = (float *)alloca(sizeof(float) * 8);
    f = (float *)((unsigned long)f & ~(16 - 1));

    f[0] = a;
    f[1] = b;
    f[2] = c;
    f[3] = d;

    return vec_ld(0, f);
}

static void floatVectorGet(vector float f, float *a, float *b, float *c, 
float *d)
{
    float *fp;

    fp = (float *)alloca(sizeof(float) * 8);
    fp = (float *)((unsigned long)fp & ~(16 - 1));

    vec_st(f, 0, fp);
    *a = fp[0];
    *b = fp[0];
    *c = fp[0];
    *d = fp[0];
}

static float vectorFloatElement(vector float f, unsigned int i)
{
    float *fp;

    fp = (float *)alloca(sizeof(float) * 8);
    fp = (float *)((unsigned long)fp & ~(16 - 1));

    vec_st(f, 0, fp);
    return fp[i];
}

static unsigned long vectorLongElement(vector unsigned long v, unsigned int 
i)
{
    unsigned long *vp;

    vp = (unsigned long *)alloca(sizeof(unsigned long) * 8);
    vp = (unsigned long *)((unsigned long)vp & ~(16 - 1));

    vec_st((vector float)v, 0, (float *)vp);
    return vp[i];
}