root/trunk/blade/Blade.d

//  Written in the D programming language 1.0
/**
* BLADE 0.6Alpha -- Basic Linear Algebra D Expressions
*
* Generate near-optimal x87/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, unary minus,
*    multiplication by a scalar,
*    cumulation via dot product, sum() and prod(), and multidimensional slicing.
*  - Generates either x87 asm code, SSE or SSE2 asm code or pure D, depending on
*    the complexity of the expression, and the availability of inline asm.
*  - When static arrays are used, mismatches in array length are detected
*    at compile time.
*  - Error messages refer to the line of user code which generated the error.
*    The library never produces a torrent of undecipherable template error messages.
*  - A 'const folding' (actually vector/scalar folding) step is performed.
*
* SPEED/ACCURACY TRADEOFF:
* Tradeoff arises because IEEE floating point multiplication and addition are not associative.
*  Assuming that overflow and underflow do not occur:
*  (a*b)*c may differ from a*(b*c) in the last bit.
*  (a+b)+c may differ from a+(b+c) by a factor a million or more.
*
* - Multiplication is assumed to be associative.
* - Addition and subtraction are not treated as associative.
*    Except: Addition inside a dot product or vector sum is treated as associative.
*
* FUTURE DIRECTIONS (in order of expected implementation):
* - nested D expressions
* - cumulative operations min, max
* - Loop unrolling for cumulative operations dot, sum, prod.
* - Dense matrix support.
* - Triangular, banded, symmetric, and sparse matrix support
*
* 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 asm instructions, which is then mixed into a function
* which accepts the tuple.
*
* COMPILER BUGS/LIMITATIONS AFFECTING THIS LIBRARY
* - Local arrays are not aligned to a 128-bit boundary, so use of aligned SSE is not
*   always possible.
* - Bugzilla #1125 -- structs in a tuple can't be used in asm.
* - Bugzilla #1382 -- CTFE strings never get deleted --> SLOOOOOW compilation. KILLER BUG.
* - Bugzilla #1768 -- in CTFE, arrays of arrays aren't initialized properly
*
* HISTORY:
* 0.1 - Used classes to make expression templates.
* 0.2 - Support for a wider variety of expressions. Dot product, imaginary numbers, etc.
* 0.3 - Based on string mixins. Most of the new features of 0.2 are gone, but SSE2 is added.
* 0.4 - Added D code generator. Nice error messages. Optimal parameter passing.
*       (passes pointers, not arrays).
* 0.5 - Expression simplification step. Slicing support.
* 0.6 - Dot product, nested expressions (asm only), intrinsics: abs, sqrt, sum.
* 0.7 - Intrinsics: prod
*/

module blade.Blade;

public import blade.SyntaxTree : AbstractSyntaxTree, mixin_tupleAndSyntaxtreeof, AST, Symbol;
private import blade.BladeUtil : enquote, itoa;
private import blade.BladeRank : isStrided, exprRank;
private import blade.BladeSimplify : simplifySyntaxTree, RevisedExpression, remapCompounds;
private import blade.CodegenX86 : generateCodeForAsmX87, MAX_X87_VECTORS,
                                 MAX_87_REALSCALARSPLUSTEMPORARIES,
                                 generateCodeForSSE,  MAX_SSE_VECTORS;
private import blade.BladeVisitor: expressionContainsAssignment;

private import blade.PostfixX86 : makePostfixForX87, makePostfixForSSE;

public:

// FOR MIXIN: Generate code to evaluate the given vector expression.
char [] vectorize(char [] expr)
{
    debug (BladeFrontEnd) {
    return `pragma(msg, \n ~ "// " __FILE__ ~ "(" ~__LINE__.stringof[0..$-1] ~ ") ` ~ enquote(expr) ~ `" ~ \n ~ ` ~ mixin_tupleAndSyntaxtreeof("makeVectorCode", expr) ~ "~\\n);"
    "mixin(" ~ mixin_tupleAndSyntaxtreeof("makeVectorCode", expr)~ ");";
    } else {
       return "mixin(" ~ mixin_tupleAndSyntaxtreeof("makeVectorCode", expr)~ ");";
    }
}

