root/trunk/blade/CodegenX86.d

//  Written in the D programming language 1.0
/**
* BLADE Alpha -- Basic Linear Algebra D Expressions
*
* Generate near-optimal x87/SSE/SSE2 asm code for BLAS1 basic vector operations
* at compile time.
* 32, 64 and 80 bit vectors are all supported.
* Uses techniques described in Agner Fog's superb Pentium optimisation manual (www.agner.org).
*
* Author:
*   Don Clugston.
* License:
*   Public domain.
*
* FEATURES:
*  - Supports any mix of vector addition, subtraction, dot product, and multiplication
*    by a scalar, with strided vector access.
*  - Generates either x87, SSE, or SSE2 asm code.
*  - If x87 code is generated, 80-bit precision is used whenever possible.
*  - Supports mixed-length operations (eg, real[] + double[] + float[]).
*
* BUGS/ FUTURE DIRECTIONS:
*  None of these support matrix operations.
* X87:
*  - Not optimal for the case of multiple real vectors (they could share a counter).
*  - Not optimal for the case where all vectors are 80-bit (two counters are used, but only is required).
*  - Doesn't take full advantage of length being known at compile time (loop unrolling
*     is possible).
*  - Doesn't use EBP register -- this would allow an extra vector in expressions.
*   (to do this, need naked asm with no stack frame).
* SSE/SSE2:
*  - SSE functions don't support unaligned data. Need to generate seperate code
*    for that case (NOTE: probably only worth doing for small expressions).
*
* THEORY:
* An expression string of the form "A+(B*C)*D" is used, together with a tuple, the
* entries of which correspond to A, B, C, ...
* This string is converted to postfix. The postfix string is converted to
* a string containing x87 asm, which is then mixed into a function which accepts the tuple.
*/

/*
 * POTENTIAL FROM RECENT INSTRUCTION SETS:
 * SSE5(AMD): fmaddpd can dramatically improve both performance and accuracy.
 * SSE4(Intel): dppd has some limited use.
 */

module blade.CodegenX86;
private import blade.BladeUtil;
private import blade.PostfixX86;

private:

// --------------
// Ranklist functions

// Count the number of vectors
int countVectors(char[] ranklist)
{
    int numVecs=0;
    for (int i=0; i<ranklist.length; ++i) {
        if (ranklist[i]=='1') ++numVecs;
    }
    return numVecs;
}

int vectorNum(char [] ranklist, char var)
{
    int numVecs=0;
    for (int i=0; i<var-'A'; ++i) {
        if (ranklist[i]=='1') ++numVecs;
    }
    return numVecs;
}

int strideVectorNum(char [] ranklist, char [] stridelist, char var)
{
    int numVecs=0;
    for (int i=0; i<var-'A'; ++i) {
        if (ranklist[i]=='1' && stridelist[i]=='1') ++numVecs;
    }
    return numVecs;
}

int scalarNum(char [] ranklist, char var)
{
    int k=0;
    for (int i=0; i<var-'A'; ++i) {
        if (ranklist[i]=='0') ++k;
    }
    return k;
}

int realScalarNum(char [][] typelist, char [] ranklist, char var)
{
    int k=0;
    for (int i=0; i<var-'A'; ++i) {
        if (ranklist[i]=='0' && typelist[i]=="real") ++k;
    }
    return k;
}
private:
// -------------------------------
//   Mixins to generate x87 ASM code
// -------------------------------

/// True if the character is an operation (everything else is an operand)
bool isInstruction(char op)
{
    return (op=='+' || op=='*' || op=='-'|| op=='_' || op=='=');
}

/// Count the number of temporaries which occur in the postfix expression.
int countTemporaries(char [] postfix)
{
// A temporary occurs whenever we load two values without an operation performed on the
// first one.
    int numTemps=0;
    for (int i=1; i<postfix.length; ++i) {
        if (!isInstruction(postfix[i-1]) && !isInstruction(postfix[i])) numTemps++;
    }
    return numTemps;
}


/// The maximum number of simultaneous temporary values in the postfix expression.
int maxActiveTemporaries(char [] postfix)
{
    int maxTemps=0;
    int numTemps=0;
    for (int i=1; i<postfix.length; ++i) {
        if (!isInstruction(postfix[i-1]) && !isInstruction(postfix[i])) numTemps++;
        if (isInstruction(postfix[i-1]) && isInstruction(postfix[i])) numTemps--;
        if (maxTemps<numTemps) maxTemps=numTemps;
    }
    return maxTemps;

}