// Simplify the expression, categorise it,
// and dispatch to the appropriate code generator.
char [] makeVectorCode(Types...)(AbstractSyntaxTree tree)
{
    RevisedExpression revised = simplifySyntaxTree(tree);
    if (revised.errorMessage.length>0)  return `static assert(0, "BLADE: ` ~ enquote(revised.errorMessage) ~ `");`;
    VecExpressionType exprType = categorizeExpression(revised);
    InvocationCode q;
    if (exprType == VecExpressionType.SSE2Expression || exprType == VecExpressionType.SSE1Expression) {
        q = invokeSSE((exprType == VecExpressionType.SSE2Expression), revised);
    } else if (exprType == VecExpressionType.X87Expression) {
        q = invokeX87(revised);
    } else {
        q = DCodeGenerator(revised);
    }
    return q.assertions ~ q.invoker ~ ";";
}

// For a compound of a different dimensionality (eg a dot product), we may need
// to calculate the result seperately.
char [] makeVectorCodeForDimensionalCompound(char [] expression, AbstractSyntaxTree tree)
{
    // TODO:
    return expression;
}

template X87RetType(char [] expr) {
    static if (expr[0]!='0' && expr[0]!='1') alias void X87RetType;
    else alias real X87RetType;
}

// These functions have the complete expression encoded in the template type.
// One of these functions is instantiated for each expression.
// A difficulty is, that the only way to transfer information from the CTFE code
// into the function, is via the template parameters. So from inside the function,
// we must re-assemble the type information, and use this to generate the asm code.

/** Function to implement BLAS1 operations using SSE/SSE2 assembler.
 * Every member of the Values tuple must only be double or double *.
 */
RetType SSEVECGEN(RetType, char [] expr, Values...)(int veclength, Values values) {
    debug(BladeBackEnd) {
       pragma(msg, generateCodeForSSE!(Values)(expr));
    }
    mixin(generateCodeForSSE!(Values)(expr));
}

/** Function to implement BLAS1 operations using X87 assembler.
 * Every member of the Values tuple must only be real,
 * float[], double [], or real[], or BladeStrided!(float), !(double), !(real)
 */
X87RetType!(expr) X87VECGEN(char [] expr, int numStrides, Values...)(int veclength, Values values) {
    debug(BladeBackEnd) {
        pragma(msg, generateCodeForAsmX87!(numStrides, Values)(expr));
    }
    mixin(generateCodeForAsmX87!(numStrides, Values)(expr));
}

private:
// Masks for setting or clearing the signbit in SSE registers.
static ulong[2] SSE_SIGNMASKpd = [0x7FFF_FFFF_FFFF_FFFFL, 0x7FFF_FFFF_FFFF_FFFFL];
static uint[4] SSE_SIGNMASKps = [0x7FFF_FFFF, 0x7FFF_FFFF, 0x7FFF_FFFF, 0x7FFF_FFFF];
static ulong[2] SSE_SIGNBITpd = [0x8000_0000_0000_0000L, 0x8000_0000_0000_0000L];
static uint[4] SSE_SIGNBITps = [0x8000_0000,0x8000_0000,0x8000_0000, 0x8000_0000];
// The value 1.0 for a parallel SSE register
static ulong[2] SSE_ONEpd = [0x3FF0_0000_0000_0000L, 0x3FF0_0000_0000_0000L];
static uint[4] SSE_ONEps = [0x3F0_000, 0x3F0_000, 0x3F0_000, 0x3F0_000];

private:

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

// ------------------------------------
// Categorize the expression type
// SSE2 is possible only if all vectors are doubles.
// SSE1 is possible only if all vectors are floats.
// X87 is possible for any mix of real, double, and float vectors.
// BUG: for X87, should also check number of temporaries (don't overflow the FPU stack)
enum VecExpressionType { SSE1Expression, SSE2Expression, X87Expression, DExpression };