unittest {
    assert(countTemporaries("AB*BC*+DE*+")==3);
    assert(maxActiveTemporaries("AB*BC*+DE*+")==2);
}

char [] operandSize(char [] typestr)
{
    switch(typestr) {
        case "real":   return "real ptr ";
        case "double": return "double ptr ";
        case "float":  return "float ptr ";
        default:
        assert(0, typestr);
    }
}

char [][char] opToX87() {
    return ['*':"fmul"[], '+': "fadd", '-': "fsub", '_': "fsubr"]; }

char [][char] opToSSE() {
    return ['*':"mulp"[], '+': "addp", '-': "subp", '/': "divp"]; }

char [][char] opToSSESingle() {
    return ['*':"muls"[], '+': "adds", '-': "subs", '/': "divs"]; }

static if (real.sizeof==10)      const char [] REALSIZE = "10";
else static if (real.sizeof==12) const char [] REALSIZE = "12";
else static if (real.sizeof==16) const char [] REALSIZE = "16";

char [] vectorSize(char [] typestr)
{
    switch (typestr) {
        case "double": return "8";
        case "float": return "4";
        case "real": return REALSIZE;
    }
}


// First, use the scratch registers (EAX, ECX, EDX). EAX is always used as
// an index register. If there are more than 2 vectors, use EBX, ESI, and EDI,
// which need to be pushed and popped.
// TODO: Finally, use the frame register EBP.
const char [][5] vectorRegister = ["ECX", "EDX", "EBX", "ESI", "EDI"];

public:

// Maximum number of vectors allowable in a expression for this code generator
const int MAX_X87_VECTORS = vectorRegister.length;
const int MAX_SSE_VECTORS = vectorRegister.length;
// Maximum number of real scalars allowable in an expression (
// (max # temporaries + max # real scalars) must be <=8, otherwise FPU stack
// will overflow).
const int MAX_87_REALSCALARSPLUSTEMPORARIES = 8;

private:

// Create code to push all used vector registors.
char [] pushRegisters(int numVectors)
{
    char [] result = "";
    for (int i=2; i<numVectors; ++i) result~= " push " ~ vectorRegister[i] ~ ";";
    return result ~ \n;
}

// Create code to pop all used vector registors.
char [] popRegisters(int numVectors)
{
    char [] result = ";  ";
    for (int i=numVectors-1; i>=2; --i) result~= "pop " ~ vectorRegister[i] ~ "; ";
    return result ~ \n;
}

// indexed by i.
char [] indexedVector(char [][] typelist, char [] ranklist, char [] stridelist, char var)
{
    if (typelist[var-'A']=="real") return " real ptr [" ~ vectorRegister[vectorNum(ranklist, var)] ~ "]";
    else if (stridelist[var-'A']=='1') return operandSize(typelist[var-'A']) ~ "[" ~ vectorRegister[vectorNum(ranklist, var)] ~ "]";
    return operandSize(typelist[var-'A']) ~ "[" ~
            vectorRegister[vectorNum(ranklist, var)] ~ " + " ~ vectorSize(typelist[var-'A']) ~ "*EAX]";
}

// indexed by i-1
char [] indexedVectorPrev(char [][] typelist, char [] ranklist, char var)
{
    char [] stride = " - " ~ vectorSize(typelist[var-'A']);
    if (typelist[var-'A'] == "real") return " real ptr [" ~ vectorRegister[vectorNum(ranklist, var)] ~ stride ~ "]";
    return operandSize(typelist[var-'A']) ~ "[" ~
            vectorRegister[vectorNum(ranklist, var)] ~ " + " ~ vectorSize(typelist[var-'A']) ~ "*EAX" ~ stride ~ "]";
}

char [] indexedSSEVector(char [] ranklist, char var, char [] vecsize)
{
    return "[" ~ vectorRegister[vectorNum(ranklist, var)] ~ " + " ~ vecsize ~"*EAX]";
}

char [] indexedSSENext(char [] ranklist, char var, char [] vecsize)
{
    return "[" ~ vectorRegister[vectorNum(ranklist, var)] ~ " + " ~ vecsize ~"*EAX+16]";
}

char [] indexedVectorWithStride(char [][] typelist, char [] ranklist, char var, int stride)
{
    char [] stridestr = " - " ~ vectorSize(typelist[var-'A']) ~ "*" ~ itoa(stride);
    if (typelist[var-'A'] == "real") return " real ptr [" ~ vectorRegister[vectorNum(ranklist, var)] ~ stridestr ~ "]";
    return operandSize(typelist[var-'A']) ~ "[" ~
            vectorRegister[vectorNum(ranklist, var)] ~ " + " ~ vectorSize(typelist[var-'A']) ~ "*EAX" ~ stridestr ~ "]";
}

// Pop N values from the FPU stack
char [] discardFromStack(int n)
{
    char [] result="";
    while (n>1) {
        result~= "  fcompp ST(0), ST;"\n; // pop two values at once
        n-=2;
    }
    if (n==1) result~= "  fstp ST(0), ST;"\n;
    return result;
}

public:

/** Generate asm code which is optimal for x87 CPUs without SSE2.
 (Pentium, PMMX, PII, PIII). It is also optimal for recent x86 CPUs
 where vector sizes are mixed.

 There are two cases:
 (A) DAXPY-style loops, where every element is independent of the other indices;
 (B) DDOT-style loops, where the result for every element is accumulated.

 For cumulative loops, best performance is achieved with loop unrolling and
 multiple accumulators, in order to break dependency chains.

The key optimisation rules for DAXPY loops are:
 1. keep the loop overhead to one clock cycle if possible.
 2. (FMUL latency) don't use the result of a multiply immediately
Techniques to address these are:
 1. Use EAX as a counter and index variable, which begins negative and counts UP to zero.
    Combine counters for all packed doubles and floats into this single counter.
 2. The latency of fmul is avoided by swapping fadd/fsub with fmul whenever possible.

The generated code is of the form:
----
 load scalars onto FPU stack
 load vector pointers into EAX, EBX, ...
L1:
 calculate result into ST(0)
 increment pointers
 goto L1 if not done
 pop scalars off FPU stack
----

*/
char [] generateCodeForAsmX87(int numStrides, Values...)(char [] postfixOperations)
{
// Because of compiler bug #1125, no structs can be stored in the 'values' tuple.
// Thus, lengths and strides must be stored seperately from vector pointers.
    char [] ranklist;
    char [][] typelist;
    char [] stridelist; // "1" = strided, "0"=unstrided
    foreach(T; Values[0..$-numStrides]) {
        static if (is(typeof(T[0]))) {
            stridelist~="0";
            ranklist~="1";
            typelist ~= typeof(T[0]).stringof;
        } else static if (is(typeof(T.data))) {
            stridelist~="1";
            ranklist~="1";
            typelist ~= typeof(T.data[0]).stringof;
        } else {
            stridelist~="0";
            ranklist~="0";
            typelist ~= T.stringof;
        }
    }
    return generateCodeForAsmX87Impl(ranklist, typelist, stridelist, postfixOperations);
}