VecExpressionType categorizeExpression(RevisedExpression tree)
{
    bool SSE2 = true;
    bool SSE1 = true;
    bool X87 = true;
    bool strided = false; // true if any strided vector or matrix operations exist
version (D_InlineAsm_X86) {} else {
    // Without an assembler, there's no chance!
    SSE2 = false;
    SSE1 = false;
    X87 = false;
}
int wholerank = exprRank(tree.expression, tree.rank);

if (wholerank==0 && expressionContainsAssignment(tree.expression)) {
    return VecExpressionType.DExpression; // Scalar assignments always use inline D.
}
    int numvectors = 0;
    int numscalars = 0;
    int numRealScalars = 0; // scalars other than float or double.
    for (int i=0; i<tree.mapping.length;++i) {
        char r = tree.rank[i];
        int x = tree.mapping[i]-'A';
        if (r=='0') {
            ++numscalars;
            if (x<tree.symbolTable.length) {
                char [] t = tree.symbolTable[x].element;
                if (t!="double" && t!="float" && t!="idouble" && t!="ifloat") {
                    ++numRealScalars;
                }
            }
            // TODO: disallow asm if any non-int/not FP types are used.
            continue;
        }
        if (r>'1') return VecExpressionType.DExpression; // can only do scalars and vectors right now.
        // At this point, all compounds are an original symbol + indexing/slicing.
        int y = x; // for compounds, get the original type
        if (x>=tree.symbolTable.length) {
            y = tree.compounds[x-tree.symbolTable.length][0]-'A';
            // Check for a stride..
            if (tree.compounds[x-tree.symbolTable.length][$-1]==']') {
                strided |= isStrided(tree.compounds[x-tree.symbolTable.length]);
            }
        }

        char [] t = tree.symbolTable[y].element;
        if (t == "double") {
            ++numvectors;
            SSE1 = false;
        } else if (t == "float") {
            ++numvectors;
            SSE2 = false;
        } else {
            SSE1 = false;
            SSE2 = false;
            if (t == "real") { ++numvectors; }
            else X87 = false;
        }
    }
    // It's not worth doing strided operations with SSE.
    if (strided) { SSE1=false; SSE2=false; }
    if (numRealScalars > MAX_87_REALSCALARSPLUSTEMPORARIES) X87 = false;
    if (numvectors > MAX_X87_VECTORS) X87 = false;
    if (numvectors > MAX_SSE_VECTORS) { SSE1=false; SSE2=false; }
    if (SSE1) return VecExpressionType.SSE1Expression;
    if (SSE2) return VecExpressionType.SSE2Expression;
    return X87 ? VecExpressionType.X87Expression : VecExpressionType.DExpression;
}


// This is mainly a workaround for compiler bug #1125. Ideally both
// pointer and stride would be stored together.
struct BladeStrided(T)
{
    T * data;     // Pointer to the first element
}

//-------------------------------------------------------
//                Invoker functions
//-------------------------------------------------------
// These are CTFE functions which, when mixed in, will call the BLAS template
// function. They ensure that all types are converted into standard
// simple forms, ensure that the vector lengths are equal, and pass in all
// of the parameters.

struct InvocationCode {
    char [] invoker; // For mixin: code to invoke the functions.
    char [] assertions; // For mixin: code to assert that everything is correct
}

/// Generate code which will call the X87 function.
InvocationCode invokeX87(RevisedExpression tree)
{
    char [] assertions = assertAllVectorLengthsEqual(tree);
    char [] result = "";
    char [] stridelist="";
    char [] alltypes="";

    char [][] typelist;

    char [] vals;
    int numstrides=0;
    for (int i=0; i<tree.mapping.length;++i) {
        char rnk = tree.rank[i];
        vals ~= ",";
        InvocationCode q = getValueForSymbol(tree.mapping[i], tree);
        char [] v = q.invoker;
        assertions ~= q.assertions;
        int x = tree.mapping[i]-'A';
        char [] t;
        bool strided = false;
        if (x<tree.symbolTable.length) {
            // it's an original symbol
            if (rnk=='0') {
                t = tree.symbolTable[x].type;
            } else t = tree.symbolTable[x].element;
        } else {
            // it's a compound
            if (rnk=='0') {
                t = "real"; // convert all compounds to real.
                // TODO (tricky): if the number is exactly representable as a double
                // or float, it could use less FPU stack space.
            } else { // for arrays, the type is the type of the original array
                t = tree.symbolTable[tree.compounds[x-tree.symbolTable.length][0]-'A'].element;
                        // Check for a stride..
                if (tree.compounds[x-tree.symbolTable.length][$-1]==']') {
                    strided = isStrided(tree.compounds[x-tree.symbolTable.length]);
                    if (strided) ++numstrides;
                }

            }
        }
        alltypes ~= ",";
        if (rnk=='0') {
            // Convert scalars into standard form.
            // long, ulong, and real must become real.
            // We convert everything else to double, since that uses less
            // FPU stack space.
            if (t == "real" || t == "double" || t == "float") alltypes ~= t;
            else if (t == "long" || t == "ulong") result ~= "real";
            else alltypes ~= "double"; // Convert all other scalars into doubles.
            vals ~= v;
        } else {
            assert (t == "real" || t=="double" || t=="float", "BLADE BUG");
            if (strided) {
                alltypes ~= "BladeStrided!(" ~ t ~ ")";
                vals ~= "BladeStrided!("~t~")(&" ~ v ~ "[0])";
                stridelist ~= "," ~ getStrideForSymbol(tree.mapping[i], tree);
            } else {
                alltypes ~= t ~ "*";
                // for vectors, we only need the pointer, not the length
                //vals ~= "&" ~  v ~ "[0]";
                vals ~= v ~ ".ptr";
            }
        }
        typelist ~= t;
    }
    char [] postfixops = makePostfixForX87(tree.expression, typelist, tree.rank);

    result ~= `X87VECGEN!("` ~ enquote(postfixops) ~ `"`;
    result ~= "," ~ itoa(numstrides);
    result ~= alltypes;
    for (int i=0; i<numstrides; ++i) result ~= ",int";
    result ~= ")(";
    int firstVector = findVectorForLength(tree);
    return InvocationCode(result ~ getValueForSymbol(tree.mapping[firstVector], tree).invoker ~ ".length"
        ~ vals ~ stridelist  ~ ")", assertions);
}