private:
// This is split off from the template to make code coverage easier.
char [] generateCodeForAsmX87Impl(char [] ranklist, char [][] typelist, char [] stridelist, char [] operations)
{
    char [] result="";
    char [] incrementRealVectors="";

    result ~= "// Operation : " ~  operations ~ \n;

    // Create local variables for pointers to vectors (avoid bug #1125)

    int vecnum = 0;
    int stridecount = 0;
    for (int i=0; i< ranklist.length;++i) {
        if (ranklist[i]=='1'){
            if (stridelist[i]=='1') {
                incrementRealVectors ~= "  add " ~ vectorRegister[vecnum] ~ ", values[" ~ itoa(ranklist.length+stridecount) ~ "];\n";
                ++stridecount;
            } else if (typelist[i]=="real") {
                incrementRealVectors ~= "  add " ~ vectorRegister[vecnum] ~ ", " ~ REALSIZE ~ ";\n";
            }
            ++vecnum;
        }
    }

    int numScalarsOnStack=0;

    result~= \n"asm {"\n ~ pushRegisters(vecnum);
    // EAX will be the counter
    result ~= "  mov EAX, veclength;"\n;

    // Load all the vector pointers into registers, and push all the scalars onto the stack

    int numvecs=0;
    int numconsts=0;
    for (int i=0; i<ranklist.length; ++i) {
      if (ranklist[i]=='1') {
          if (typelist[i]=="real" || stridelist[i]=='1') {
              result ~= "  mov " ~ vectorRegister[numvecs] ~ ", values[" ~ itoa(i) ~ "];";
          } else  {
            result ~= "  lea " ~ vectorRegister[numvecs]
              ~ ", [" ~ vectorSize(typelist[i]) ~ "*EAX];   "
              ~ "  add " ~ vectorRegister[numvecs] ~ ", values[" ~ itoa(i) ~ "];";
         }
         result ~= "  //" ~ cast(char)('A'+i) ~ \n;
        ++numvecs;
      } else if (typelist[i]=="real") {
          result ~= "  fld real ptr values["~ itoa(i) ~"];";
          ++numconsts;
          ++numScalarsOnStack;
         result ~= "  //" ~ cast(char)('A'+i) ~ \n;
      }
    }
    int done=0;

    // We need to keep track of how many things are on the FPU stack.
    // Every time something is pushed, the indices of our variables change!
    int numOnStack = 0; // How much of the FP stack is being used?

    bool isCumulative = (operations[0]=='0' || operations[0]=='1');
    if (operations[0]=='0') {
        result ~= "  fldz;"\n; // dot product
        ++numOnStack;
        done = 1;
    } else if (operations[0]=='1') {
        result ~= "  fld1;"\n; // prod
        ++numOnStack;
        done = 1;
    }
    result ~= "  xor EAX, EAX; "\n
        "  sub EAX, veclength; // counter=-length"\n
        "  jz short L3; // test for length==0"\n;

    // Construct the main body of the loop (the main body does not include
    // the final storage instruction, because of the FST latency).
    char [] mainbody = "";

    while(done<operations.length) {
        char [] next;
        if (isInstruction(operations[done])) {
            // Perform an arithmetic operation on the top two FPU stack items.
            next = "  " ~ opToX87[operations[done]] ~ "p ST(1), ST;  //" ~ operations[done] ~ \n;
            mainbody ~= next;
            ++done;
            numOnStack--;
        } else if (operations[done]=='a') {
            mainbody ~= "  fabs;"\n;
            ++done;
        } else if (operations[done]=='n') {
            mainbody ~= "  fchs;"\n;
            ++done;
        } else if (operations[done]=='q') {
            mainbody ~= "  fsqrt;"\n;
            ++done;
        } else if (!isInstruction(operations[done+1])){
            // load a vector onto the FPU stack, to begin a new subexpression.
            int u  = operations[done]-'A';
            if (ranklist[operations[done]-'A']=='1') {
                next = "  fld "  ~ indexedVector(typelist, ranklist, stridelist, operations[done] ) ~ ";  //" ~ operations[done] ~\n;
            } else { // load constant. Will never be a real
                next = "  fld " ~ operandSize(typelist[operations[done]-'A']) ~ "values[" ~ itoa(operations[done]-'A') ~"]; // * " ~ operations[done..done+1] ~ "\n";
            }
            mainbody ~= next;
            ++done;
            numOnStack++;
        } else if (operations[done]==',') {
            mainbody ~= "  " ~ opToX87[operations[done+1]] ~ " ST, ST(0);    // dup " ~ operations[done+1] ~ \n;
            done+=2;
        } else if (ranklist[operations[done]-'A']=='1') {
             // An operation will be performed between the stack top and a vector.
             // If it's a float or double, we can combine the load+arithmetic op
             // into a single instruction. Stores of reals can also be done in one instr.
            char [] comment = ";  // " ~ operations[done..done+2] ~ \n;
            if (operations[done+1]=='=') {
                // If it's the last operation, pop it from the stack; otherwise,
                // it chains.
                next = ((done+2 == operations.length)? "  fstp " : "  fst ")
                    ~ indexedVector(typelist, ranklist, stridelist, operations[$-2] ) ~ comment;
            } else if (typelist[operations[done]-'A']=="real") {
                 // 80-bit vectors must be loaded onto the FPU stack first
                next = "  fld real ptr ["  ~ vectorRegister[vectorNum(ranklist, operations[done])] ~ "]; //" ~ operations[done] ~ \n
                    ~ "  " ~ opToX87[operations[done+1]] ~ "p ST(1), ST; //" ~ operations[done+1] ~\n;
            } else { // floats and doubles can be used directly
                next = "  " ~ opToX87[operations[done+1]] ~ " "
                  ~ indexedVector(typelist, ranklist, stridelist, operations[done] ) ~ comment;
            }
            mainbody ~= next;
            done+=2;
        } else { // multiply by scalar.
          if (typelist[operations[done]-'A']=="real") {
            // Multiply by real scalar, which is already on the stack.
            next = "  fmul ST, ST(" ~ itoa(numOnStack + numScalarsOnStack - realScalarNum(typelist, ranklist, operations[done]-'A')-1) ~ "); // * " ~ operations[done] ~ \n;
            mainbody ~= next;
          } else {
            // For scalar float or double values, we can multiply directly, saving one slot on the FP stack.
            next = "  fmul " ~ operandSize(typelist[operations[done]-'A']) ~ "values[" ~ itoa(operations[done]-'A') ~"]; // * " ~ operations[done..done+1] ~ "\n";
            mainbody ~= next; //firstbody ~= next;
          }
            done +=2;
        }
    }

    result ~= \n
        ~ "  align 4;\n"
        ~ "L1:\n" ~ mainbody;

//    if (cumulatingOp) result ~= "  " ~ opToX87[cumulatingOp] ~ "p ST(2), ST;"\n;

    result ~= incrementRealVectors // Update the counters
           ~ "  inc EAX;\n  jnz L1;\n";

    // Discard any scalars that are left on the stack
    if (isCumulative && numScalarsOnStack>0) {
        // Preserve the result of the dot product
        result ~= "  fxch ST(" ~ itoa(numScalarsOnStack) ~ "), ST;"\n;
    }
    result ~= discardFromStack(numScalarsOnStack);

    result~= "L3:" \n ~ popRegisters(vecnum) ~ "}\r\n";

    return result;
}

//-----------------------------

char [] XMM(int k) { return "XMM"~ itoa(k); }

public:

/** Generate BLAS1 asm code using SSE or SSE2.
 * For SSE2, all scalars are double, vectors are double*; for SSE1, all are float.
 * At entry, all vector parameters are aligned.
 */
char [] generateCodeForSSE(Values...)(char [] operations)
{
    char [] ranklist;
    bool usingDoubles=false;
    foreach(T; Values) {
        static if (is(typeof(T[0]))) ranklist~="1"; else ranklist~="0";
        static if (is(T == double) || is(T == double *)) {
            usingDoubles = true;
        } //else assert(is(T==float)|| is(T==float*));
    }
    return generateCodeForSSEImpl(usingDoubles, ranklist, operations);
//    makePostfixForSSE(infixOperations, ranklist));
}

// Note: If SSE4 is available, the dppd and dpps instructions could be
// used to replace the final *+ in dot-product operations. This would allow
// an in-order dot product to be performed, but otherwise doesn't give much speed
// improvement.

private:

// split off from the template to make code coverage work
char [] generateCodeForSSEImpl(bool usingDoubles, char [] ranklist, char [] operations, char cumulatingOp=0)
{
    char [] result="";

    result ~= "// Operation : " ~  operations ~ \n;

    int numvecs = countVectors(ranklist);
    int numScalarsOnStack=0;
    bool isCumulative = (operations[0]=='0') || (operations[0]=='1');
    if (isCumulative) result ~= (usingDoubles? "  double" : "  float") ~" sum;"\n;

    result~= \n"asm {"\n ~ pushRegisters(numvecs);
    // EAX will be the counter
    result ~= "  mov EAX, veclength;"\n;
    // Load all the vector pointers into registers

    char [] vectorsize = usingDoubles? "8" :"4"; // size of a double
    char [] suffix = usingDoubles? "d " :"s ";

    int vecregnum = 0;
    int numconsts=0;
    for (int i=0; i<ranklist.length; ++i) {
      if (ranklist[i]=='1') {
        result ~= "  lea " ~ vectorRegister[vecregnum]
          ~ ", [" ~ vectorsize ~ "*EAX];   "