char [] generateAsserts(RevisedExpression tree, bool checkAlignment)
{
    char [] result = assertAllVectorLengthsEqual(tree);
    if (checkAlignment) result ~= assertAllVectorsAlign128(tree);
    return result;
}

/// Generate code which will call the SSE/SSE2 code generation function
InvocationCode invokeSSE(bool SSE2, RevisedExpression tree)
{
    char [] assertions = assertAllVectorLengthsEqual(tree)
    ~ assertAllVectorsAlign128(tree);

    char [] postfix = makePostfixForSSE(tree.expression, tree.rank);
    char [] retType = "void";
    if (postfix[0]=='0' || postfix[0]=='1') retType = (SSE2? "double" : "float");

    char [] result = "SSEVECGEN!(" ~ retType ~ `,"` ~ enquote(postfix) ~ `"`;
    // For SSE2, everything must be implicitly convertible to double.
    char [] vals;
    for (int i=0; i<tree.mapping.length;++i) {
        char rnk = tree.rank[i];
        if (rnk=='0') result ~= SSE2? ",double" : ",float";
        else result ~= SSE2? ",double*" : ",float*";
        vals ~= ",";
//        if (rnk=='1') vals ~= "&";
        InvocationCode q = getValueForSymbol(tree.mapping[i], tree);
        vals ~= q.invoker;
        assertions ~= q.assertions;
        // for vectors, we only need the pointer, not the length
        if (rnk=='1') vals ~= ".ptr";
    }

    result ~= ")(";
    int firstVector = findVectorForLength(tree);
    result ~= getValueForSymbol(tree.mapping[firstVector], tree).invoker ~ ".length";
    result ~= vals;

    return InvocationCode(result ~ ")", assertions);
}

/** Generates an assert which ensures that all vectors are of equal length.
 * If possible, the error will be detected at compile time.
 */
char [] assertAllVectorLengthsEqual(RevisedExpression tree)
{
    char [] result ="";
    int firstVector = findVectorForLength(tree);
//    bool known = arrayLengthIsStatic(tree.symbolTable[firstVector].type);
    for (int i=0; i<tree.mapping.length;++i) {
        if (tree.rank[i]=='1') {
            if (firstVector != i) {
//                if (known && arrayLengthIsStatic(tree.symbolTable[i].type)) {
//                    // both lengths are known at compile time - make it a
//                    // compile-time static assert
//                    result ~= "static ";
//                }
                result ~= "assert("
                 ~ getDimensionLengthForSymbol(tree.mapping[i], tree, 1)
                    ~ "==" ~ getDimensionLengthForSymbol(tree.mapping[firstVector], tree, 1)
                    ~ ", `Vector length mismatch`);"\n;
            }
        }
    }
    return result;
}

char [] assertAllVectorsAlign128(RevisedExpression tree)
{
    char [] result ="";
    for (int i=0; i<tree.mapping.length;++i) {
        if (tree.rank[i]=='1'){
            result ~= "assert( (cast(size_t)(" ~ getValueForSymbol(tree.mapping[i], tree).invoker
                    ~ ".ptr)& 0x0F) == 0, `SSE Vector misalignment: " ~ getValueForSymbol(tree.mapping[i], tree).invoker ~ "`);"\n;
        }
    }
    return result;
}

// Return true if the type has a length which is known at compile time
bool arrayLengthIsStatic(char [] type)
{
    if (type.length<3) return false;
    // BUG: this returns true for AA's.
    return type[$-1]==']' && type[$-2]!='[';
}

// Return a vector which contains the length of the expression.
// If possible, a vector with known (static) length will be chosen.
// If this is not possible, a normal dynamic array will be used.
// If all else fails, a sliced vector will be used.
int findVectorForLength(RevisedExpression tree)
{
    int dynamic = -1; // last dynamic vector
    int strided = 0; // last unstrided vector
    for (int i = 0; i < tree.mapping.length; ++i) {
        if (tree.rank[i]!='1') continue;
        int x = tree.mapping[i]-'A';
        strided = i;
        if (x < tree.symbolTable.length) {
            if (arrayLengthIsStatic(tree.symbolTable[x].type)) return i;
            dynamic = i;
        } else {
            // Check for a stride.
            if (tree.compounds[x-tree.symbolTable.length][$-1]==']') {
                if (!isStrided(tree.compounds[x-tree.symbolTable.length])) {
                    dynamic = i;
                }
            }
        }
    }
    return dynamic>=0? dynamic : strided;
}

bool hasDollar(char [] s)
{
    foreach(c; s) { if (c=='$') return true; }
    return false;
}

char [] getDimensionLengthForSymbol(char c, RevisedExpression tree, int dimension)
{
    int numSlicesRemaining = 1;
    assert(dimension == 1);
    char [] v = "";
    // is it an original symbol?
    if (c-'A'<tree.symbolTable.length) {
        v = tree.symbolTable[c-'A'].value;
        return v ~ ".length";
    } else {  // else it's a compound or an indexed array
        char [] comp = tree.compounds[c-'A'-tree.symbolTable.length];

        if (comp[$-1]!=']') { // simple compound expression
            foreach(d; comp) {
                if (d=='{') assert(0, "BLADE ICE");
                if (d>='A' && d<='Z') v ~= tree.symbolTable[d-'A'].value;
                else v ~= d;
            }
            return v ~ ".length";
        } else {
            // indexed array, possibly involving slicing
            int numbrack=0;
            bool hasSliced=false;
            // it's easier if we go backwards
            // Replace the last slice [a..b] operation with [a+firstIndexExpr]
            // or if no slices exist, append [firstIndexExpr] to the end.
            int numbracks = 0;
            bool isSlice = false;
            char [] nextIndex;
            char [] sliceTo;

            for (int k = comp.length-1;k>=1; --k) {
                char d = comp[k];
                if (d == ']') { ++numbracks; }
                if (d == '[') { --numbracks; }

                if (d == ']' && numbracks == 1) { nextIndex = ""; }
                else if (numbracks == 1 && comp[k-1..k+1]=="..") {
                    isSlice = true;
                    sliceTo = nextIndex;
                    nextIndex = "";
                    --k;
                } else if ((d == '[' && numbracks==0) || (d==',' && numbracks==1)) {
                    if (isSlice && numSlicesRemaining>0) {
                        if (numSlicesRemaining==1) {
                            if (!hasDollar(nextIndex) && !hasDollar(sliceTo)) {
                                return "(" ~ sliceTo ~ "-" ~ nextIndex ~ ")";
                            }
                        }
                        v = "[" ~ nextIndex ~ ".." ~ sliceTo
                         ~ "].length";
                        --numSlicesRemaining;
                    } else {
                        if (isSlice)
                            v = "[" ~ nextIndex ~ ".." ~ sliceTo ~ "]" ~ v;
                        else v = "[" ~ nextIndex ~ "]" ~ v;
                    }
                    nextIndex = "";
                    isSlice = false;
                } else {
                    if (d>='A' && d<='Z') nextIndex = tree.symbolTable[d-'A'].value ~ nextIndex;
                    else nextIndex = d ~ nextIndex;
                }
            }
            if (numSlicesRemaining>0) v ~= ".length";
            return tree.symbolTable[comp[0]-'A'].