Merge remote-tracking branch 'origin/master' into hydro

16 files changed, 1418 insertions, 509 deletions

diff --git a/Tools/CodeGen/CactusBoundary.m b/Tools/CodeGen/CactusBoundary.m
index 090151c..610987b 100644
--- a/Tools/CodeGen/CactusBoundary.m
+++ b/Tools/CodeGen/CactusBoundary.m
@@ -69,7 +69,7 @@ GetScheduledGroups[thornName_] :=
  {Name          -> "ApplyBCs",
   Language      -> "None", (* groups do not have a language *)
   SchedulePoint -> "as " <> thornName <> "_ApplyBCs in MoL_PostStep "
-                   <> " after " <> boundariesName[thornName], 
+                   <> "after " <> boundariesName[thornName], 
   Comment       -> "Apply boundary conditions "
                    <> "controlled by thorn Boundary"
  }
@@ -95,10 +95,10 @@ GetScheduledFunctions[thornName_, evolvedGroups_] :=
   }
 };
 
-createBoundTypeParam[groupOrGF_] := {
+createBoundTypeParam[groupOrGF_, def_] := {
                  Name          ->  ToString@groupOrGF <> "_bound",
                  Type          ->  "KEYWORD",
-                 Default       ->  "skip",
+                 Default       ->  def,
                  Description   ->  "Boundary condition to implement",
                  Visibility    ->  "private",
                  AllowedValues ->  {
@@ -108,8 +108,9 @@ createBoundTypeParam[groupOrGF_] := {
         {Value -> "radiative", Description -> "Radiation boundary condition"},
         {Value -> "scalar",    Description -> "Dirichlet boundary condition"},
         {Value -> "newrad",    Description -> "Improved radiative boundary condition"},
-        {Value -> "skip",      Description -> "skip boundary condition code"}
-}};
+        {Value -> "skip",      Description -> "skip boundary condition code"}},
+                 Steerable -> Always
+};
 
 
 createBoundSpeedParam[groupOrGF_] := {
@@ -119,7 +120,8 @@ createBoundSpeedParam[groupOrGF_] := {
                  Description   ->  "characteristic speed at boundary",
                  Visibility    ->  "private",
                  AllowedValues ->  {{Value -> "0:*" ,
-                      Description -> "outgoing characteristic speed > 0"}}
+                      Description -> "outgoing characteristic speed > 0"}},
+                 Steerable -> Always
 };
 
 createBoundLimitParam[groupOrGF_] := {
@@ -129,7 +131,8 @@ createBoundLimitParam[groupOrGF_] := {
                  Description   ->  "limit value for r -> infinity",
                  Visibility    ->  "private",
                  AllowedValues ->  {{Value -> "*:*" ,
-                      Description -> "value of limit value is unrestricted"}}
+                      Description -> "value of limit value is unrestricted"}},
+                 Steerable -> Always
 };
 
 createBoundScalarParam[groupOrGF_] := {
@@ -139,12 +142,13 @@ createBoundScalarParam[groupOrGF_] := {
                  Description   ->  "Dirichlet boundary value",
                  Visibility    ->  "private",
                  AllowedValues ->  {{Value -> "*:*" ,
-                      Description -> "unrestricted"}}
+                      Description -> "unrestricted"}},
+                 Steerable -> Always
 };
 
 GetParameters[evolvedGFs_, evolvedGroups_] :=
-  Join[Map[createBoundTypeParam, evolvedGFs],
-       Map[createBoundTypeParam, Map[unqualifiedGroupName,evolvedGroups]],
+  Join[Map[createBoundTypeParam[#,"skip"] &, evolvedGFs],
+       Map[createBoundTypeParam[#,"none"] &, Map[unqualifiedGroupName,evolvedGroups]],
 
        Map[createBoundSpeedParam, evolvedGFs],
        Map[createBoundSpeedParam, Map[unqualifiedGroupName,evolvedGroups]], 
@@ -169,7 +173,7 @@ GetSources[evolvedGroups_, groups_, implementation_, thornName_] :=
       ExcisionGFs -> evolvedGFs
     };
 
-    Return[{{Filename -> "Boundaries.c", 
+    Return[{{Filename -> "Boundaries.cc", 
              Contents -> CreateMoLBoundariesSource[boundarySpec]}}]];
 
 End[];
diff --git a/Tools/CodeGen/CalculationFunction.m b/Tools/CodeGen/CalculationFunction.m
index 58e7a5f..4f236ef 100644
--- a/Tools/CodeGen/CalculationFunction.m
+++ b/Tools/CodeGen/CalculationFunction.m
@@ -20,7 +20,7 @@
 
 BeginPackage["CalculationFunction`", {"CodeGen`",
   "MapLookup`", "KrancGroups`", "Differencing`", "Errors`",
-  "Helpers`", "Kranc`"}];
+  "Helpers`", "Kranc`", "Optimize`", "Jacobian`"}];
 
 CreateCalculationFunction::usage = "";
 VerifyCalculation::usage = "";
@@ -159,6 +159,26 @@ removeUnusedShorthands[calc_] :=
       removeUnusedShorthands[newCalc],
       newCalc]];
 
+(* Return all the groups that are used in a given calculation *)
+groupsInCalculation[calc_, imp_] :=
+  Module[{groups,gfs,eqs,gfsUsed, groupNames},
+    groups = lookup[calc, Groups];
+    gfs = allGroupVariables[groups];
+    eqs = lookup[calc, Equations];
+    gfsUsed = Union[Cases[eqs, _ ? (MemberQ[gfs,#] &), Infinity]];
+    groupNames = containingGroups[gfsUsed, groups];
+    Map[qualifyGroupName[#, imp] &, groupNames]];
+
+CheckGroupStorage[groupNames_, calcName_] :=
+  Module[{ignoreGroups, groupsNames2},
+    ignoreGroups = {"TmunuBase::stress_energy_scalar", "TmunuBase::stress_energy_vector",
+      "TmunuBase::stress_energy_tensor"};
+    groupNames2 = Select[groupNames, !MemberQ[ignoreGroups, #] &];
+    {"\nconst char *groups[] = {",
+    Riffle[Map[Quote,groupNames2], ","],
+    "};\n",
+    "GenericFD_AssertGroupStorage(cctkGH, ", Quote[calcName],", ", Length[groupNames2], ", groups);\n"}];
+
 (* --------------------------------------------------------------------------
    Variables
    -------------------------------------------------------------------------- *)
@@ -180,7 +200,7 @@ localName[x_] :=
 
 definePreDefinitions[pDefs_] :=
   CommentedBlock["Initialize predefined quantities",
-    Map[DeclareAssignVariable["CCTK_REAL", #[[1]], #[[2]]] &, pDefs]];
+    Map[DeclareAssignVariable[DataType[], #[[1]], #[[2]]] &, pDefs]];
 
 (* --------------------------------------------------------------------------
    Equations
@@ -203,7 +223,7 @@ printEq[eq_] :=
     split[ x_ + y___] := { x, " + ..."};
     split[-x_ + y___] := {-x, " + ..."};
     split[ x_       ] := { x, ""};
-    rhsSplit = split[Expand@ReplacePowers@rhs];
+    rhsSplit = split[Expand@ReplacePowers[rhs,False]];
     rhsString = ToString@CForm[rhsSplit[[1]]] <> rhsSplit[[2]];
     InfoMessage[InfoFull, " " <> ToString@lhs <> " -> " <> rhsString]];
 
@@ -228,11 +248,11 @@ simpCollect[collectList_, eqrhs_, localvar_, debug_] :=
     all];
 
 (* Return a CodeGen block which assigns dest by evaluating expr *)
-assignVariableFromExpression[dest_, expr_, declare_] :=
+assignVariableFromExpression[dest_, expr_, declare_, vectorise_] :=
   Module[{tSym, type, cleanExpr, code},
     tSym = Unique[];
-    type = If[StringMatchQ[ToString[dest], "dir*"], "int", "CCTK_REAL"];
-    cleanExpr = ReplacePowers[expr] /. Kranc`t -> tSym;
+    type = If[StringMatchQ[ToString[dest], "dir*"], "ptrdiff_t", DataType[]];
+    cleanExpr = ReplacePowers[expr, vectorise] /. Kranc`t -> tSym;
 
     If[SOURCELANGUAGE == "C",
       code = If[declare, type <> " ", ""] <> ToString[dest] <> " = " <>
@@ -258,6 +278,34 @@ assignVariableFromExpression[dest_, expr_, declare_] :=
 
     {code}];
 
+(* This and assignVariableFromExpression should be combined *)
+generateCodeFromExpression[expr_, vectorise_] :=
+  Module[{tSym, type, cleanExpr, code},
+    tSym = Unique[];
+    cleanExpr = ReplacePowers[expr, vectorise] /. Kranc`t -> tSym;
+
+    If[SOURCELANGUAGE == "C",
+      code =
+        ToString[cleanExpr, CForm,         PageWidth -> Infinity],
+      code = ToString[cleanExpr, FortranForm, PageWidth -> 120]];
+
+    If[SOURCELANGUAGE != "C",
+      code = StringReplace[code, "\n  "      -> " &\n"];
+      code = StringReplace[code, "   -  "    -> " &  "];
+      code = StringReplace[code, ".eq."      -> " = "];
+      code = StringReplace[code, "=        " -> "="];
+      code = StringReplace[code, "\\"        -> ""];
+      code = StringReplace[code, "(index)"   -> "(i,j,k)"]];
+
+    code = StringReplace[code, "normal1"     -> "normal[0]"];
+    code = StringReplace[code, "normal2"     -> "normal[1]"];
+    code = StringReplace[code, "normal3"     -> "normal[2]"];
+    code = StringReplace[code, "BesselJ"-> "gsl_sf_bessel_Jn"];
+    code = StringReplace[code, ToString@tSym -> "cctk_time"];
+    code = StringReplace[code, "\"" -> ""];
+
+    {code}];
+
 (* --------------------------------------------------------------------------
    Shorthands
    -------------------------------------------------------------------------- *)
@@ -326,24 +374,32 @@ pdCanonicalOrdering[name_[inds___] -> x_] :=
 
 Options[CreateCalculationFunction] = ThornOptions;
 
-CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] :=
+CreateCalculationFunction[calcp_, debug_, imp_, opts:OptionsPattern[]] :=
   Module[{gfs, allSymbols, knownSymbols,
-          shorts, eqs, parameters,
+          shorts, eqs, parameters, parameterRules,
           functionName, dsUsed, groups, pddefs, cleancalc, eqLoop, where,
-          addToStencilWidth, pDefs, haveCondTextuals, condTextuals},
+          addToStencilWidth, pDefs, haveCondTextuals, condTextuals, calc},
+
+  calc = If[OptionValue[UseJacobian], InsertJacobian[calcp, opts], calcp];
 
   cleancalc = removeUnusedShorthands[calc];
+  If[OptionValue[CSE],
+    cleancalc = EliminateCommonSubexpressions[cleancalc, opts]];
+
   shorts = lookupDefault[cleancalc, Shorthands, {}];
 
   eqs    = lookup[cleancalc, Equations];
   parameters = lookupDefault[cleancalc, Parameters, {}];
   groups = lookup[cleancalc, Groups];
+  If[OptionValue[UseJacobian], groups = Join[groups, JacobianGroups[]]];
   pddefs = lookupDefault[cleancalc, PartialDerivatives, {}];
   where = lookupDefault[cleancalc, Where, Everywhere];
   addToStencilWidth = lookupDefault[cleancalc, AddToStencilWidth, 0];
   pDefs = lookup[cleancalc, PreDefinitions];
   haveCondTextuals = mapContains[cleancalc, ConditionalOnTextuals];
 
+  SetDataType[If[OptionValue[UseVectors],"CCTK_REAL_VEC", "CCTK_REAL"]];
+
   VerifyCalculation[cleancalc];
 
   gfs = allGroupVariables[groups];
@@ -368,10 +424,14 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] :=
 
     If[!lookupDefault[cleancalc, NoSimplify, False],
        InfoMessage[InfoFull, "Simplifying equations", eqs];
-       eqs = Simplify[eqs, {r>0}]]];
+       eqs = Simplify[eqs, {r>=0}]]];
 
   InfoMessage[InfoFull, "Equations:"];
 
+  (* Wrap parameters with ToReal *)
+  parameterRules = Map[(#->ToReal[#])&, parameters];
+  eqs = eqs /. parameterRules;
+
   Map[printEq, eqs];
 
   (* Check all the function names *)
@@ -393,7 +453,8 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] :=
   allSymbols = calculationSymbols[cleancalc];
   knownSymbols = Join[lookupDefault[cleancalc, AllowedSymbols, {}], gfs, shorts, parameters,
     {dx,dy,dz,idx,idy,idz,t, Pi, E, Symbol["i"], Symbol["j"], Symbol["k"], normal1, normal2,
-    normal3, tangentA1, tangentA2, tangentA3, tangentB1, tangentB2, tangentB3}];
+    normal3, tangentA1, tangentA2, tangentA3, tangentB1, tangentB2, tangentB3},
+    If[OptionValue[UseJacobian], JacobianSymbols[], {}]];
 
   unknownSymbols = Complement[allSymbols, knownSymbols];
 
@@ -402,7 +463,7 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] :=
        "Calculation is:", cleancalc]];
 
   {
-  DefineFunction[bodyFunctionName, "void", "cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]",
+  DefineFunction[bodyFunctionName, "static void", "cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]",
   {
     "DECLARE_CCTK_ARGUMENTS;\n",
     "DECLARE_CCTK_PARAMETERS;\n\n",
@@ -415,14 +476,21 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] :=
     ConditionalOnParameterTextual["cctk_iteration % " <> functionName <> "_calc_every != " <>
       functionName <> "_calc_offset", "return;\n"],
 
+    CheckGroupStorage[groupsInCalculation[cleancalc, imp], functionName],
+
+    "\n",
+    CheckStencil[pddefs, eqs, functionName],
+
     If[haveCondTextuals, Map[ConditionalOnParameterTextual["!(" <> # <> ")", "return;\n"] &,condTextuals], {}],
 
     CommentedBlock["Include user-supplied include files",
       Map[IncludeFile, lookupDefault[cleancalc, DeclarationIncludes, {}]]],
 
-    InitialiseFDVariables[],
+    InitialiseFDVariables[OptionValue[UseVectors]],
     definePreDefinitions[pDefs],
 
+    If[OptionValue[UseJacobian], CreateJacobianVariables[], {}],
+
     If[Cases[{pddefs}, SBPDerivative[_], Infinity] != {},
       CommentedBlock["Compute Summation By Parts derivatives",
         IncludeFile["sbp_calc_coeffs.h"]],
@@ -431,7 +499,7 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] :=
     If[gfs != {},
       {
 	eqLoop = equationLoop[eqs, cleancalc, gfs, shorts, {}, groups,
-          pddefs, where, addToStencilWidth, useCSE, opts]},
+          pddefs, where, addToStencilWidth, opts]},
 
        ConditionalOnParameterTextual["verbose > 1",
          "CCTK_VInfo(CCTK_THORNSTRING,\"Leaving " <> bodyFunctionName <> "\");\n"],
@@ -464,12 +532,13 @@ CreateCalculationFunction[calc_, debug_, useCSE_, opts:OptionsPattern[]] :=
 Options[equationLoop] = ThornOptions;
 
 equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
-             where_, addToStencilWidth_, useCSE_,
+             where_, addToStencilWidth_,
              opts:OptionsPattern[]] :=
   Module[{rhss, lhss, gfsInRHS, gfsInLHS, gfsOnlyInRHS, localGFs,
           localMap, eqs2, derivSwitch, code, functionName, calcCode,
           loopFunction, gfsInBoth, gfsDifferentiated,
-          gfsDifferentiatedAndOnLHS, declare, eqsReplaced},
+          gfsDifferentiatedAndOnLHS, declare, eqsReplaced,
+          generateEquationsCode},
 
     rhss = Map[#[[2]] &, eqs];
     lhss = Map[#[[1]] &, eqs];
@@ -520,9 +589,6 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
     (* Replace grid functions with their local forms *)
     eqsReplaced = eqs2 /. localMap;
 
-    If[useCSE,
-      eqsReplaced = CSE[eqsReplaced]];
-
     (* Construct a list, corresponding to the list of equations,
        marking those which need their LHS variables declared.  We
        declare variables at the same time as assigning to them as it
@@ -531,22 +597,80 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
        functions which appear in the RHSs have been declared and set
        already (DeclareMaybeAssignVariableInLoop below), so assignments
        to these do not generate declarations here. *)
-    declare = markFirst[First /@ eqsReplaced, Map[localName, gfsInRHS]];
+    declare = Block[{$RecursionLimit=Infinity},markFirst[First /@ eqsReplaced, Map[localName, gfsInRHS]]];
+
+    (* Generate the code for the equations, factoring out any
+       sequential IfThen statements with the same condition.  Try to
+       declare variables as they are assigned, but it is only possible
+       to do this outside all if(){} statements. "declare2" is a list
+       of booleans indicating if the LHS variables should be declared
+       as they are assigned.  *)
+    generateEquationsCode[eqs2_, declare2_] :=
+      Module[{ifs, ind, cond, num, preDeclare},
+        ifs = Position[eqs2, _ -> IfThen[__], {1}];
+        If[Length[ifs] <= 1,
+          Riffle[MapThread[assignVariableFromExpression[#1[[1]], #1[[2]], #2, OptionValue[UseVectors]] &,
+            {eqs2, declare2}],"\n"],
+        (* else *)
+          ind = ifs[[1,1]];
+          If[ind > 1,
+            {generateEquationsCode[Take[eqs2, ind-1], Take[declare2, ind-1]], "\n",
+             generateEquationsCode[Drop[eqs2, ind-1], Drop[declare2, ind-1]]},
+          (* else *)
+            cond = eqs2[[1,2,1]];
+            num = LengthWhile[eqs2, MatchQ[#, _ -> IfThen[cond, _, _]] &];
+            If[num == 1,
+              {generateEquationsCode[Take[eqs2,1], Take[declare2,1]], "\n",
+               generateEquationsCode[Drop[eqs2,1], Drop[declare2,1]]},
+            (* else *)
+              e1 = Take[eqs2, num];
+              e2 = Drop[eqs2, num];
+              preDeclare = #[[2,1]] & /@ Select[MapThread[List, {Take[declare2,num], e1}], #[[1]] &];
+              {Map[DeclareVariableNoInit[#, DataType[]] &, Complement[Union[preDeclare], localName/@gfsInRHS]], {"\n"},
+              {Conditional[generateCodeFromExpression[cond, OptionValue[UseVectors]],
+                generateEquationsCode[Map[#[[1]] -> #[[2,2]]&, e1], ConstantArray[False, Length[e1]]],
+                generateEquationsCode[Map[#[[1]] -> #[[2,3]]&, e1], ConstantArray[False, Length[e1]]]],
+              "\n",
+              generateEquationsCode[e2,Drop[declare2,Length[e1]]]}}]]]];
+
+    calcCode = generateEquationsCode[eqsReplaced, declare];
+      
+
+    assignLocalGridFunctions[gs_, useVectors_, useJacobian_] :=
+      Module[{conds, varPatterns, varsInConds, simpleVars, code},
+        conds =
+          {{"eT" ~~ _ ~~ _, "*stress_energy_state", "ToReal(0.0)"}}; (* This should be passed as an option *)
+        If[useJacobian,
+          conds = Append[conds, JacobianConditionalGridFunctions[]]];
+
+        varPatterns = Map[First, conds];
+        varsInConds = Map[Function[pattern, Select[gs,StringMatchQ[ToString[#], pattern] &]], varPatterns];
+        simpleVars = Complement[gs, Flatten[varsInConds]];
+        code = {"\n",
+         Map[DeclareMaybeAssignVariableInLoop[
+              DataType[], localName[#], GridName[#],
+              False,"", useVectors] &, simpleVars],
+         {"\n",
+         Riffle[
+           MapThread[
+             If[Length[#2] > 0,
+               {DeclareVariables[localName/@#2, DataType[]],"\n",
+                Conditional[#1,
+                  Table[AssignVariableInLoop[localName[var], GridName[var], useVectors], {var, #2}],
+                  Sequence@@If[#3 =!= None, {Table[AssignVariableInLoop[localName[var], #3, False (*useVectors*)], {var, #2}]}, {}]]},
+               (* else *)
+               {}] &,
+             {Map[#[[2]]&, conds], varsInConds, Map[#[[3]]&, conds]}], "\n"]}};
+         code
+      ];
 
-    calcCode =
-      MapThread[{assignVariableFromExpression[#1[[1]], #1[[2]], #2], "\n"} &,
-       {eqsReplaced, declare}];
 
     GenericGridLoop[functionName,
     {
-      DeclareDerivatives[defsWithoutShorts, eqsOrdered],
+      (* DeclareDerivatives[defsWithoutShorts, eqsOrdered], *)
 
       CommentedBlock["Assign local copies of grid functions",
-        Map[DeclareMaybeAssignVariableInLoop[
-              "CCTK_REAL", localName[#], GridName[#],
-              StringMatchQ[ToString[GridName[#]], "eT" ~~ _ ~~ _ ~~ "[" ~~ __ ~~ "]"],
-                "*stress_energy_state"] &,
-            gfsInRHS]],
+        assignLocalGridFunctions[gfsInRHS, OptionValue[UseVectors], OptionValue[UseJacobian]]],
 
       CommentedBlock["Include user supplied include files",
         Map[IncludeFile, incs]],
@@ -560,8 +684,34 @@ equationLoop[eqs_, cleancalc_, gfs_, shorts_, incs_, groups_, pddefs_,
         Map[InfoVariable[#[[1]]] &, (eqs2 /. localMap)],
         ""],
 
+      If[OptionValue[UseVectors], {
+        CommentedBlock["If necessary, store only partial vectors after the first iteration",
+          ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 2 && CCTK_BUILTIN_EXPECT(i < lc_imin && i+CCTK_REAL_VEC_SIZE > lc_imax, 0)",
+            {
+              DeclareAssignVariable["ptrdiff_t", "elt_count_lo", "lc_imin-i"],
+              DeclareAssignVariable["ptrdiff_t", "elt_count_hi", "lc_imax-i"],
+              Map[StoreMiddlePartialVariableInLoop[GridName[#], localName[#], "elt_count_lo", "elt_count_hi"] &,
+                  gfsInLHS],
+              "break;\n"
+            }]],
+        CommentedBlock["If necessary, store only partial vectors after the first iteration",
+          ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 1 && CCTK_BUILTIN_EXPECT(i < lc_imin, 0)",
+            {
+              DeclareAssignVariable["ptrdiff_t", "elt_count", "lc_imin-i"],
+              Map[StoreHighPartialVariableInLoop[GridName[#], localName[#], "elt_count"] &,
+                  gfsInLHS],
+              "continue;\n"
+            }]],
+        CommentedBlock["If necessary, store only partial vectors after the last iteration",
+          ConditionalOnParameterTextual["CCTK_REAL_VEC_SIZE > 1 && CCTK_BUILTIN_EXPECT(i+CCTK_REAL_VEC_SIZE > lc_imax, 0)",
+            {
+              DeclareAssignVariable["ptrdiff_t", "elt_count", "lc_imax-i"],
+              Map[StoreLowPartialVariableInLoop[GridName[#], localName[#], "elt_count"] &,
+                  gfsInLHS],
+              "break;\n"
+            }]]}, {}],
       CommentedBlock["Copy local copies back to grid functions",
-        Map[AssignVariableInLoop[GridName[#], localName[#]] &,
+        Map[(If[OptionValue[UseVectors], StoreVariableInLoop, AssignVariableInLoop][GridName[#], localName[#]]) &,
             gfsInLHS]],
 
       If[debugInLoop, Map[InfoVariable[GridName[#]] &, gfsInLHS], ""]}, opts]];
diff --git a/Tools/CodeGen/CodeGen.m b/Tools/CodeGen/CodeGen.m
index 7cae5fb..5e8cd6a 100644
--- a/Tools/CodeGen/CodeGen.m
+++ b/Tools/CodeGen/CodeGen.m
@@ -43,6 +43,9 @@ IncludeSystemFile::usage = "IncludeFile[name] returns a block of code" <>
 DeclareVariable::usage = "DeclareVariable[name, type] returns a block of code " <>
   "that declares a variable of given name and type.  'name' and 'type' should be " <>
   "strings.";
+DeclareVariableNoInit::usage = "DeclareVariableNoInit[name, type] returns a block of code " <>
+  "that declares a variable of given name and type without initialising it.  'name' and 'type' should be " <>
+  "strings.";
 DeclareVariables::usage = "DeclareVariables[names, type] returns a block of code " <>
   "that declares a list of variables of given name and type.  'names' should be a list" <>
   " of strings and 'type' should be a string string.";
@@ -60,14 +63,20 @@ DeclareAssignVariable::usage = "DeclareAssignVariable[type_, dest_, src_] return
   "that declares and sets a constant variable of given name and type.";
 AssignVariableInLoop::usage = "AssignVariableInLoop[dest_, src_] returns a block of code " <>
   "that assigns 'src' to 'dest'.";
+StoreVariableInLoop::usage = "StoreVariableInLoop[dest_, src_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+StoreLowPartialVariableInLoop::usage = "StoreLowPartialVariableInLoop[dest_, src_, count_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+StoreHighPartialVariableInLoop::usage = "StoreHighPartialVariableInLoop[dest_, src_, count_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
+StoreMiddlePartialVariableInLoop::usage = "StoreMiddlePartialVariableInLoop[dest_, src_, countLow_, countHigh_] returns a block of code " <>
+  "that assigns 'src' to 'dest'.";
 DeclareAssignVariableInLoop::usage = "DeclareAssignVariableInLoop[type_, dest_, src_] returns a block of code " <>
   "that assigns 'src' to 'dest'.";
 MaybeAssignVariableInLoop::usage = "MaybeAssignVariableInLoop[dest_, src_, cond_] returns a block of code " <>
   "that assigns 'src' to 'dest'.";
 DeclareMaybeAssignVariableInLoop::usage = "DeclareMaybeAssignVariableInLoop[type_, dest_, src_, cond_] returns a block of code " <>
   "that assigns 'src' to 'dest'.";
-DeclareVariablesInLoopVectorised::usage = "";
-AssignVariablesInLoopVectorised::usage = "";
 TestForNaN::usage = "TestForNaN[expr_] returns a block of code " <>
   "that tests 'expr' for nan.";
 CommentedBlock::usage = "CommentedBlock[comment, block] returns a block consisting " <>
@@ -117,7 +126,6 @@ CommaNewlineSeparated::usage = ""; (* This should not really be in CodeGen *)
 CommaSeparated::usage = "";
 ReplacePowers::usage = "";
 CFormHideStrings::usage = "";
-CSE::usage = "";
 BoundaryLoop::usage = "";
 BoundaryWithGhostsLoop::usage = "";
 GenericGridLoop::usage = "";
@@ -125,15 +133,18 @@ GenericGridLoop::usage = "";
 NameRoot::usage = "";
 PartitionVarList::usage = "";
 Quote::usage = "Quote[x] returns x surrounded by quotes";
+DataType::usage = "DataType[] returns a string for the grid function data type (e.g. CCTK_REAL)";
+SetDataType::usage = "SetDataType[type] sets a string for the grid function data type (e.g. CCTK_REAL)";
+Conditional;
 
 Begin["`Private`"];
 
 SOURCELANGUAGE =  "C";
-SOURCESUFFIX   = ".c";
+SOURCESUFFIX   = ".cc";
 
 setSourceSuffix[lang_] :=
 If[ (lang == "C"),
-  SOURCESUFFIX = ".c";
+  SOURCESUFFIX = ".cc";
 ,
   SOURCESUFFIX = ".F90";
 ];
@@ -146,10 +157,18 @@ If[ (lang == "C" || lang == "Fortran"),
   InfoMessage[Terse, "User set source language to " <> lang],
 
   SOURCELANGUAGE = "C";
-  setSourceSuffix[".c"];
+  setSourceSuffix[".cc"];
   InfoMessage[Terse, "Setting Source Language to C"];
 ];
 
+SetDataType[type_String] :=
+  dataType = type;
+
+DataType[] :=
+  If[dataType === Symbol["datatype"],
+    Throw["DataType: Have not set a data type"],
+    dataType];
+
 (* Code generation utilities; not specific to any language *)
 
 FlattenBlock[b_] := Apply[StringJoin,Map[ToString,If[! AtomQ[b], Flatten[b, Infinity], b]]];
@@ -234,6 +253,12 @@ If[SOURCELANGUAGE == "C",
   {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *)
   ];
 
+DeclareVariableNoInit[name_, type_] :=
+If[SOURCELANGUAGE == "C",
+  {type, " ",    name, EOL[]},
+  {type, " :: ", name, EOL[]} (* no value init here to avoid implicit SAVE attribute *)
+  ];
+
 
 DeclareVariables[names_?ListQ, type_] := 
 If[SOURCELANGUAGE == "C",
@@ -273,69 +298,41 @@ AssignVariable[dest_, src_] :=
 DeclareAssignVariable[type_, dest_, src_] :=
   {type, " const ", dest, " = ", src, EOL[]};
 
-AssignVariableInLoop[dest_, src_] :=
-  {dest, " = ", src, EOL[]};
+AssignVariableInLoop[dest_, src_, vectorise_:False] :=
+  Module[{loader},
+    loader[x_] := If[vectorise, {"vec_load(", x, ")"}, x];
+    {dest, " = ", loader[src], EOL[]}];
 (*
   {dest, " = ", src, EOL[],
    TestForNaN[dest]};
 *)
 
-DeclareAssignVariableInLoop[type_, dest_, src_] :=
-  {type, " const ", dest, " = ", src, EOL[]};
+StoreVariableInLoop[dest_, src_] :=
+  {"vec_store_nta(", dest, ",", src, ")", EOL[]};
 
-MaybeAssignVariableInLoop[dest_, src_, cond_] :=
-  If [cond,
-      {dest, " = useMatter ? ", src, " : 0.0", EOL[]},
-      {dest, " = ", src, EOL[]}];
+StoreLowPartialVariableInLoop[dest_, src_, count_] :=
+  {"vec_store_nta_partial_lo(", dest, ",", src, ",", count, ")", EOL[]};
 
-DeclareMaybeAssignVariableInLoop[type_, dest_, src_, mmaCond_, codeCond_] :=
-  If [mmaCond,
-      {type, "  ", dest, " = (", codeCond, ") ? (", src, ") : 0.0", EOL[]},
-      {type, "  ", dest, " = ", src, EOL[]}];
+StoreHighPartialVariableInLoop[dest_, src_, count_] :=
+  {"vec_store_nta_partial_hi(", dest, ",", src, ",", count, ")", EOL[]};
 
-(* TODO: move these into OpenMP loop *)
-DeclareVariablesInLoopVectorised[dests_, temps_, srcs_] :=
-  {
-   {"#undef LC_PRELOOP_STATEMENTS", "\n"},
-   {"#define LC_PRELOOP_STATEMENTS", " \\\n"},
-   {"int const GFD_imin = lc_imin + ((lc_imin + cctk_lsh[0] * (j + cctk_lsh[1] * k)) & (CCTK_REAL_VEC_SIZE-1))", "; \\\n"},
-   {"int const GFD_imax = lc_imax + ((lc_imax + cctk_lsh[0] * (j + cctk_lsh[1] * k)) & (CCTK_REAL_VEC_SIZE-1)) - CCTK_REAL_VEC_SIZE", "; \\\n"},
-   Map[Function[x, Module[{dest, temp, src},
-                          {dest, temp, src} = x;
-                          {"CCTK_REAL_VEC ", temp, "; \\\n"}]],
-       Transpose[{dests, temps, srcs}]],
-   {"\n"}
-  };
+StoreMiddlePartialVariableInLoop[dest_, src_, countLow_, countHigh_] :=
+  {"vec_store_nta_partial_mid(", dest, ",", src, ",", countLow, ",", countHigh, ")", EOL[]};
 
-AssignVariablesInLoopVectorised[dests_, temps_, srcs_] :=
-  {
-   {"{\n"},
-   {"  if (i < GFD_imin || i >= GFD_imax) {\n"},
-   Map[Function[x, Module[{dest, temp, src},
-                          {dest, temp, src} = x;
-                          {"    ", dest, "[index] = ", src, EOL[]}]],
-       Transpose[{dests, temps, srcs}]],
-   {"  } else {\n"},
-   {"    size_t const index0 = index & (CCTK_REAL_VEC_SIZE-1)", EOL[]},
-   Map[Function[x, Module[{dest, temp, src},
-                          {dest, temp, src} = x;
-                          {"    ((CCTK_REAL*)&", temp, ")[index0] = ",
-                           src, EOL[]}]],
-       Transpose[{dests, temps, srcs}]],
-   {"    if (index0 == CCTK_REAL_VEC_SIZE-1) {\n"},
-   {"      size_t const index1 = index - (CCTK_REAL_VEC_SIZE-1)", EOL[]},
-   Map[Function[x, Module[{dest, temp, src},
-                          {dest, temp, src} = x;
-                          {"      _mm_stream_pd (&", dest, "[index1], ",
-                           temp, ")", EOL[]}]],
-       Transpose[{dests, temps, srcs}]],
-   {"    }\n"},
-   {"  }\n"},
-   {"}\n"}
-  };
+DeclareAssignVariableInLoop[type_, dest_, src_] :=
+  {type, " const ", dest, " = vec_load(", src, ")", EOL[]};
+
+MaybeAssignVariableInLoop[dest_, src_, cond_] :=
+  If [cond,
+      {dest, " = useMatter ? vec_load(", src, ") : ToReal(0.0)", EOL[]},
+      {dest, " = vec_load(", src, ")", EOL[]}];
 
-AssignVariableInLoopsVectorised[dest_, temp_, src_] :=
-  {"GFD_save_and_store(", dest, ",", "index", ",", "&", temp, ",", src, ")", EOL[]};
+DeclareMaybeAssignVariableInLoop[type_, dest_, src_, mmaCond_, codeCond_, vectorise_:False] :=
+  Module[{loader},
+    loader[x_] := If[vectorise, {"vec_load(", x, ")"}, x];
+    If [mmaCond,
+        {type, " ", dest, " = (", codeCond, ") ? ", loader[src], " : ToReal(0.0)", EOL[]},
+        {type, " ", dest, " = ", loader[src], EOL[]}]];
 
 TestForNaN[expr_] :=
   {"if (isnan(", expr, ")) {\n",
@@ -343,7 +340,7 @@ TestForNaN[expr_] :=
    "  CCTK_VInfo(CCTK_THORNSTRING, \"ipos: %d %d %d\", i, j, k);\n",
    "  CCTK_VInfo(CCTK_THORNSTRING, \"lbnd: %d %d %d\", cctk_lbnd[0], cctk_lbnd[1], cctk_lbnd[2]);\n",
    "  CCTK_VInfo(CCTK_THORNSTRING, \"lsh: %d %d %d\", cctk_lsh[0], cctk_lsh[1], cctk_lsh[2]);\n",
-   "  CCTK_VInfo(CCTK_THORNSTRING, \"lssh: %d %d %d\", cctk_lssh[CCTK_LSSH_IDX(0,0)], cctk_lssh[CCTK_LSSH_IDX(0,1)], cctk_lssh[CCTK_LSSH_IDX(0,2)]);\n",
+   "  CCTK_VInfo(CCTK_THORNSTRING, \"LSSH: %d %d %d\", CCTK_LSSH(0,0), CCTK_LSSH(0,1), CCTK_LSSH(0,2));\n",
    "  CCTK_VInfo(CCTK_THORNSTRING, \"", expr, ": %.17g\", (double)", expr, ");\n",
    "}\n"};
 
@@ -406,7 +403,7 @@ defineSubroutine[name_, args_, contents_] :=
 
 defineSubroutineC[name_, args_, contents_] :=
   SeparatedBlock[
-    {"void ", name, "(", args, ")", "\n",
+    {"extern \"C\" void ", name, "(", args, ")", "\n",
      CBlock[contents]}];
 
 defineSubroutineF[name_, args_, contents_] :=
@@ -426,7 +423,7 @@ defineSubroutineF[name_, args_, contents_] :=
 
 (* This is a Cactus-callable function *)
 DefineCCTKFunction[name_, type_, contents_] :=
-  DefineFunction[name, type, "CCTK_ARGUMENTS", 
+  DefineFunction[name, "extern \"C\" " <> type, "CCTK_ARGUMENTS", 
     {
       "DECLARE_CCTK_ARGUMENTS;\n",
       "DECLARE_CCTK_PARAMETERS;\n\n",
@@ -451,47 +448,65 @@ DeclareFDVariables[] :=
                                               {khalf,kthird,ktwothird,kfourthird,keightthird},
                                               {"hdxi", "hdyi", "hdzi"}}],
      "\n"},
-    {Map[DeclareVariables[#, "int"] &, {{"di", "dj", "dk"}}],
+    {Map[DeclareVariables[#, "ptrdiff_t"] &, {{"di", "dj", "dk"}}],
      "\n"}];
 *)
   CommentedBlock["Declare finite differencing variables", {}];
 
-InitialiseFDSpacingVariablesC[] := 
+InitialiseFDSpacingVariablesC[vectorise_:False] := 
   {
-    DeclareAssignVariable["CCTK_REAL", "dx", "CCTK_DELTA_SPACE(0)"],
-    DeclareAssignVariable["CCTK_REAL", "dy", "CCTK_DELTA_SPACE(1)"],
-    DeclareAssignVariable["CCTK_REAL", "dz", "CCTK_DELTA_SPACE(2)"],
-    (* DeclareAssignVariable["int", "di", "CCTK_GFINDEX3D(cctkGH,1,0,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], *)
-    DeclareAssignVariable["int", "di", "1"],
-    DeclareAssignVariable["int", "dj", "CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"],
-    DeclareAssignVariable["int", "dk", "CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0)"]
+    (* DeclareAssignVariable["ptrdiff_t", "di", "CCTK_GFINDEX3D(cctkGH,1,0,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"], *)
+    DeclareAssignVariable["ptrdiff_t", "di", "1"],
+    DeclareAssignVariable["ptrdiff_t", "dj", "CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0)"],
+    DeclareAssignVariable["ptrdiff_t", "dk", "CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0)"],
+    If[vectorise,
+    {
+      DeclareAssignVariable["ptrdiff_t", "cdi", "sizeof(CCTK_REAL) * di"],
+      DeclareAssignVariable["ptrdiff_t", "cdj", "sizeof(CCTK_REAL) * dj"],
+      DeclareAssignVariable["ptrdiff_t", "cdk", "sizeof(CCTK_REAL) * dk"]
+    },
+    ""],
+    DeclareAssignVariable[DataType[], "dx", "ToReal(CCTK_DELTA_SPACE(0))"],
+    DeclareAssignVariable[DataType[], "dy", "ToReal(CCTK_DELTA_SPACE(1))"],
+    DeclareAssignVariable[DataType[], "dz", "ToReal(CCTK_DELTA_SPACE(2))"],
+    DeclareAssignVariable[DataType[], "dt", "ToReal(CCTK_DELTA_TIME)"]
   };
 
 InitialiseFDSpacingVariablesFortran[] := 
   {
+    AssignVariable["dt", "CCTK_DELTA_TIME"],
     AssignVariable["dx", "CCTK_DELTA_SPACE(1)"],
     AssignVariable["dy", "CCTK_DELTA_SPACE(2)"],
     AssignVariable["dz", "CCTK_DELTA_SPACE(3)"]
   }
 
 
-InitialiseFDVariables[] :=
+InitialiseFDVariables[vectorise_] :=
   CommentedBlock["Initialise finite differencing variables",
   { If[SOURCELANGUAGE == "Fortran",
        InitialiseFDSpacingVariablesFortran[],
-       InitialiseFDSpacingVariablesC[]],
+       InitialiseFDSpacingVariablesC[vectorise]],
     
-    DeclareAssignVariable["CCTK_REAL", "dxi", "1.0 / dx"],
-    DeclareAssignVariable["CCTK_REAL", "dyi", "1.0 / dy"],
-    DeclareAssignVariable["CCTK_REAL", "dzi", "1.0 / dz"],
-    DeclareAssignVariable["CCTK_REAL", "khalf", "0.5"],
-    DeclareAssignVariable["CCTK_REAL", "kthird", "1/3.0"],
-    DeclareAssignVariable["CCTK_REAL", "ktwothird", "2.0/3.0"],
-    DeclareAssignVariable["CCTK_REAL", "kfourthird", "4.0/3.0"],
-    DeclareAssignVariable["CCTK_REAL", "keightthird", "8.0/3.0"],
-    DeclareAssignVariable["CCTK_REAL", "hdxi", "0.5 * dxi"],
-    DeclareAssignVariable["CCTK_REAL", "hdyi", "0.5 * dyi"],
-    DeclareAssignVariable["CCTK_REAL", "hdzi", "0.5 * dzi"]}];
+    DeclareAssignVariable[DataType[], "dxi", "INV(dx)"],
+    DeclareAssignVariable[DataType[], "dyi", "INV(dy)"],
+    DeclareAssignVariable[DataType[], "dzi", "INV(dz)"],
+    If[vectorise,
+     {DeclareAssignVariable[DataType[], "khalf", "ToReal(0.5)"],
+      DeclareAssignVariable[DataType[], "kthird", "ToReal(1.0/3.0)"],
+      DeclareAssignVariable[DataType[], "ktwothird", "ToReal(2.0/3.0)"],
+      DeclareAssignVariable[DataType[], "kfourthird", "ToReal(4.0/3.0)"],
+      DeclareAssignVariable[DataType[], "keightthird", "ToReal(8.0/3.0)"],
+      DeclareAssignVariable[DataType[], "hdxi", "kmul(ToReal(0.5), dxi)"],
+      DeclareAssignVariable[DataType[], "hdyi", "kmul(ToReal(0.5), dyi)"],
+      DeclareAssignVariable[DataType[], "hdzi", "kmul(ToReal(0.5), dzi)"]},
+     {DeclareAssignVariable[DataType[], "khalf", "0.5"],
+      DeclareAssignVariable[DataType[], "kthird", "1/3.0"],
+      DeclareAssignVariable[DataType[], "ktwothird", "2.0/3.0"],
+      DeclareAssignVariable[DataType[], "kfourthird", "4.0/3.0"],
+      DeclareAssignVariable[DataType[], "keightthird", "8.0/3.0"],
+      DeclareAssignVariable[DataType[], "hdxi", "0.5 * dxi"],
+      DeclareAssignVariable[DataType[], "hdyi", "0.5 * dyi"],
+      DeclareAssignVariable[DataType[], "hdzi", "0.5 * dzi"]}]}];
 
 GridName[x_] := If[SOURCELANGUAGE == "C",
                    ToExpression[ToString[x] <> "[index]"],
@@ -543,9 +558,9 @@ InitialiseGridLoopVariables[derivativesUsedSwitch_, addToStencilWidth_] :=
   AssignVariable["kstart", arrayIndex["index_offset_z"]],
 
   "\n",
-  AssignVariable["iend", {arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,0)"], " - index_offset_x"}],
-  AssignVariable["jend", {arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,1)"], " - index_offset_y"}],
-  AssignVariable["kend", {arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,2)"], " - index_offset_z"}]
+  AssignVariable["iend", {"CCTK_LSSH(0,0) - index_offset_x"}],
+  AssignVariable["jend", {"CCTK_LSSH(0,1) - index_offset_y"}],
+  AssignVariable["kend", {"CCTK_LSSH(0,2) - index_offset_z"}]
   },
 
   {
@@ -554,9 +569,9 @@ InitialiseGridLoopVariables[derivativesUsedSwitch_, addToStencilWidth_] :=
   AssignVariable["kstart", arrayIndex[0]],
 
   "\n",
-  AssignVariable["iend", arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,0)"]],
-  AssignVariable["jend", arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,1)"]],
-  AssignVariable["kend", arrayElement["cctk_lssh", "CCTK_LSSH_IDX(0,2)"]]
+  AssignVariable["iend", "CCTK_LSSH(0,0)"],
+  AssignVariable["jend", "CCTK_LSSH(0,1)"],
+  AssignVariable["kend", "CCTK_LSSH(0,2)"]
   }]
 ];
 
@@ -629,7 +644,7 @@ Options[GenericGridLoop] = ThornOptions;
 
 GenericGridLoop[functionName_, block_, opts:OptionsPattern[]] :=
   If[OptionValue[UseLoopControl],
-    GenericGridLoopUsingLoopControl[functionName, block],
+    GenericGridLoopUsingLoopControl[functionName, block, OptionValue[UseVectors]],
     GenericGridLoopTraditional[block]];
 
 GenericGridLoopTraditional[block_] :=
@@ -647,24 +662,26 @@ GenericGridLoopTraditional[block_] :=
        }
         ]]]];
 
-GenericGridLoopUsingLoopControl[functionName_, block_] :=
+GenericGridLoopUsingLoopControl[functionName_, block_, vectorise_] :=
   If[SOURCELANGUAGE == "C",  
     CommentedBlock["Loop over the grid points",
       {
         "#pragma omp parallel\n",
-        "LC_LOOP3 (", functionName, ",\n",
-        "          i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],\n",
-        "          cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])\n",
+        If[vectorise, "LC_LOOP3VEC", "LC_LOOP3"] <> " (", functionName, ",\n",
+        "  i,j,k, min[0],min[1],min[2], max[0],max[1],max[2],\n",
+        "  cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]", If[vectorise, {",\n",
+        "  CCTK_REAL_VEC_SIZE"},""] <> ")\n",
         "{\n",
         indentBlock[
           {
-             DeclareVariable["index", "// int"],
-             DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"],
+             (* DeclareVariable["index", "// int"], *)
+             (* DeclareAssignVariable["int", "index", "CCTK_GFINDEX3D(cctkGH,i,j,k)"], *)
+             DeclareAssignVariable["ptrdiff_t", "index", "di*i + dj*j + dk*k"],
              block
           }
         ],
         "}\n",
-        "LC_ENDLOOP3 (", functionName, ");\n"
+        If[vectorise, "LC_ENDLOOP3VEC", "LC_ENDLOOP3"] <> " (", functionName, ");\n"
       }
     ],
     ""
@@ -691,9 +708,9 @@ BoundaryLoop[block_] :=
   AssignVariable[arrayElement["bmin", 0], "is_physbnd[0*2+0] ? 0 : imin[0]"],
   AssignVariable[arrayElement["bmin", 1], "is_physbnd[1*2+0] ? 0 : imin[1]"],
   AssignVariable[arrayElement["bmin", 2], "is_physbnd[2*2+0] ? 0 : imin[2]"],
-  AssignVariable[arrayElement["bmax", 0], "is_physbnd[0*2+1] ? cctk_from[CCTK_LSSH_IDX(0,0)] : imax[0]"],
-  AssignVariable[arrayElement["bmax", 1], "is_physbnd[1*2+1] ? cctk_from[CCTK_LSSH_IDX(0,1)] : imax[1]"],
-  AssignVariable[arrayElement["bmax", 2], "is_physbnd[2*2+1] ? cctk_from[CCTK_LSSH_IDX(0,2)] : imax[2]"]}], 
+  AssignVariable[arrayElement["bmax", 0], "is_physbnd[0*2+1] ? CCTK_LSSH(0,0) : imax[0]"],
+  AssignVariable[arrayElement["bmax", 1], "is_physbnd[1*2+1] ? CCTK_LSSH(0,1) : imax[1]"],
+  AssignVariable[arrayElement["bmax", 2], "is_physbnd[2*2+1] ? CCTK_LSSH(0,2) : imax[2]"]}], 
 
   CommentedBlock["Loop over all faces",
    loopOverInteger["dir", "0", "3",
@@ -704,7 +721,7 @@ BoundaryLoop[block_] :=
         {0,  {AssignVariable[arrayElement["bmax", "dir"], {arrayElement["imin", "dir"], ""}], 
               AssignVariable[arrayElement["bmin", "dir"], {0, ""}]}},
         {1,  {AssignVariable[arrayElement["bmin", "dir"], {arrayElement["imax", "dir"], "" }],
-              AssignVariable[arrayElement["bmax", "dir"], {"cctk_lssh[CCTK_LSSH_IDX(0,dir)]", ""}]}}]],
+              AssignVariable[arrayElement["bmax", "dir"], {"CCTK_LSSH(0,dir)", ""}]}}]],
        conditional[arrayElement["is_physbnd", "dir * 2 + face"],
          loopOverInteger["k", arrayElement["bmin",2], arrayElement["bmax",2],
            loopOverInteger["j", arrayElement["bmin",1], arrayElement["bmax",1],
@@ -728,9 +745,9 @@ BoundaryWithGhostsLoop[block_] :=
   AssignVariable[arrayElement["bmin", 0], "0"],
   AssignVariable[arrayElement["bmin", 1], "0"],
   AssignVariable[arrayElement["bmin", 2], "0"],
-  AssignVariable[arrayElement["bmax", 0], "cctk_lssh[CCTK_LSSH_IDX(0,0)]"],
-  AssignVariable[arrayElement["bmax", 1], "cctk_lssh[CCTK_LSSH_IDX(0,1)]"],
-  AssignVariable[arrayElement["bmax", 2], "cctk_lssh[CCTK_LSSH_IDX(0,2)]"]}], 
+  AssignVariable[arrayElement["bmax", 0], "CCTK_LSSH(0,0)"],
+  AssignVariable[arrayElement["bmax", 1], "CCTK_LSSH(0,1)"],
+  AssignVariable[arrayElement["bmax", 2], "CCTK_LSSH(0,2)"]}], 
 
   CommentedBlock["Loop over all faces",
    loopOverInteger["dir", "0", "3",
@@ -741,7 +758,7 @@ BoundaryWithGhostsLoop[block_] :=
         {0,  {AssignVariable[arrayElement["bmax", "dir"], {arrayElement["imin", "dir"], ""}], 
               AssignVariable[arrayElement["bmin", "dir"], {0, ""}]}},
         {1,  {AssignVariable[arrayElement["bmin", "dir"], {arrayElement["imax", "dir"], "" }],
-              AssignVariable[arrayElement["bmax", "dir"], {"cctk_lssh[CCTK_LSSH_IDX(0,dir)]", ""}]}}]],
+              AssignVariable[arrayElement["bmax", "dir"], {"CCTK_LSSH(0,dir)]", ""}]}}]],
        conditional[arrayElement["is_physbnd", "dir * 2 + face"],
          loopOverInteger["k", arrayElement["bmin",2], arrayElement["bmax",2],
            loopOverInteger["j", arrayElement["bmin",1], arrayElement["bmax",1],
@@ -757,12 +774,28 @@ BoundaryWithGhostsLoop[block_] :=
       ]}
      ]]]};
 
-conditional[condition_, block_] :=
- {"if (", condition, ")\n",
+(* Remove call to ToReal from condtion. This should really instead be
+   done much earlier. Sometimes Conditional is called with a
+   conditional that is an expression, sometimes it is a list of
+   strings, so we need to handle both cases. *)
+(* One approach to remove calls to ToReal is to wrap the condition
+  into a function call, such as e.g. Cond[x], which is then handled by
+  the vectorisation code in the same way the IfThen[x,y,z] expression
+  is handled there. *)
+removeToRealFunc[condition_] := Replace[condition, {ToReal[x_] -> x}]
+removeToRealString[condition_] := If[StringQ[condition], StringReplace[condition, "ToReal(" ~~ x__ ~~ ")" :> x], condition]
+removeToReal[condition_] := Map[removeToRealString, removeToRealFunc[condition]]
+
+Conditional[condition_, block_] :=
+ {"if (", removeToReal[condition], ")\n",
   CBlock[block]};
 
+Conditional[condition_, block1_, block2_] :=
+ {"if (", removeToReal[condition], ")\n",
+  CBlock[block1], "else\n", CBlock[block2]};
+
 onceInGridLoop[block_] :=
-  conditional["i == 5 && j == 5 && k == 5",
+  Conditional["i == 5 && j == 5 && k == 5",
               block];
 
 InfoVariable[name_] :=
@@ -789,47 +822,166 @@ insertFile[name_] :=
     Close[istream];
     contents];
 
+vectoriseExpression[exprp_] :=
+  Module[{isNotMinusOneQ, isNotTimesMinusOneQ, fmaRules, isNotKneg, arithRules, undoRules, expr},
+    expr = exprp;
+
+    (* Constants *)
+    (* expr = expr /. xx_Integer/; xx!=-1 :> ToReal[xx]; *)
+    expr = expr /. xx_Integer -> ToReal[xx];
+    expr = expr /. xx_Real -> ToReal[xx];
+    removeToRealRules = {
+      - ToReal[xx_] -> ToReal[- xx],
+      ToReal[xx_] + ToReal[yy_] -> ToReal[xx + yy],
+      ToReal[xx_] * ToReal[yy_] -> ToReal[xx * yy],
+      ToReal[xx_] == ToReal[yy_] -> ToReal[xx == yy],
+      ToReal[xx_] != ToReal[yy_] -> ToReal[xx != yy],
+      pow[xx_, ToReal[power_]] -> pow[xx, power],
+      (* keep the conditional expression scalar *)
+      IfThen[ToReal[xx_], yy_, zz_] -> IfThen[xx, yy, zz]
+    };
+    expr = expr //. removeToRealRules;
+
+    (* Replace all operators and functions *)
+    (* kneg, kadd, ksub, kmul, kdiv *)
+    isNotKneg[n_] := ! MatchQ[n,kneg[_]];
+    arithRules = {
+      - xx_ -> kneg[xx],
+      xx_ * yy_ -> kmul[xx,yy],
+      xx_ / yy_ -> kdiv[xx,yy],
+      xx_ + yy_ -> kadd[xx,yy],
+      xx_ - yy_ -> ksub[xx,yy],
+      kmul[-1,xx_]            -> kneg[xx],
+      kmul[xx_,-1]            -> kneg[xx],
+      kmul[ToReal[-1],xx_]    -> kneg[xx],
+      kmul[xx_,ToReal[-1]]    -> kneg[xx],
+      (* kmul[xx_,INV[yy_]]      -> kdiv[xx,yy], *)
+      (* kmul[INV[xx_],yy_]      -> kdiv[yy,xx], *)
+      kdiv[xx_,kdiv[yy_,zz_]] -> kdiv[kmul[xx,zz],yy],
+      kdiv[kdiv[xx_,yy_],zz_] -> kdiv[xx,kmul[yy,zz]],
+      kmul[kneg[xx_],yy_]     -> kneg[kmul[xx,yy]],
+      kmul[xx_,kneg[yy_]]     -> kneg[kmul[xx,yy]],
+      kdiv[kneg[xx_],yy_]     -> kneg[kdiv[xx,yy]],
+      kdiv[xx_,kneg[yy_]]     -> kneg[kdiv[xx,yy]],
+      kadd[kneg[xx_],yy_]     -> ksub[yy,xx],
+      ksub[kneg[xx_],yy_]     -> kneg[kadd[xx,yy]],
+      kadd[xx_,kneg[yy_]]     -> ksub[xx,yy],
+      ksub[xx_,kneg[yy_]]     -> kadd[xx,yy],
+      kneg[ksub[xx_,yy_]]     -> ksub[yy,xx],
+      Abs[xx_]      -> kfabs[xx],
+      Log[xx_]      -> klog[xx],
+      fabs[xx_]     -> kfabs[xx],
+      fmax[xx_,yy_] -> kfmax[xx,yy],
+      fmin[xx_,yy_] -> kfmin[xx,yy],
+      sqrt[xx_]     -> ksqrt[xx],
+      exp[xx_]      -> kexp[xx],
+      log[xx_]      -> klog[xx],
+      pow[xx_,yy_]  -> kpow[xx,yy],
+      kneg[kneg[xx_]]   -> xx,
+      kfabs[kneg[xx_]]  -> kfabs[xx],
+      kfnabs[kneg[xx_]] -> kfnabs[xx],
+      kneg[kfabs[xx_]]  -> kfnabs[xx],
+      kneg[kfnabs[xx_]] -> kfabs[xx]
+    };
+    expr = expr //. arithRules;
+
+    (* Undo some transformations *)
+    undoRules = {
+      IfThen[kmul[xx_,yy_], aa_, bb_] -> IfThen[xx*yy, aa, bb],
+      IfThen[kmul[xx_,yy_] != zz_, aa_, bb_] -> IfThen[xx*yy!=zz, aa, bb],
+      ToReal[kneg[xx_]] -> ToReal[-xx],
+      ToReal[kadd[xx_,yy_]] -> ToReal[xx+yy],
+      ToReal[ksub[xx_,yy_]] -> ToReal[xx-yy],
+      ToReal[kmul[xx_,yy_]] -> ToReal[xx*yy],
+      ToReal[xx_*kadd[yy_,zz_]] -> ToReal[xx*(yy+zz)],
+      kpow[xx_, kneg[power_]] -> kpow[xx, -power]
+    };
+    expr = expr //. undoRules;
+    
+    (* FMA (fused multiply-add) instructions *)
+    (* kmadd (x,y,z) =   xy+z
+       kmsub (x,y,z) =   xy-z
+       knmadd(x,y,z) = -(xy+z)
+       knmsub(x,y,z) = -(xy-z) *)
+    fmaRules = {
+      kadd[kmul[xx_,yy_],zz_] -> kmadd[xx,yy,zz],
+      kadd[zz_,kmul[xx_,yy_]] -> kmadd[xx,yy,zz],
+      ksub[kmul[xx_,yy_],zz_] -> kmsub[xx,yy,zz],
+      ksub[zz_,kmul[xx_,yy_]] -> knmsub[xx,yy,zz],
+      kneg[kmadd [xx_,yy_,zz_]] -> knmadd[xx,yy,zz],
+      kneg[kmsub [xx_,yy_,zz_]] -> knmsub[xx,yy,zz],
+      kneg[knmadd[xx_,yy_,zz_]] -> kmadd [xx,yy,zz],
+      kneg[knmsub[xx_,yy_,zz_]] -> kmsub [xx,yy,zz]
+      (* we could match this and similar patterns
+      kmul[xx_, kadd[yy_, ToReal[+1]]] -> kmadd[xx, yy, xx],
+      kmul[xx_, kadd[yy_, ToReal[-1]]] -> kmsub[xx, yy, xx],
+      *)
+    };
+    expr = expr //. fmaRules;
+    
+    Return[expr]];
+
 (* Take an expression x and replace occurrences of Powers with the C
 macros SQR, CUB, QAD *)
-ReplacePowers[x_] :=
+ReplacePowers[expr_, vectorise_] :=
   Module[{rhs},
-    rhs = x /. Power[xx_, -1] -> INV[xx];
+    rhs = expr /. Power[xx_, -1] -> INV[xx];
 
     If[SOURCELANGUAGE == "C",
            Module[{},
-	     rhs = rhs //. Power[xx_, 2] -> SQR[xx];
-	     rhs = rhs //. Power[xx_, 3] -> CUB[xx];
-	     rhs = rhs //. Power[xx_, 4] -> QAD[xx];
+	     rhs = rhs /. Power[xx_,  2  ] -> SQR[xx];
+	     rhs = rhs /. Power[xx_,  3  ] -> CUB[xx];
+	     rhs = rhs /. Power[xx_,  4  ] -> QAD[xx];
+	     rhs = rhs /. Power[xx_, -2  ] -> INV[SQR[xx]];
+             rhs = rhs /. Power[xx_,  1/2] -> sqrt[xx];
+             rhs = rhs /. Power[xx_, -1/2] -> INV[sqrt[xx]];
+             rhs = rhs /. Power[xx_,  0.5] -> sqrt[xx];
+             rhs = rhs /. Power[xx_, -0.5] -> INV[sqrt[xx]];
+
+             (*
+	     rhs = rhs /.  1/2 ->  khalf
+	     rhs = rhs /. -1/2 -> -khalf;
 
-	     rhs = rhs //. xx_/2 -> khalf xx;
-	     rhs = rhs //. (-1/2) -> -khalf;
+	     rhs = rhs /.  1/3 ->  kthird;
+	     rhs = rhs /. -1/3 -> -kthird;
 
-	     rhs = rhs //. xx_/3 -> kthird xx;
-	     rhs = rhs //. (-1/3) -> -kthird;
+	     rhs = rhs /.  2/3 ->  ktwothird;
+	     rhs = rhs /. -2/3 -> -ktwothird;
 
-	     rhs = rhs //. 2/3 -> ktwothird;
-	     rhs = rhs //. (-2/3) -> -ktwothird;
+	     rhs = rhs /.  4/3 ->  kfourthird;
+	     rhs = rhs /. -4/3 -> -kfourthird;
 
-	     rhs = rhs //. 4/3 -> kfourthird;
-	     rhs = rhs //. (-4/3) -> -kfourthird;
+	     rhs = rhs /.  8/3 ->  keightthird;
+	     rhs = rhs /. -8/3 -> -keightthird;
+             *)
 
-	     rhs = rhs //. 8/3 -> keightthird;
-	     rhs = rhs //. (-8/3) -> -keightthird;
+             (* Avoid rational numbers *)
+	     rhs = rhs /. Rational[xx_,yy_] :> N[xx/yy, 30];
 
-	     rhs = rhs //. xx_ y_ + xx_ z_ -> xx(y+z); 
+             rhs = rhs //.      IfThen[cond1_,xx1_,yy1_] +      IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[    xx1 +     xx2], Simplify[    yy1 +     yy2]];
+             rhs = rhs //. ff1_ IfThen[cond1_,xx1_,yy1_] +      IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[ff1 xx1 +     xx2], Simplify[ff1 yy1 +     yy2]];
+             rhs = rhs //.      IfThen[cond1_,xx1_,yy1_] + ff2_ IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[    xx1 + ff2 xx2], Simplify[    yy1 + ff2 yy2]];
+             rhs = rhs //. ff1_ IfThen[cond1_,xx1_,yy1_] + ff2_ IfThen[cond2_,xx2_,yy2_] /; cond1==cond2 :> IfThen[cond1, Simplify[ff1 xx1 + ff2 xx2], Simplify[ff1 yy1 + ff2 yy2]];
 
-             rhs = rhs //. Power[E, power_] -> exp[power];
-             rhs = rhs //. Power[xx_, 0.5] -> sqrt[xx];
+             (* Is this still a good idea when FMA instructions are used? *)
+	     rhs = rhs //. xx_ yy_ + xx_ zz_ -> xx (yy+zz);
+	     rhs = rhs //. xx_ yy_ - xx_ zz_ -> xx (yy-zz);
+
+             rhs = rhs /. Power[E, power_] -> exp[power];
 
              (* there have been some problems doing the Max/Min
                 replacement via the preprocessor for C, so we do it
                 here *)
-             rhs = rhs //. Max[xx_, yy_] -> fmax[xx, yy];
-             rhs = rhs //. Min[xx_, yy_] -> fmin[xx, yy];
+             rhs = rhs /. Max[xx_, yy_] -> fmax[xx, yy];
+             rhs = rhs /. Min[xx_, yy_] -> fmin[xx, yy];
+
+             rhs = rhs /. Power[xx_, power_] -> pow[xx, power];
 
-             rhs = rhs //. Power[xx_, power_] -> pow[xx, power]],
+             If[vectorise === True,
+              rhs = vectoriseExpression[rhs]];
+           ],
 
-           rhs = rhs //. Power[xx_, power_] -> xx^power
+           rhs = rhs /. Power[xx_, power_] -> xx^power
        ];
 (*       Print[rhs//FullForm];*)
     rhs
@@ -841,180 +993,6 @@ CFormHideStrings[x_, opts___] := StringReplace[ToString[CForm[x,opts]], "\"" ->
 
 
 
-(* Debug output *)
-DebugCSE = True;
-
-(* Eliminate common subexpressions in a code sequence *)
-CSE[code_] := Module[
-  {expr, optexpr,
-   decomposed, locals, block,
-   block1, block2, temps1, stmts1, stmts2, stmts3,
-   replacevar,
-   stmts4,
-   stmts5, stmts6, stmts7},
-  if [DebugCSE, Print["code\n", code, "\nendcode\n"]];
-  
-  (* The code is passed in as list of {lhs,rhs} tuples.  Turn this
-     list into a single expression, so that it can be optimised.  *)
-  expr = code //. {a_, b__} -> CSequence[a, {b}]
-              //. {a_} -> a
-              //. (a_ -> b_) -> CAssign[a, b];
-  If [DebugCSE, Print["expr\n", expr, "\nendexpr\n"]];
-  
-  (* Optimise this expression *)
-  optexpr = Experimental`OptimizeExpression[expr];
-  If [DebugCSE, Print["optexpr\n", optexpr, "\nendoptexpr\n"]];
-  
-  (* This expression is a Mathematica expression.  Decompose it into
-     the set of newly introduced local variables and the optimised
-     expression itself.  *)
-  decomposed =
-    ReleaseHold[(Hold @@ optexpr)
-                /. Verbatim[Block][vars_, seq_] :> {vars, Hold[seq]}];
-  
-  If[decomposed[[0]] =!= List,
-    (* If the optimiser didn't create a Block expression, we assume it
-       didn't do anything useful and return the original.  *)
-    code,
-    
-    {locals, block} = decomposed;
-    If [DebugCSE, Print["locals\n", locals, "\nendlocals\n"]];
-    If [DebugCSE, Print["block\n", block, "\nendblock\n"]];
-    
-    block1 = block /. Hold[CompoundExpression[seq__]] :> Hold[{seq}];
-    If [DebugCSE, Print["block1\n", block1, "\nendblock1\n"]];
-    block2 = First[block1 //. Hold[{a___Hold, b_, c___}]
-                              /; Head[Unevaluated[b]] =!= Hold
-                              :> Hold[{a, Hold[b], c}]];
-    If [DebugCSE, Print["block2\n", block2, "\nendblock2\n"]];
-    
-    (* Temporaries, including a fake declaration for them *)
-    temps1 = Most[block2] //. Hold[lhs_ = rhs_] -> CAssign[CDeclare[lhs], rhs];
-    If [DebugCSE, Print["temps1\n", temps1, "\nendtemps1\n"]];
-    
-    (* Expression *)
-    stmts1 = ReleaseHold[Last[block2]];
-    If [DebugCSE, Print["stmts1\n", stmts1, "\nendstmts1\n"]];
-    
-    (* Turn CSequence back into a list *)
-    stmts2 = Flatten[{stmts1} //. CSequence[a_,b_] -> {a,b}];
-    If [DebugCSE, Print["stmts2\n", stmts2, "\nendstmts2\n"]];
-    
-    (* Combine temporaries and expression *)
-    stmts3 = Join[temps1, stmts2];
-    If [DebugCSE, Print["stmts3\n", stmts3, "\nendstmts3\n"]];
-    
-    (* Replace the internal names of the newly generated temporaries
-       with legal C names *)
-    replacevar =
-      Rule @@@ Transpose[{(*ToString[CForm[#]] & /@*) locals,
-        Symbol[
-          StringReplace[StringReplace[ToString[#], {__ ~~ "`" ~~ a_ :> a}],
-            "$" -> "T"]] & /@ locals}];
-    If [DebugCSE, Print["replacevar\n", replacevar, "\nendreplacevar\n"]];
-    
-    stmts4 = stmts3 //. replacevar;
-    If [DebugCSE, Print["stmts4\n", stmts4, "\nendstmts4\n"]];
-    
-    (* Sort statements topologically *)
-(*
-    stmts5 = stmts4;
-*)
-    If [DebugCSE, Print["A\n"]];
-    stmts5 =
-      Module[{debug,
-              tmpVars, newVars, i,
-              stmtsLeft, stmtsDone,
-              lhs, rhs, any, contains, containsAny,
-              canDoStmts, cannotDoStmts,
-              selfStmts, selfVars, allVars, nonSelfVars},
-        debug = False;
-        stmtsLeft = stmts4;
-        If [DebugCSE, Print["B\n"]];
-        stmtsDone = {};
-        If [DebugCSE, Print["C\n"]];
-        (* lhs[x_] := x[[1]]; *)
-        lhs[x_] := x /. (CAssign[lhs_, rhs_] -> lhs);
-        If [DebugCSE, Print["D\n"]];
-        (* rhs[x_] := x[[2]]; *)
-        rhs[x_] := x /. (CAssign[lhs_, rhs_] -> rhs);
-        If [DebugCSE, Print["E\n"]];
-        (* any[xs_] := Fold[Or, False, xs]; *)
-        any[xs_] := MemberQ[xs, True];
-        If [DebugCSE, Print["F\n"]];
-        (* contains[e_, x_] := (e /. x -> {}) =!= e; *)
-        (* contains[e_, x_] := Count[{e}, x, Infinity] > 0; *)
-        contains[e_, x_] := MemberQ[{e}, x, Infinity];
-        If [DebugCSE, Print["G\n"]];
-        containsAny[e_, xs_] := any[Map[contains[e,#]&, xs]];
-        If [DebugCSE, Print["H\n"]];
-        getVars[stmts_] := Map[lhs, stmts] //. (CDeclare[lhs_] -> lhs);
-
-        (* Rename temporary variables deterministically *)
-        tmpVars = Select[getVars[stmtsLeft],
-                         StringMatchQ[ToString[#], "TT"~~__]&];
-        newVars = Table[Symbol["T"<>ToString[1000000+i]],
-                        {i, 1, Length[tmpVars]}];
-        stmtsLeft = stmtsLeft /. MapThread[(#1->#2)&, {tmpVars, newVars}];
-
-        allVars = getVars[stmtsLeft];
-        While[stmtsLeft =!= {},
-          If[debug, Print["stmtsLeft = \n", stmtsLeft]];
-          If[debug, Print["stmtsDone = \n", stmtsDone]];
-          allVars = getVars[stmtsLeft];
-          If[debug, Print["allVars = \n", allVars]];
-          canDoStmts =
-            Select[stmtsLeft, Not[containsAny[rhs[#], allVars]] &];
-          cannotDoStmts =
-            Select[stmtsLeft,     containsAny[rhs[#], allVars]  &];
-          If[debug, Print["canDoStmts = \n", canDoStmts]];
-          If[debug, Print["cannotDoStmts = \n", cannotDoStmts]];
-          If[False && canDoStmts == {},
-            (* Handle assignment where LHS and RHS access the same variables
-               (hopefully without taking derivatives!) *)
-            selfStmts = Select[stmtsLeft, contains[rhs[#], lhs[#]]];
-            selfVars = getVars[selfStmts];
-            nonSelfVars = Select[allVars, Not[contains[selfVars, #]] &];
-            canDoStmts =
-              Select[stmtsLeft, Not[containsAny[rhs[#], nonSelfVars]] &];
-            cannotDoStmts =
-              Select[stmtsLeft,     containsAny[rhs[#], nonSelfVars]  &];
-            If[debug, Print["nonself/canDoStmts = \n", canDoStmts]];
-            If[debug, Print["nonself/cannotDoStmts = \n", cannotDoStmts]];
-          ];
-          If[canDoStmts == {},
-            (* Accept the first statement *)
-            canDoStmts = {First[stmtsLeft]};
-            cannotDoStmts = Rest[stmtsLeft];
-            If[debug, Print["takeone/canDoStmts = \n", canDoStmts]];
-            If[debug, Print["takeone/cannotDoStmts = \n", cannotDoStmts]];
-          ];
-          If[canDoStmts == {}, ThrowError["canDoStmts == {}"]];
-          stmtsDone = Join[stmtsDone, canDoStmts];
-          If [DebugCSE, Print["I\n"]];
-          stmtsLeft = cannotDoStmts;
-          If [DebugCSE, Print["J\n"]];
-        ];
-        If[debug, Print["stmtsLeft\n", stmtsLeft]];
-        If[debug, Print["stmtsDone\n", stmtsDone]];
-      stmtsDone];
-    If [DebugCSE, Print["Z\n"]];
-    
-    (* Turn CAssign statements back into (->) tuples *)
-    stmts6 = stmts5 //. CAssign[lhs_,rhs_] -> (lhs -> rhs);
-    If [DebugCSE, Print["stmts6\n", stmts6, "\nendstmts6\n"]];
-    
-    (* Turn CDeclare statements into "faked" declarations *)
-    stmts7 = stmts6
-      //. CDeclare[var_]
-          :> "CCTK_REAL const " <>
-             StringReplace[ToString[var], __ ~~ "`" -> ""];
-    If [DebugCSE, Print["stmts7\n", stmts7, "\nendstmts7\n"]];
-    
-    stmts7
-  ]
-];
-
 Quote[x_] := {"\"", x, "\""};
 
 End[];
diff --git a/Tools/CodeGen/Differencing.m b/Tools/CodeGen/Differencing.m
index ff6db9c..ab3d478 100644
--- a/Tools/CodeGen/Differencing.m
+++ b/Tools/CodeGen/Differencing.m
@@ -142,6 +142,7 @@ DZero::usage = "";
 shift::usage = "";
 spacing::usage = "";
 ComponentDerivativeOperatorStencilWidth::usage = "";
+CheckStencil::usage = "";
 
 Begin["`Private`"];
 
@@ -155,7 +156,7 @@ DZero[n_] := (DPlus[n] + DMinus[n])/2;
 (* User API *)
 (*************************************************************)
 
-CreateDifferencingHeader[derivOps_, zeroDims_] :=
+CreateDifferencingHeader[derivOps_, zeroDims_, vectorise_] :=
   Module[{componentDerivOps, dupsRemoved, expressions, componentDerivOps2, zeroDimRules, derivOps2, pDefs},
     Map[DerivativeOperatorVerify, derivOps];
 
@@ -167,14 +168,12 @@ CreateDifferencingHeader[derivOps_, zeroDims_] :=
 
     dupsRemoved = RemoveDuplicateRules[componentDerivOps2];
 
-    mDefPairs = Map[ComponentDerivativeOperatorMacroDefinition, dupsRemoved];
+    mDefPairs = Map[ComponentDerivativeOperatorMacroDefinition[#, vectorise] &, dupsRemoved];
 
     pDefs = Union[Flatten[Map[First, mDefPairs]]];
     expressions = Flatten[Map[#[[2]]&, mDefPairs]];
 
-(*    expressions = Flatten[Map[ComponentDerivativeOperatorInlineDefinition, dupsRemoved]];*)
-
-    {pDefs,Map[{#, "\n"} &, expressions]}];
+    {pDefs, Map[{#, "\n"} &, expressions]}];
 
 ordergfds[_[v1_,___], _[v2_,___]] := 
   Order[v1,v2] != -1;
@@ -205,6 +204,16 @@ ReplaceDerivatives[derivOps_, expr_, precompute_] :=
       rules = Map[# :> evaluateDerivative[#] &, gfds]];
     expr /. rules];
 
+(* Generate code to ensure that there are sufficient ghost and
+   boundary points for the passed derivative operators used in eqs *)
+CheckStencil[derivOps_, eqs_, name_] :=
+  Module[{gfds, rgzList, rgz},
+    gfds = Map[GridFunctionDerivativesInExpression[{#}, eqs] &, derivOps];
+    rgzList = MapThread[If[Length[#2] > 0, DerivativeOperatorStencilWidth[#1], {0,0,0}] &, {derivOps, gfds}];
+    If[Length[rgzList] === 0, Return[{}]];
+    rgz = Map[Max, Transpose[rgzList]];
+    If[Max[rgz] == 0, {},
+    {"GenericFD_EnsureStencilFits(cctkGH, ", Quote@name, ", ", Riffle[rgz,", "], ");\n"}]];
 
 (*************************************************************)
 (* Misc *)
@@ -212,15 +221,18 @@ ReplaceDerivatives[derivOps_, expr_, precompute_] :=
 
 PrecomputeDerivative[d:pd_[gf_, inds___]] :=
   Module[{},
-    DeclareAssignVariable["CCTK_REAL", GridFunctionDerivativeName[d], evaluateDerivative[d]]];
+    DeclareAssignVariable[DataType[], GridFunctionDerivativeName[d], evaluateDerivative[d]]];
 
 evaluateDerivative[d:pd_[gf_, inds___]] :=
   Module[{macroname},
     macroName = ComponentDerivativeOperatorMacroName[pd[inds] -> expr];
-    Return[ToString[macroName] <> "(" <> ToString[gf] <> ", i, j, k)"]];
+    (* Return[ToString[macroName] <> "(" <> ToString[gf] <> ", i, j, k)"] *)
+    (* Return[ToString[macroName] <> "(" <> ToString[gf] <> ")"] *)
+    Return[ToString[macroName] <> "(&" <> ToString[gf] <> "[index])"]
+  ];
 
 DeclareDerivative[d:pd_[gf_, inds___]] :=
-  DeclareVariable[GridFunctionDerivativeName[d], "// CCTK_REAL"];
+  DeclareVariable[GridFunctionDerivativeName[d], "// CCTK_REAL_VEC"];
 
 
 (*************************************************************)
@@ -251,8 +263,8 @@ sbpMacroDefinition[macroName_, d_] :=
     FlattenBlock[{"#define ", macroName, "(u,i,j,k) (sbp_deriv_" <> ds
     <> "(i,j,k,sbp_" <> l <> "min,sbp_" <> l <> "max,d" <> ds <> ",u,q" <> ds <> ",cctkGH))"}]    ];
 
-ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> expr_)] :=
-  Module[{macroName, rhs, rhs2, i = "i", j = "j", k = "k", spacings, spacings2, pat, ss, num, den, newnum, signModifier, quotient, liName, rhs3, rhs4},
+ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> expr_), vectorise_] :=
+  Module[{macroName, rhs, i = "i", j = "j", k = "k", spacings, spacings2, pat, ss, num, den, newnum, signModifier, quotient, liName, finalDef},
   
     macroName = ComponentDerivativeOperatorMacroName[componentDerivOp];
 
@@ -265,23 +277,24 @@ ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> e
     If[expr === SBPDerivative[3],
       Return[sbpMacroDefinition[macroName, 3]]];
 
-    rhs = DifferenceGF[expr, i, j, k];
+    rhs = DifferenceGF[expr, i, j, k, vectorise];
+(*    Print["rhs1 == ", FullForm[rhs]];*)
     spacings = {spacing[1] -> 1/"dxi", spacing[2] -> 1/"dyi", spacing[3] -> 1/"dzi"};
     spacings2 = {spacing[1] -> "dx", spacing[2] -> "dy", spacing[3] -> "dz"};
 
-    rhs2 = FullSimplify[rhs];
+    rhs = FullSimplify[rhs];
 
-(*    Print["rhs2 == ", FullForm[rhs2]];*)
+(*    Print["rhs2 == ", FullForm[rhs]];*)
 
     pat = (Times[spInExpr:(Power[spacing[_],_]...), (Rational[x_,y_])..., rest__]) | (rest__);
     (* Print["pat == ", pat//FullForm]; *)
 
-    If[MatchQ[rhs2, pat],
+    If[MatchQ[rhs, pat],
 (*       Print["matches!"];*)
-       ss = Times[rhs2 /. pat -> spInExpr];
+       ss = Times[rhs /. pat -> spInExpr];
 (*       Print["ss == ", ss];*)
-       num = rhs2 /. pat -> x;
-       den = rhs2 /. pat -> y;
+       num = rhs /. pat -> x;
+       den = rhs /. pat -> y;
 (*       Print["num == ", num];
        Print["den == ", den];*)
        If[{num, 1, 2} === {1, 2},(* Print["SEQ!"]; *) newnum = 1; den=1; signModifier = "",
@@ -307,84 +320,117 @@ ComponentDerivativeOperatorMacroDefinition[componentDerivOp:(name_[inds___] -> e
        liName = "p" <> signModifier <> quotient <> ToString[Apply[SequenceForm,Simplify[1/(ss /. spacings2)],{0,Infinity}]];
 (*       Print["liName == ", liName];*)
 
-       rhs3 = rhs2 /. pat -> Times[liName, rest],
-       ThrowError["Partial derivative operator definition ", rhs2, 
-         " does not match pattern ", pat];
-       rhs3 = rhs2];
-
-(*    Print["rhs3 == ", rhs3];*)
+       rhs = rhs /. pat -> Times[liName, rest],
+(*
+       rhs = (rhs /. pat -> Times[liName, rest]) / If[vectorise, liName, 1], (* horrible *)
+*)
+(*       Print["!!!!!!!!DOES NOT MATCH!!!!!!!!!"];*)
+       rhs = rhs];
 
-    pDefs = {{liName -> CFormHideStrings[ReplacePowers[num / den ss /. spacings2]]}};
+(*    Print["rhs3 == ", FullForm[rhs]];*)
 
-(*    rhs4 = Factor[rhs3];*)
+    pDefs = {{liName -> CFormHideStrings[ReplacePowers[num / den ss /. spacings2, vectorise]]}};
 
-    rhs4 = rhs3 //. (x_ a_ + x_ b_) -> x(a+b);
-    rhs5 = rhs4 //. (x_ a_ - x_ b_) -> x(a-b);
+(*    rhs = Factor[rhs];*)
+    rhs = rhs //. (x_ a_ + x_ b_) -> x (a+b);
+    rhs = rhs //. (x_ a_ - x_ b_) -> x (a-b);
     
 (*    Print[componentDerivOp, ": "];
-    Print[FullForm[rhs5]];
+    Print[FullForm[rhs]];
     Print[""];*)
 
-    rhs6 = CFormHideStrings[ReplacePowers[rhs5 /. spacings]];
-    {pDefs, FlattenBlock[{"#define ", macroName, "(u,i,j,k) ", "(", rhs6, ")"}]}];
-
-ComponentDerivativeOperatorInlineDefinition[componentDerivOp:(name_[inds___] -> expr_)] :=
-  Module[{inlineName, rhs, rhs2, i = "i", j = "j", k = "k", spacings},
-  
-    inlineName = ComponentDerivativeOperatorMacroName[componentDerivOp];
-
-    rhs = DifferenceGF[expr, i, j, k];
-(*    rhs = DifferenceGFInline[expr, i, j, k];*)
-    spacings = {spacing[1] -> 1/"dxi", spacing[2] -> 1/"dyi", spacing[3] -> 1/"dzi"};
-    rhs2 = CFormHideStrings[FullSimplify[ReplacePowers[rhs /. spacings]]];
-
-    DefineFunction[inlineName, "static inline CCTK_REAL", 
-      "CCTK_REAL *u, int i, int j, int k",
-      {"return ", rhs2, ";\n"}]];
+    rhs = CFormHideStrings[ReplacePowers[rhs /. spacings, vectorise]];
+    (* Print["rhs=",FullForm[rhs]]; *)
+    (* {pDefs, FlattenBlock[{"#define ", macroName, "(u,i,j,k) ", "(", rhs, ")"}]} *)
+    finalDef =
+      If[vectorise,
+    {pDefs, FlattenBlock[{
+      "#ifndef KRANC_DIFF_FUNCTIONS\n",
+       (* default, differencing operators are macros *)
+(*
+      "#  define ", macroName, "(u) ", "(kmul(", liName, ",", rhs, "))\n",
+*)
+      "#  define ", macroName, "(u) ", "(", rhs, ")\n",
+      "#else\n",
+       (* new, differencing operators are static functions *)
+(*
+      "#  define ", macroName, "(u) ", "(kmul(", liName, ",", macroName, "_impl((u),dj,dk)))\n",
+      "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, ptrdiff_t const dj, ptrdiff_t const dk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;\n",
+      "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, ptrdiff_t const dj, ptrdiff_t const dk)\n",
+      If[StringMatchQ[rhs, RegularExpression[".*\\bdir\\d\\b.*"]],
+         { "{ return ToReal(1e30); /* ERROR */ }\n" },
+         { "{ return ", rhs, "; }\n" }],
+*)
+(*
+      "#  define ", macroName, "(u) ", "(kmul(", liName, ",", macroName, "_impl(u,cdj,cdk)))\n",
+      "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;\n",
+      "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, ptrdiff_t const cdj, ptrdiff_t const cdk)\n",
+      (* We cannot handle dirN,
+         so we punt on all expressions that contain dirN *)
+      If[StringMatchQ[rhs, RegularExpression[".*\\bdir\\d\\b.*"]],
+         { "{ return ToReal(1e30); /* ERROR */ }\n" },
+         { "{\n",
+           "  ptrdiff_t const cdi=sizeof(CCTK_REAL);\n",
+           "  return ", rhs, ";\n",
+           "}\n" }],
+*)
+      "#  define ", macroName, "(u) ", "(", macroName, "_impl(u,", liName, ",cdj,cdk))\n",
+      "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, CCTK_REAL_VEC const ", liName, ", ptrdiff_t const cdj, ptrdiff_t const cdk) CCTK_ATTRIBUTE_NOINLINE CCTK_ATTRIBUTE_UNUSED;\n",
+      "static CCTK_REAL_VEC ", macroName, "_impl(CCTK_REAL const* restrict const u, CCTK_REAL_VEC const ", liName, ", ptrdiff_t const cdj, ptrdiff_t const cdk)\n",
+      (* We cannot handle dirN,
+         so we punt on all expressions that contain dirN *)
+      If[StringMatchQ[rhs, RegularExpression[".*\\bdir\\d\\b.*"]],
+         { "{ return ToReal(1e30); /* ERROR */ }\n" },
+         { "{\n",
+           "  ptrdiff_t const cdi=sizeof(CCTK_REAL);\n",
+           "  return ", rhs, ";\n",
+           "}\n" }],
+      "#endif\n"
+    }]},
+
+    {pDefs, FlattenBlock[{"#define ", macroName, "(u) ", "(", rhs, ")"}]}];
+    finalDef
+];
 
 ComponentDerivativeOperatorMacroName[componentDerivOp:(name_[inds___] -> expr_)] :=
   Module[{stringName},
     stringName = StringJoin[Map[ToString, Join[{name}, {inds}]]];
     stringName];
 
-
-
+DerivativeOperatorStencilWidth[derivOp_] :=
+  Map[Max, Transpose[Map[ComponentDerivativeOperatorStencilWidth, DerivativeOperatorToComponents[derivOp]]]];
 
 ComponentDerivativeOperatorStencilWidth[componentDerivOp:(name_[inds___] -> expr_)] :=
-  Module[{cases, nx, ny, nz},
-    cases = Union[Flatten[Cases[{expr}, shift[_] | Power[shift[_],_], Infinity]]];
-    Print[cases];
-
-    nx = Exponent[op, shift[1]];
-    ny = Exponent[op, shift[2]];
-    nz = Exponent[op, shift[3]];
-
-  ];
-
-
-
+  Module[{cases, nx, ny, nz, result},
+    result = Table[
+    cases = Union[Flatten[Cases[{expr}, shift[d] | Power[shift[d],_], Infinity]]];
+    ns = Map[Exponent[#, shift[d]] &, cases];
+    If[Length[ns] == 0, 0, Max[Abs[ns]]], {d, 1, 3}];
+
+    (* We do not know the run-time value of any shorthands used in
+       operator definitions.  In all the current known cases, this
+       will be a "direction" which is +/- 1.  In future, the
+       differencing mechanism will support shorthand arguments to
+       operators and this hack can be removed. *)
+    result = Replace[result, _Symbol -> 1, {-1}];
+
+    If[!And@@Map[NumericQ, result],
+      Throw["Stencil width is not numeric in "<>ToString[componentDerivOp]]];
+    result];
 
 (* Farm out each term of a difference operator *)
-DifferenceGF[op_, i_, j_, k_] :=
-  Module[{expanded},
-    expanded = Expand[op];
-    
-    If[Head[expanded] === Plus,
-      Apply[Plus, Map[DifferenceGFTerm[#, i, j, k] &, expanded]],
-      DifferenceGFTerm[expanded, i, j, k]]];
-
-DifferenceGFInline[op_, i_, j_, k_] :=
+DifferenceGF[op_, i_, j_, k_, vectorise_] :=
   Module[{expanded},
     expanded = Expand[op];
     
     If[Head[expanded] === Plus,
-      Apply[Plus, Map[DifferenceGFTermInline[#, i, j, k] &, expanded]],
+      Apply[Plus, Map[DifferenceGFTerm[#, i, j, k, vectorise] &, expanded]],
       DifferenceGFTerm[expanded, i, j, k]]];
 
 
 (* Return the fragment of a macro definition for defining a derivative
    operator *)
-DifferenceGFTerm[op_, i_, j_, k_] :=
+DifferenceGFTerm[op_, i_, j_, k_, vectorise_] :=
   Module[{nx, ny, nz, remaining},
 
     If[op === 0,
@@ -409,10 +455,52 @@ DifferenceGFTerm[op_, i_, j_, k_] :=
       "(int)(" <> ToString[CFormHideStrings[j+ny]] <> ")," <>
       "(int)(" <> ToString[CFormHideStrings[k+nz]] <> "))]",
 *)
-    remaining "(u)[index" <>
-      "+di*(" <> ToString[CFormHideStrings[nx]] <> ")" <>
+(*
+    remaining "vec_loadu_maybe" <>
+      "(" <> ToString[CFormHideStrings[nx]] <> "," <>
+             "(u)[index" <>
+                  "+di*(" <> ToString[CFormHideStrings[nx]] <> ")" <>
+                  "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <>
+                  "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")])",
+*)
+(*
+    remaining "vec_loadu_maybe(" <> ToString[CFormHideStrings[nx]] <> "," <>
+      "(u)[(" <> ToString[CFormHideStrings[nx]] <> ")" <>
+      "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <>
+      "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")])",
+*)
+
+  If[vectorise,
+(*
+    remaining "vec_loadu_maybe3" <>
+      "(" <> ToString[CFormHideStrings[nx /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
+             ToString[CFormHideStrings[ny /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
+             ToString[CFormHideStrings[nz /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
+      "(u)[(" <> ToString[CFormHideStrings[nx]] <> ")" <>
+      "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <>
+      "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")])",
+*)
+    remaining "vec_loadu_maybe3" <>
+      "(" <> ToString[CFormHideStrings[nx /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
+             ToString[CFormHideStrings[ny /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
+             ToString[CFormHideStrings[nz /. {dir1->1, dir2->1, dir3->1}]] <> "," <>
+      "*(CCTK_REAL const*)&((char const*)(u))" <>
+        "[cdi*(" <> ToString[CFormHideStrings[nx]] <> ")" <>
+        "+cdj*(" <> ToString[CFormHideStrings[ny]] <> ")" <>
+        "+cdk*(" <> ToString[CFormHideStrings[nz]] <> ")])",
+
+    remaining "(u)[" <>
+      "di*(" <> ToString[CFormHideStrings[nx]] <> ")" <>
       "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <>
-      "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")]",
+      "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")]"],
+
+
+(*
+    remaining "vec_loadu" <>
+      "(u[(" <> ToString[CFormHideStrings[nx]] <> ")" <>
+                "+dj*(" <> ToString[CFormHideStrings[ny]] <> ")" <>
+                "+dk*(" <> ToString[CFormHideStrings[nz]] <> ")])",
+*)
 (*
     remaining "(u)[CCTK_GFINDEX3D(cctkGH,floor((" <>
       ToString[CFormHideStrings[i+nx]] <> ")+0.5),floor((" <>
@@ -422,27 +510,6 @@ DifferenceGFTerm[op_, i_, j_, k_] :=
     remaining "u(" <> ToString[FortranForm[i+nx]] <> "," <> 
       ToString[FortranForm[j+ny]] <> "," <> ToString[FortranForm[k+nz]] <> ")"] ];
 
-(* Return the fragment of a function definition for defining a derivative
-   operator *)
-DifferenceGFTermInline[op_, i_, j_, k_] :=
-  Module[{nx, ny, nz, remaining},
-
-    If[op === 0,
-      Return[0]];
-
-    nx = Exponent[op, shift[1]];
-    ny = Exponent[op, shift[2]];
-    nz = Exponent[op, shift[3]];
-
-    remaining = op / (shift[1]^nx) / (shift[2]^ny) / (shift[3]^nz);
-
-    If[Cases[{remaining}, shift[_], Infinity] != {},
-      ThrowError["Could not parse difference operator:", op]];
-    
-    remaining "(u)[CCTK_GFINDEX3D(cctkGH," <> ToString[CFormHideStrings[i+nx]] <> "," <>
-      ToString[CFormHideStrings[j+ny]] <> "," <> ToString[CFormHideStrings[k+nz]] <> ")]"
-     ];
-
 
 DerivativeOperatorGFDs[gf_];
 
diff --git a/Tools/CodeGen/Interface.m b/Tools/CodeGen/Interface.m
index 29c198f..5f8f613 100644
--- a/Tools/CodeGen/Interface.m
+++ b/Tools/CodeGen/Interface.m
@@ -73,7 +73,7 @@ CreateKrancInterface[nonevolvedGroups_, evolvedGroups_, rhsGroups_, groups_,
 
   Module[{registerEvolved, (*registerConstrained,*)
     nonevolvedGroupStructures, evolvedGroupStructures, rhsGroupStructures,
-    groupStructures, interface},
+    groupStructures, interface, getMap},
     VerifyGroupNames[nonevolvedGroups];
     VerifyGroupNames[evolvedGroups];
     VerifyGroupNames[rhsGroups];
@@ -105,6 +105,13 @@ CreateKrancInterface[nonevolvedGroups_, evolvedGroups_, rhsGroups_, groups_,
       ArgString -> "CCTK_POINTER_TO_CONST IN cctkGH, CCTK_INT IN dir, CCTK_INT IN nsize, CCTK_INT OUT ARRAY imin, CCTK_INT OUT ARRAY imax, CCTK_REAL OUT ARRAY q, CCTK_INT IN table_handle"
     };
 
+    getMap = 
+    {
+      Name -> "MultiPatch_GetMap",
+      Type -> "CCTK_INT",
+      ArgString -> "CCTK_POINTER_TO_CONST IN cctkGH"
+    };
+
 
     (* For each group declared in this thorn, we need an entry in the
         interface file.  Each evolved group needs an associated rhs
@@ -127,10 +134,11 @@ CreateKrancInterface[nonevolvedGroups_, evolvedGroups_, rhsGroups_, groups_,
 
     interface = CreateInterface[implementation, inheritedImplementations, 
       Join[includeFiles, {CactusBoundary`GetIncludeFiles[]},
-           If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}]],
+           If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}],
+           If[OptionValue[UseVectors], {"vectors.h"}, {}]],
       groupStructures,
       UsesFunctions ->
-        Join[{registerEvolved, (*registerConstrained,*) diffCoeff}, 
+        Join[{registerEvolved, (*registerConstrained,*) diffCoeff, getMap}, 
              CactusBoundary`GetUsedFunctions[]]];
     Return[interface]];
 
diff --git a/Tools/CodeGen/Jacobian.m b/Tools/CodeGen/Jacobian.m
new file mode 100644
index 0000000..31f6b67
--- /dev/null
+++ b/Tools/CodeGen/Jacobian.m
@@ -0,0 +1,140 @@
+
+(*  Copyright 2011 Ian Hinder
+
+    This file is part of Kranc.
+
+    Kranc is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    Kranc is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Kranc; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*)
+
+BeginPackage["Jacobian`", {"Errors`", "Helpers`", "Kranc`", "Differencing`", "MapLookup`", "CodeGen`", "KrancGroups`"}];
+
+JacobianQ;
+InsertJacobian;
+CreateJacobianVariables;
+JacobianGenericFDParameters;
+JacobianSymbols;
+JacobianGroups;
+JacobianCheckGroups;
+JacobianConditionalGridFunctions;
+
+Begin["`Private`"];
+
+Options[JacobianQ] = ThornOptions;
+JacobianQ[opts:OptionsPattern[]] :=
+  Length[OptionValue[Jacobian]] > 0;
+
+(* Assign a shorthand containing the Jacobian multiplied by the passed
+   1st derivative *)
+jacobianShorthand[d:(deriv_[var_, i_])] :=
+  Module[{},
+     derivToJacDeriv[d] ->
+      IfThen["use_jacobian", Sum[Symbol["J"<>ToString[j]<>ToString[i]] deriv[var, j], {j, 1 3}], deriv[var, i]]
+  ];
+
+(* Assign a shorthand containing the Jacobian multiplied by the passed
+   2nd derivative *)
+jacobianShorthand[d:(deriv_[var_, i_,j_])] :=
+  Module[{ip,jp},
+     {ip,jp} = Sort[{i,j}]; (* dJ is symmetric in the last two indices *)
+     derivToJacDeriv[d] ->
+      IfThen["use_jacobian", Sum[Symbol["dJ"<>ToString[a]<>ToString[ip]<>ToString[jp]] deriv[var, a], {a, 1 3}] + 
+      Sum[Symbol["J"<>ToString[a]<>ToString[i]] Symbol["J"<>ToString[b]<>ToString[j]] deriv[var, a, b], {a, 1 3}, {b, 1, 3}],
+      deriv[var, i, j]]
+  ];
+
+(* Convert a 1st derivative to a Jacobian-multiplied derivative *)
+derivToJacDeriv[deriv_[var_, i_]] :=
+  Symbol["Global`Jac"<>ToString[deriv]<>ToString[i]<>ToString[var]];
+
+(* Convert a 2nd derivative to a Jacobian-multiplied derivative *)
+derivToJacDeriv[deriv_[var_, i_, j_]] :=
+  Symbol["Global`Jac"<>ToString[deriv]<>ToString[i]<>ToString[j]<>ToString[var]];
+
+(* Given a calculation containing partial derivatives, return a
+   version of the calculation with all the partial derivatives multiplied
+   by the Jacobian *)
+Options[InsertJacobian] = ThornOptions;
+InsertJacobian[calc_List, opts:OptionsPattern[]] :=
+  Module[{pdDefs, derivs, newShortDefs, newShorts, combinedShorts, combinedEqs, combinedCalc, eqs, newEqs},
+    pdDefs = OptionValue[PartialDerivatives];
+    derivs = GridFunctionDerivativesInExpression[pdDefs, lookup[calc, Equations]];
+    If[Length[derivs] === 0, Return[calc]];
+    newShortDefs = Map[jacobianShorthand, derivs];
+    newShorts = Map[First, newShortDefs];
+    combinedShorts = Join[lookupDefault[calc, Shorthands, {}], newShorts];
+    eqs = lookup[calc, Equations];
+    newEqs = eqs /. (x_?(MemberQ[derivs, #] &) :> derivToJacDeriv[x]);
+    combinedEqs = Join[newShortDefs, newEqs];
+    combinedCalc = mapReplace[mapReplace[mapEnsureKey[calc, Shorthands, {}], Shorthands, combinedShorts], Equations, combinedEqs];
+    combinedCalc];
+
+(* Define local pointers to the members of the Jacobian and Jacobian
+   derivatives groups *)
+CreateJacobianVariables[] :=
+CommentedBlock["Jacobian variable pointers",
+  {"bool const use_jacobian = (!CCTK_IsFunctionAliased(\"MultiPatch_GetMap\") || MultiPatch_GetMap(cctkGH) != jacobian_identity_map)\n                     && strlen(jacobian_group) > 0;\n",
+   "if (use_jacobian && strlen(jacobian_derivative_group) == 0)\n",
+   "{\n",
+   "  CCTK_WARN (1, \"GenericFD::jacobian_group and GenericFD::jacobian_derivative_group must both be set to valid group names\");\n",
+   "}\n\n",
+   "CCTK_REAL const *restrict jacobian_ptrs[9];\n",
+   "if (use_jacobian) GenericFD_GroupDataPointers(cctkGH, jacobian_group,\n",
+   "                                              9, jacobian_ptrs);\n",
+    "\n",
+    Table[{"CCTK_REAL const *restrict const J",i,j," = use_jacobian ? jacobian_ptrs[",(i-1)*3+j-1,"] : 0;\n"},{i,1,3},{j,1,3}],
+   "\n",
+   "CCTK_REAL const *restrict jacobian_derivative_ptrs[18];\n",
+   "if (use_jacobian) GenericFD_GroupDataPointers(cctkGH, jacobian_derivative_group,\n",
+   "                                              18, jacobian_derivative_ptrs);\n",
+    "\n",
+    Module[{syms = Flatten[Table[{"dJ",i,j,k},{i,1,3},{j,1,3},{k,j,3}],2]},
+      MapIndexed[{"CCTK_REAL const *restrict const ", #1, " = use_jacobian ? jacobian_derivative_ptrs[", #2-1, "] : 0;\n"} &, syms]]}];
+
+(* List of symbols which should be allowed in a calculation *)
+JacobianSymbols[] :=
+  Map[Symbol, Join[Flatten[Table[FlattenBlock[{"dJ",i,j,k}],{i,1,3},{j,1,3},{k,j,3}],2],
+    Flatten[Table[FlattenBlock[{"J",i,j}],{i,1,3},{j,1,3}],1]]];
+
+(* Parameters to inherit from GenericFD *)
+JacobianGenericFDParameters[] :=
+  {{Name -> "jacobian_group",     Type -> "CCTK_STRING"},
+   {Name -> "jacobian_derivative_group",     Type -> "CCTK_STRING"},
+   {Name -> "jacobian_identity_map",     Type -> "CCTK_INT"}};
+
+(* The symbols which are used for the Jacobian variables in the
+   generated source code.  These do not have to coincide with the
+   actual variable names, as the variable pointers are read using
+   CCTK_VarDataPtr. *)
+JacobianGroups[] :=
+  {{"unknown::unknown",  {Global`J11, Global`J12, Global`J13, Global`J21, Global`J22, Global`J23, Global`J31, Global`J32, Global`J33}},
+   {"unknown::unknown",  {Global`dJ111, Global`dJ112, Global`dJ113, Global`dJ122, Global`dJ123, Global`dJ133,
+                               Global`dJ211, Global`dJ212, Global`dJ213, Global`dJ222, Global`dJ223, Global`dJ233,
+                               Global`dJ311, Global`dJ312, Global`dJ313, Global`dJ322, Global`dJ323, Global`dJ333}}};
+
+JacobianCheckGroups[groups_] :=
+  Module[{int},
+    int = Intersection[allGroupVariables[groups], allGroupVariables[JacobianGroups[]]];
+    If[int =!= {},
+      Throw["Error: Some group variables conflict with reserved Jacobian variable names: " <> ToString[int]]]];
+
+(* These gridfunctions are only given local variable copies if the use_jacobian variable is true *)
+JacobianConditionalGridFunctions[] :=
+  {("dJ" ~~ DigitCharacter ~~ DigitCharacter ~~ DigitCharacter) | ("J" ~~ DigitCharacter ~~ DigitCharacter),
+   "use_jacobian",
+   None};
+
+End[];
+
+EndPackage[];
diff --git a/Tools/CodeGen/Kranc.m b/Tools/CodeGen/Kranc.m
index 2224d18..ffa98e7 100644
--- a/Tools/CodeGen/Kranc.m
+++ b/Tools/CodeGen/Kranc.m
@@ -22,7 +22,11 @@ BeginPackage["Kranc`"];
 
 (* CodeGen.m *)
 
-{INV, SQR, CUB, QAD, exp, pow, fmax, fmin, dx, dy, dz, khalf, kthird, ktwothird, kfourthird, keightthird};
+{INV, SQR, CUB, QAD, IfThen, ToReal, sqrt, exp, pow, fmax, fmin,
+ kmadd, kmsub, knmadd, knmsub, kpos, kneg, kadd, ksub, kmul, kdiv,
+ kfabs, kfmax, kfmin, ksqrt, kexp, klog, kpow,
+ dir1, dir2, dir3, dt, dx, dy, dz,
+ khalf, kthird, ktwothird, kfourthird, keightthird};
 
 (* Helpers.m *)
 
@@ -34,7 +38,7 @@ dummy;
 LoopPreIncludes, GroupImplementations, PartialDerivatives, NoSimplify,
 Boundary, Interior, InteriorNoSync, Where, AddToStencilWidth,
 Everywhere, normal1, normal2, normal3, INV, SQR, CUB, QAD, dot, pow,
-exp, dx, dy, dz, idx, idy, idz, MinMod, VanLeer}
+exp, dt, dx, dy, dz, idx, idy, idz, MinMod, VanLeer}
 
 {ConditionalOnKeyword, ConditionalOnKeywords, CollectList, Interior,
 InteriorNoSync, Boundary, BoundaryWithGhosts, Where, PreDefinitions,
@@ -67,9 +71,11 @@ ThornOptions =
   ReflectionSymmetries -> {},
   ZeroDimensions -> {},
   UseLoopControl -> False,
-  UseCSE -> False,
+  UseVectors -> False,
   ProhibitAssignmentToGridFunctionsRead -> False,
-  IncludeFiles -> {}};
+  IncludeFiles -> {},
+  CSE -> False,
+  UseJacobian -> False};
 
 (* Thorn.m *)
 
diff --git a/Tools/CodeGen/KrancGroups.m b/Tools/CodeGen/KrancGroups.m
index 05ebdea..13d53d2 100644
--- a/Tools/CodeGen/KrancGroups.m
+++ b/Tools/CodeGen/KrancGroups.m
@@ -52,6 +52,7 @@ AddGroupExtra;
 GroupTimelevels;
 allGroupVariables;
 NonevolvedTimelevels;
+CheckGroups;
 
 Begin["`Private`"];
 
@@ -245,6 +246,18 @@ qualifyGFName[gfname_, allgroups_, defaultImp_] :=
 allGroupVariables[groups_] :=
   Flatten[Map[groupVariables, groups], 1];
 
+CheckGroups[groups_] :=
+  Module[{vs, names},
+(*    If[!MatchQ[{_String, {_Symbol ...}, {_Symbol -> _} ...}],
+        ThrowError["Groups structure should be of the form {name, {vars, ...}, extras}"]]; *)
+
+    vs = Map[ToString, Union[Flatten[Map[groupVariables, groups]]]];
+    names = Map[groupName, groups];
+
+    If[(int = Intersection[vs,names]) =!= {},
+      ThrowError["Variable names and group names must be distinct.  Group names which are also variable names:", int]];
+  ];
+
 End[];
 
 EndPackage[];
diff --git a/Tools/CodeGen/KrancTensor.m b/Tools/CodeGen/KrancTensor.m
new file mode 100644
index 0000000..204bc60
--- /dev/null
+++ b/Tools/CodeGen/KrancTensor.m
@@ -0,0 +1,90 @@
+(* ::Package:: *)
+
+(*  Copyright 2004-2010
+    Sascha Husa, Ian Hinder, Christiane Lechner, Barry Wardell
+
+    This file is part of Kranc.
+
+    Kranc is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    Kranc is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Kranc; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*)
+
+(****************************************************************************)
+(* Wrapper providing tensor support to Kranc (from TensorTools or xTensor)  *)
+(****************************************************************************)
+If[!ValueQ[$KrancTensorPackage], $KrancTensorPackage = "TensorToolsKranc`"];
+
+BeginPackage["KrancTensor`", {"Errors`", "KrancThorn`", "MapLookup`", "KrancGroups`", "Kranc`", $KrancTensorPackage}];
+
+CreateKrancThornTT::usage = "Construct a Kranc thorn using tensor expressions.";
+
+(* FIXME: Move CreateGroupFromTensor here *)
+
+Begin["`Private`"];
+
+(* --------------------------------------------------------------------------
+   Tensor Tools
+   -------------------------------------------------------------------------- *)
+
+CreateKrancThornTT[groups_, parentDirectory_, thornName_, opts___] :=
+  Module[{calcs, expCalcs, expGroups, options, derivs, expDerivs, reflectionSymmetries, declaredGroups},
+    InfoMessage[Terse, "Processing tensorial arguments"];
+    calcs = lookup[{opts}, Calculations];
+    derivs = lookupDefault[{opts}, PartialDerivatives, {}];
+    Map[CheckCalculationTensors, calcs];
+    expCalcs = Map[makeCalculationExplicit, calcs];
+
+    InfoMessage[Info, "Group definitions:", groups];
+    VerifyGroups[groups];
+
+    expDerivs = Flatten[Map[ExpandComponents,derivs],1];
+    expGroups = Map[makeGroupExplicit, groups];
+    options = Join[DeleteCases[{opts}, Calculations -> _], {Calculations -> expCalcs}];
+    options = Join[DeleteCases[options, PartialDerivatives -> _], {PartialDerivatives -> expDerivs}];
+
+    declaredGroups = lookupDefault[{opts}, DeclaredGroups, {}];
+    evolutionTimelevels = lookupDefault[{opts}, EvolutionTimelevels, 3];
+    defaultEvolutionTimelevels = lookupDefault[{opts}, DefaultEvolutionTimelevels, evolutionTimelevels];
+    InfoMessage[Info, "Declared groups: " <> ToString[declaredGroups]];
+    InfoMessage[Terse, "Computing reflection symmetries"];
+    reflectionSymmetries = computeReflectionSymmetries[declaredGroups, groups];
+    InfoMessage[Info, "Reflection symmetries: ", reflectionSymmetries];
+
+    InfoMessage[Terse, "Creating (component-based) Kranc thorn"];
+
+    (* It is necessary to include the KrancThorn context here due to some annoying Needs[] dependency issue *)
+    KrancThorn`CreateKrancThorn[expGroups, parentDirectory, thornName,
+      Apply[Sequence, options], ReflectionSymmetries -> reflectionSymmetries]];
+
+computeReflectionSymmetries[declaredGroups_, groups_] :=
+  Module[{variables, syms},
+    variables = variablesFromGroups[declaredGroups, groups];
+    syms = Flatten[Map[ReflectionSymmetries, variables], 1];
+    syms];
+
+makeCalculationExplicit[calc_] :=
+  mapValueMapMultiple[calc, 
+    {Shorthands -> ExpandComponents,
+     CollectList -> ExpandComponents,
+     Equations -> ExpandComponents}];
+
+makeGroupExplicit[g_] :=
+  Module[{variables, newVariables, newGroup},
+    variables = groupVariables[g];
+    newVariables = DeleteDuplicates[ExpandComponents[variables]];
+    newGroup = SetGroupVariables[g, newVariables];
+    newGroup];
+
+End[];
+EndPackage[];
diff --git a/Tools/CodeGen/KrancThorn.m b/Tools/CodeGen/KrancThorn.m
index d41fb51..4d53ccb 100644
--- a/Tools/CodeGen/KrancThorn.m
+++ b/Tools/CodeGen/KrancThorn.m
@@ -1,3 +1,4 @@
+(* ::Package:: *)
 
 (* $Id$ *)
 
@@ -27,12 +28,10 @@
 BeginPackage["KrancThorn`", {"CodeGen`", "Thorn`",
  "MapLookup`", "KrancGroups`", "Differencing`",
  "CalculationFunction`", "Errors`", "Helpers`", "CactusBoundary`",
- "TensorTools`", "Param`", "Schedule`", "Interface`", "Kranc`",
+ "KrancTensor`", "Param`", "Schedule`", "Interface`", "Kranc`", "Jacobian`",
  "ConservationCalculation`"}];
 
 CreateKrancThorn::usage = "Construct a Kranc thorn";
-CreateKrancThornTT::usage = "Construct a Kranc thorn using TensorTools";
-CreateGroupFromTensor::usage = "";
 
 Begin["`Private`"];
 
@@ -94,7 +93,7 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
     interface, evolvedGroupDefinitions, rhsGroupDefinitions, thornspec,
     allParams, boundarySources, reflectionSymmetries,
     realParamDefs, intParamDefs,
-    pDefs, useCSE, consCalcs, consCalcsIn, consGroups},
+    pDefs, consCalcs, consCalcsIn, consGroups},
 
     (* Parse named arguments *)
 
@@ -125,9 +124,11 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
     partialDerivs = OptionValue[PartialDerivatives] ~Join~
       ConservationDifferencingOperators[];
     reflectionSymmetries = OptionValue[ReflectionSymmetries];
-    useCSE = OptionValue[UseCSE];
 
     coordGroup = {"grid::coordinates", {Kranc`x,Kranc`y,Kranc`z,Kranc`r}};
+
+    CheckGroups[groupsOrig];
+
     groups = Join[groupsOrig, {coordGroup}];
     includeFiles = Join[includeFiles, {"GenericFD.h", "Symmetry.h", "sbp_calc_coeffs.h"}];
 
@@ -143,6 +144,8 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
     VerifyString[implementation];
     VerifyGroupNames[declaredGroups];
 
+    If[OptionValue[UseJacobian], JacobianCheckGroups[groups]];
+
     InfoMessage[Terse, "Creating startup file"];
     startup = CreateStartupFile[thornName, thornName];
 
@@ -196,7 +199,7 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
       inheritedRealParams, inheritedIntParams, inheritedKeywordParams,
       extendedRealParams, extendedIntParams, extendedKeywordParams,
       evolutionTimelevels, defaultEvolutionTimelevels,
-      calcs];
+      calcs, opts];
 
     (* Construct the schedule file *)
     InfoMessage[Terse, "Creating schedule file"];
@@ -221,7 +224,10 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
 
     (* Write the differencing header file *)
     InfoMessage[Terse, "Creating differencing header file"];
-    {pDefs, diffHeader} = CreateDifferencingHeader[partialDerivs, OptionValue[ZeroDimensions]];
+    {pDefs, diffHeader} = CreateDifferencingHeader[partialDerivs, OptionValue[ZeroDimensions], OptionValue[UseVectors]];
+    diffHeader = Join[
+        If[OptionValue[UseVectors], {"#include \"vectors.h\"\n", "\n"}, {}],
+        diffHeader];
 
     (* Add the predefinitions into the calcs *)
     calcs = Map[Join[#, {PreDefinitions -> pDefs}] &, calcs];
@@ -238,12 +244,12 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
     InfoMessage[Terse, "Creating calculation source files"];
     calcSources = Map[CreateSetterSource[
       {Join[#, {Parameters -> allParams, PartialDerivatives -> partialDerivs}]},
-      False, useCSE, {}, implementation, opts] &, calcs];
+      False, {}, implementation, opts] &, calcs];
     calcFilenames = Map[lookup[#, Name] <> ext &, calcs];
 
     (* Makefile *)
     InfoMessage[Terse, "Creating make file"];
-    make = CreateMakefile[Join[{"Startup.c", "RegisterMoL.c", "RegisterSymmetries.c"}, calcFilenames, 
+    make = CreateMakefile[Join[{"Startup.cc", "RegisterMoL.cc", "RegisterSymmetries.cc"}, calcFilenames, 
       Map[lookup[#, Filename] &, boundarySources]]];
 
     (* Put all the above together and generate the Cactus thorn *)
@@ -255,10 +261,10 @@ CreateKrancThorn[groupsOrig_, parentDirectory_, thornName_, opts:OptionsPattern[
                  Param         -> param,
                  Makefile      -> make,
                  Sources       -> Join[{
-                  {Filename -> "Startup.c", Contents -> startup}, 
-                  {Filename -> "RegisterMoL.c", Contents -> molregister},
-                  {Filename -> "RegisterSymmetries.c", Contents -> symregister},
-                  {Filename -> "Differencing.h", Contents -> diffHeader}},
+                  {Filename -> "Startup.cc", Contents -> startup}, 
+                  {Filename -> "RegisterMoL.cc", Contents -> molregister},
+                  {Filename -> "RegisterSymmetries.cc", Contents -> symregister},
+                 {Filename -> "Differencing.h", Contents -> diffHeader}},
                   MapThread[{Filename -> #1, Contents -> #2} &, 
                             {calcFilenames, calcSources}], boundarySources]};
     InfoMessage[Terse, "Creating thorn"];
@@ -327,6 +333,7 @@ createKrancMoLRegister[evolvedGroupNames_, nonevolvedGroupNames_, groups_, imple
     molregister = CreateMoLRegistrationSource[molspec, False];
     Return[molregister]];
 
+<<<<<<< HEAD
 (* --------------------------------------------------------------------------
    Tensors
    -------------------------------------------------------------------------- *)
@@ -457,5 +464,7 @@ CheckCalculationTensors[calc_] :=
   eqs = lookup[calc, Equations];
   Map[CheckEquationTensors, eqs]];
 
+=======
+>>>>>>> origin/master
 End[];
 EndPackage[];
diff --git a/Tools/CodeGen/Optimize.m b/Tools/CodeGen/Optimize.m
new file mode 100644
index 0000000..2f791da
--- /dev/null
+++ b/Tools/CodeGen/Optimize.m
@@ -0,0 +1,190 @@
+(* ::Package:: *)
+
+(*  Copyright 2011 Barry Wardell
+
+  This file is part of Kranc.
+
+  Kranc is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2 of the License, or
+  (at your option) any later version.
+
+  Kranc is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with Kranc; if not, write to the Free Software
+  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*)
+
+BeginPackage["Optimize`", {"Kranc`", "Errors`"}];
+
+EliminateCommonSubexpressions::usage = "EliminateCommonSubexpressions[calc] identifies common subexpressions in calc and introduces new shorthands for them.";
+TopologicallySortEquations::usage = "TopologicallySortEquations[eqs, v] sorts eqs topologically, including Head v as a possible vertex."
+
+Begin["`Private`"];
+
+CSEPrint[___] = null;
+(* CSEPrint = Print; *)
+
+Options[EliminateCommonSubexpressions] = ThornOptions;
+EliminateCommonSubexpressions[calc_List, OptionsPattern[]] :=
+ Module[{eqs, shorts, name, pdDefs, derivs, newShorts, newEqs, allShorts, newCalc},
+  name   = (Name /. calc);
+
+  InfoMessage[InfoFull, "Doing common subexpression elimination for "<>name];
+
+  eqs    = (Equations /. calc) /. Equations -> {};
+  shorts = (Shorthands /. calc) /. Shorthands -> {};
+
+  (* Get a list of symbols used for derivatives. We will not eliminate these as subexpressions. *)
+  pdDefs = OptionValue[PartialDerivatives];
+  derivs = DeleteDuplicates[Head/@(First/@pdDefs)];
+
+  (* Generate new equations with subexpressions eliminated. *)
+  {newShorts, newEqs} = cse[eqs, Symbol["csetemp"], derivs];
+
+  If[Length[newShorts]>0,
+    InfoMessage[Info, "Extracted "<>ToString[Length[newShorts]]<>" common subexpressions from "<>name];
+  ];
+
+  allShorts = Join[shorts, newShorts];
+  newCalc = Join[calc /. {(Shorthands->_) -> Sequence[], (Equations->_) -> Sequence[]},
+                 {Shorthands->allShorts}, {Equations->newEqs}];
+
+  newCalc
+];
+
+cse[eqs_, v_, exceptions_, minSaving_:0] :=
+ Module[{subexprs, replacements, replace, newEqs, defs, newDefs, i, relabelVars, allEqs, sortedEqs, newVars},
+  (* Find all possible subexpressions and how many times they occur *)
+  CSEPrint["CSE"];
+  CSEPrint["CSE: eqs=", eqs];
+  subexprs = Reap[Scan[If[! AtomQ[#], Sow[#]] &, eqs[[All,2]], Infinity]];
+  CSEPrint["CSE: subexprs=", subexprs];
+  If[subexprs[[2]]=={}, Return[{{}, eqs}]];
+  subexprs = Tally[subexprs[[2, 1]]];
+
+  (* Discard subexpressions which only appear once *)
+  subexprs = Select[subexprs, #[[2]] >= 2 &];
+
+  (* Sort subexpressions in ascending order by size (LeafCount) *)
+  subexprs = Sort[subexprs, LeafCount[#1] < LeafCount[#2] &];
+
+  (* Ony keep subexpressions larger than minSaving=(numoccurances-1)*size *)
+  subexprs = Select[subexprs, (#[[2]]-1) LeafCount[#[[1]]] >= minSaving &][[All,1]];
+
+  (* Discard some specific cases *)
+  subexprs = Cases[subexprs, Except[_?AtomQ]];
+  subexprs = Cases[subexprs, Except[Times[-1, _?AtomQ]]]; (* -x *)
+  subexprs = Cases[subexprs, Except[Alternatives@@(Blank/@exceptions)]]; (* specified exceptions *)
+  subexprs = Cases[subexprs, Except[Times[-1, Alternatives@@(Blank/@exceptions)]]]; (* -exceptions *)
+
+  (* Get the list of replacements for our original expression *)
+  replacements = Thread[subexprs -> Table[v[i], {i, Length[subexprs]}]];
+
+  (* Replace common subexpressions with new variables *)
+  (* Do not replace certain terms, e.g. the first argument of IfThen. *)
+  (* newEqs = eqs //. replacements; *)
+  CSEPrint["CSE: eqs=", eqs];
+  CSEPrint["CSE: replacements=", replacements];
+  replace[expr_] := Replace[Switch[expr,
+                                   IfThen[_,_,_], IfThen[expr[[1]], replace[expr[[2]]], replace[expr[[3]]]],
+                                   (* ToReal[_],     ToReal[expr[[1]]], *)
+                                   _?AtomQ,       expr,
+                                   _,             Map[replace, expr]],
+                            replacements];
+  newEqs = FixedPoint[replace, eqs];
+  CSEPrint["CSE: newEqs=", newEqs];
+
+  (* Build up definitions for the new variables *)
+  defs = Reverse/@replacements;
+  CSEPrint["CSE: defs=", defs];
+  For[i = 2, i <= Length[subexprs], i++,
+    defs[[i,2]] = defs[[i,2]] /. replacements[[1;;i-1]];
+  ];
+  CSEPrint["CSE: defs=", defs];
+
+  (* Select only the definitions which are needed for the new expressions.
+     This accounts for cases where a subexpression appears multiple times,
+     but always as part of the same larger subexpression. For example, in
+       expr = Sqrt[(a+b)(a-b)c]+(a+b)(a-b)c+(a+b)d+Sqrt[(a+b)d+(a+b)c];
+     we would identify the subexpressions
+       {v[1]->a+b,v[2]->d v[1],v[3]->a-b,v[4]->c v[1] v[3]};
+     whereas all we really want it to identify is
+       {v[1]->a+b,v[2]->d v[1],v[4]->(a-b) c v[1]};
+     and the introduction of v[3] is unnecessary. To achieve this, we only
+     keep temporary variables which appear in the expression after substitution
+     or which appear more than once in the definition of the temporary variables.
+  *)
+  newDefs = Select[defs, (Count[newEqs, #[[1]], Infinity] > 0) ||
+                         (Count[defs[[All,2]], #[[1]], Infinity] > 1) &];
+  CSEPrint["CSE: newDefs=", newDefs];
+
+  (* Replace any temporaries eliminated by the previous procedure with their definition *)
+  newDefs = newDefs //. Complement[defs, newDefs];
+  CSEPrint["CSE: newDefs2=", newDefs];
+
+  (* Check we actually have subexpressions to eliminate. Otherwise just return the original expression *)
+  If[Length[newDefs]==0, Return[{{}, eqs}]];
+
+  (* This is our new system of equations *)
+  allEqs = Join[newDefs, newEqs];
+  CSEPrint["CSE: allEqs=", allEqs];
+
+  sortedEqs = Fold[InsertNewEquation, newEqs, Reverse[newDefs]];
+  CSEPrint["CSE: sortedEqs=", sortedEqs];
+
+  (* Relabel new temporary variables so that they are sequential and C friendly *)
+  newVars = Select[sortedEqs[[All,1]], MemberQ[newDefs[[All, 1]],#]&];
+  CSEPrint["CSE: newVars=", newVars];
+  i = 0;
+  relabelVars = (# -> Symbol[ToString[v] <> ToString[i++]]) & /@ newVars;
+
+  (* Return the list of new variables and the new equations *)
+  {newVars, sortedEqs} /. relabelVars
+];
+
+TopologicallySortEquations[eqs_] := Module[{lhs, rhs, lhsInrhs, dag, sortedVars, indVars, allVars, sortedEqs},
+  lhs = eqs[[All,1]];
+  rhs = eqs[[All,2]];
+
+  (* Generate an directed acyclic graph for the system of equations *)
+  lhsInrhs = DeleteDuplicates[Cases[{#}, _?(MemberQ[lhs, #] &), Infinity]] & /@ rhs;
+  dag = Graph[ Flatten[MapThread[Thread[Rule[#1, #2]] &, {lhsInrhs, lhs}]] ];
+
+  (* Topologically sort the DAG *)
+  sortedVars = Quiet[TopologicalSort[dag], TopologicalSort::argx];
+
+  (* Check if the topological sorting failed. This can happen if the graph for
+     the equations is cyclic. For example, we could have a->a+1 or {b->a*a, a->b} *)
+  If[SameQ[Head[sortedVars], TopologicalSort],
+    InfoMessage[Info, "Failed to topologically sort equations."];
+    Return[$Failed]
+  ];
+
+  (* Some variables might be independent. Add them back in. *)
+  indVars = Complement[lhs, sortedVars];
+  allVars = Join[indVars, sortedVars];
+
+  sortedEqs = Thread[allVars -> (allVars/.eqs)];
+  sortedEqs
+];
+
+InsertNewEquation[oldEqs_, newEq_] := Module[{before},
+  CSEPrint["InsertNewEquation oldEqs=", oldEqs, " newEq=", newEq];
+  (* For some reason, we can be asked to insert an equation that is
+     not actually needed. This should not be the case. However, handle
+     it gracefully for now. *)
+  (* before = Position[oldEqs[[All,2]], newEq[[1]]][[1,1]]; *)
+  before = Position[oldEqs[[All,2]], newEq[[1]]];
+  If[before=={},
+     oldEqs,
+     Insert[oldEqs, newEq, before[[1,1]]]]
+];
+
+End[];
+
+EndPackage[];
diff --git a/Tools/CodeGen/Param.m b/Tools/CodeGen/Param.m
index 7821f40..1b4ad37 100644
--- a/Tools/CodeGen/Param.m
+++ b/Tools/CodeGen/Param.m
@@ -18,7 +18,7 @@
     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *)
 
-BeginPackage["Param`", {"Thorn`", "Errors`", "Helpers`", "MapLookup`", "KrancGroups`", "Kranc`"}];
+BeginPackage["Param`", {"Thorn`", "Errors`", "Helpers`", "MapLookup`", "KrancGroups`", "Kranc`", "Jacobian`"}];
 
 CreateKrancParam;
 MakeFullParamDefs;
@@ -72,14 +72,16 @@ krancParamStructExtended[definition_, type_] :=
      AllowedValues -> Map[{Value -> #, Description -> ""} &, allowedValues]}];
 
 krancKeywordParamStruct[struct_] :=
-{
-  Name -> lookup[struct, Name],
-  Type -> "KEYWORD",
-  Default -> lookup[struct, Default],
-  Description -> lookupDefault[struct, Description, lookup[struct, Name]],
-  Visibility -> lookupDefault[struct, Visibility, "private"],
-  AllowedValues -> Map[{Value -> #, Description -> #} &, lookup[struct, AllowedValues]]
-};
+  Join[
+  {Name -> lookup[struct, Name],
+   Type -> "KEYWORD",
+   Default -> lookup[struct, Default],
+   Description -> lookupDefault[struct, Description, lookup[struct, Name]],
+   Visibility -> lookupDefault[struct, Visibility, "private"]},
+  If[mapContains[struct, Steerable],
+    {Steerable -> lookup[struct, Steerable]},
+    {}],
+  {AllowedValues -> Map[{Value -> #, Description -> #} &, lookup[struct, AllowedValues]]}];
 
 MakeFullParamDefs[params_] :=
   Module[{p},
@@ -116,12 +118,13 @@ extendParameters[imp_, reals_, ints_, keywords_] :=
       Return[{Name -> imp, ExtendedParameters -> Join[realStructs, intStructs, keywordStructs]}],
       Return[{}]]];
 
+Options[CreateKrancParam] = ThornOptions;
 CreateKrancParam[evolvedGroups_, nonevolvedGroups_, groups_, thornName_, 
   reals_, ints_, keywords_,
   inheritedReals_, inheritedInts_, inheritedKeywords_,
   extendedReals_, extendedInts_, extendedKeywords_,
   evolutionTimelevels_, defaultEvolutionTimelevels_,
-  calcs_] :=
+  calcs_, opts:OptionsPattern[]] :=
   Module[{nEvolved, evolvedMoLParam, evolvedGFs,
     (*constrainedMoLParam,*) genericfdStruct, realStructs, intStructs,
     allInherited, allExtended, implementationNames, molImplementation,
@@ -146,7 +149,8 @@ CreateKrancParam[evolvedGroups_, nonevolvedGroups_, groups_, thornName_,
       Visibility -> "restricted",
       AccumulatorBase -> "MethodofLines::MoL_Num_Evolved_Vars",
       AllowedValues -> {{Value -> ToString[nEvolved] <> ":" <> ToString[nEvolved] , 
-                         Description -> "Number of evolved variables used by this thorn"}}
+                         Description -> "Number of evolved variables used by this thorn"}},
+      Steerable -> Recover
     };
 
     (*
@@ -171,7 +175,8 @@ CreateKrancParam[evolvedGroups_, nonevolvedGroups_, groups_, thornName_,
       Description -> "Number of active timelevels",
       Visibility -> "restricted",
       AllowedValues -> {{Value -> ToString[0] <> ":" <> ToString[evolutionTimelevels],
-                         Description -> ""}}
+                         Description -> ""}},
+      Steerable -> Recover
     };
 
     rhsTimelevelsParam =
@@ -182,25 +187,22 @@ CreateKrancParam[evolvedGroups_, nonevolvedGroups_, groups_, thornName_,
       Description -> "Number of active RHS timelevels",
       Visibility -> "restricted",
       AllowedValues -> {{Value -> ToString[0] <> ":" <> ToString[evolutionTimelevels],
-                         Description -> ""}}
+                         Description -> ""}},
+      Steerable -> Recover
     };
 
     genericfdStruct =
     {
       Name -> "GenericFD",
       UsedParameters -> 
-        {{Name -> "stencil_width",    Type -> "CCTK_INT"},
-         {Name -> "stencil_width_x",  Type -> "CCTK_INT"},
-         {Name -> "stencil_width_y",  Type -> "CCTK_INT"},
-         {Name -> "stencil_width_z",  Type -> "CCTK_INT"},
-         {Name -> "boundary_width",   Type -> "CCTK_INT"}}
+        If[OptionValue[UseJacobian], JacobianGenericFDParameters[], {}]
     };
 
     realStructs = Map[krancParamStruct[#, "CCTK_REAL", False] &, reals];
-    verboseStruct = krancParamStruct[{Name -> "verbose", Default -> 0}, "CCTK_INT", False];
+    verboseStruct = krancParamStruct[{Name -> "verbose", Default -> 0, Steerable -> Always}, "CCTK_INT", False];
     intStructs = Map[krancParamStruct[#, "CCTK_INT", False] &, ints];
-    calcEveryStructs = Map[krancParamStruct[{Name -> lookup[#, Name] <> "_calc_every", Default -> 1}, "CCTK_INT", False] &, calcs];
-    calcOffsetStructs = Map[krancParamStruct[{Name -> lookup[#, Name] <> "_calc_offset", Default -> 0}, "CCTK_INT", False] &, calcs];
+    calcEveryStructs = Map[krancParamStruct[{Name -> lookup[#, Name] <> "_calc_every", Default -> 1, Steerable -> Always}, "CCTK_INT", False] &, calcs];
+    calcOffsetStructs = Map[krancParamStruct[{Name -> lookup[#, Name] <> "_calc_offset", Default -> 0, Steerable -> Always}, "CCTK_INT", False] &, calcs];
     keywordStructs = Map[krancKeywordParamStruct, keywords];
 
     allInherited = Join[inheritedReals, inheritedInts, inheritedKeywords];
diff --git a/Tools/CodeGen/TensorTools.m b/Tools/CodeGen/TensorTools.m
index c6f4cfd..6579c63 100644
--- a/Tools/CodeGen/TensorTools.m
+++ b/Tools/CodeGen/TensorTools.m
@@ -217,10 +217,10 @@ IndexIsLower[TensorIndex[_, "u"]] := False;
    -------------------------------------------------------------------------- *)
 
 Format[TensorIndex[label_, "u"], OutputForm] :=
-  Superscript[null,label];
+  "u"<>ToString[label];
 
 Format[TensorIndex[label_, "l"], OutputForm] :=
-  Subscript[null,label];
+  "l"<>ToString[label];
 
 Format[TensorIndex[label_, "u"], StandardForm] :=
   Superscript[null,label];
@@ -251,12 +251,12 @@ DefineTensor[T_] :=
 
     Format[Tensor[T, is:((TensorIndex[_,_] | _Integer) ..) ], StandardForm] := 
       PrecedenceForm[
-        SequenceForm[T,is],
+        SequenceForm[T,"[",Sequence@@Riffle[{is},","],"]"],
         10000];
 
     Format[Tensor[T, is:((TensorIndex[_,_] | _Integer) ..) ], OutputForm] := 
       PrecedenceForm[
-        SequenceForm[T,is],
+        SequenceForm[T,"[",Sequence@@Riffle[{is},","],"]"],
         10000];
 
 (* Cannot get InputForm to work *)
@@ -265,7 +265,7 @@ DefineTensor[T_] :=
       HoldForm[T[is]];*)
 
     T[is:((TensorIndex[_,_] | _Integer) ..)] := Tensor[T, is];
-    TensorAttributes[T] = {TensorWeight -> 1, Symmetries -> {}};
+    TensorAttributes[T] = {TensorWeight -> 0, Symmetries -> {}};
     T];
 
 (* -------------------------------------------------------------------------- 
@@ -1235,10 +1235,11 @@ CheckTensors[x_, y_] :=
       ];
 
 CheckTensors[t:Tensor[k_, is__]] :=
-  Module[{},
+  Module[{is2},
 (*    Print["CheckTensors: Tensor: ", t];*)
-    If[!(Union[{is}] === Sort[{is}]),
-       ThrowError["Tensor has repeated indices: ", t, {is}]];
+    is2 = Select[{is}, !NumericQ[#]&];
+    If[!(Union[is2] === Sort[is2]),
+       ThrowError["Tensor has repeated indices: ", t, is2]];
     True];
 
 CheckTensors[t:f_[TensorIndex[__]..]] :=
diff --git a/Tools/CodeGen/TensorToolsKranc.m b/Tools/CodeGen/TensorToolsKranc.m
new file mode 100644
index 0000000..1abda8b
--- /dev/null
+++ b/Tools/CodeGen/TensorToolsKranc.m
@@ -0,0 +1,106 @@
+(* ::Package:: *)
+
+(*  Copyright 2004-2010
+    Sascha Husa, Ian Hinder, Christiane Lechner, Barry Wardell
+
+    This file is part of Kranc.
+
+    Kranc is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    Kranc is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Kranc; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+*)
+
+(* This package provides a TensorTools wrapper for Kranc *)
+
+BeginPackage["TensorToolsKranc`", {"Kranc`", "KrancGroups`", "TensorTools`"}];
+
+CreateGroupFromTensor::usage = "";
+ReflectionSymmetries::usage = "Produce a list of reflection symmetries of a tensor.";
+ExpandComponents::usage = "ExpandComponents[expr] converts an expression containing abstract indices into one containing components instead."
+IncludeCharacter::usage = "IncludeCharacter is an option for makeExplicit which specifies whether the character should also be included in the generated variable names."
+TensorCharacterString::usage = "TensorCharacterString[tensor[inds]] returns a string consisting of a sequence of U's and D's representing the character of tensor."
+
+Begin["`Private`"];
+
+ReflectionSymmetries[x___] := ReflectionSymmetriesOfTensor[x];
+ExpandComponents[expr_] := MakeExplicit[expr];
+
+TensorCharacterString[k_, inds_] :=
+  Module[{},
+  InfoMessage[InfoFull, "Tensor attributes of " <> ToString[k], TensorAttributes[k]];
+  If[HasTensorAttribute[k, TensorManualCartesianParities],
+    "ManualCartesian",
+    If[Length[inds] == 0, 
+      "Scalar",
+      Apply[StringJoin, Map[If[IndexIsUpper[#], "U", "D"] &, inds]]]]];
+
+CreateGroupFromTensor[T:Tensor[k_, is__]] :=
+  CreateGroupFromTensor[k, {is}];
+
+reflectionParityString[l_] :=
+  Module[{chars},
+    If[!ListQ[l] || !Length[l] == 3,
+      ThrowError["Expecting a list of three parities for TensorManualCartesianParities, must be 1 or -1"]];
+      
+      chars= Map[Switch[#, -1, "-", +1, "+", _, 
+                        ThrowError["Expecting a list of three parities for TensorManualCartesianParities, must be 1 or -1"]] &,
+                 l];
+
+      Apply[StringJoin, chars]];
+
+CreateGroupFromTensor[k_, inds_] :=
+  Module[{ttypeString, nInds, tags, group, vars},
+    InfoMessage[InfoFull, "Creating group from tensor with kernel " <> ToString[k] <> " and indices " <> ToString[inds]];
+    ttypeString = TensorCharacterString[k, inds];
+    InfoMessage[InfoFull, "Tensor type string: ", ttypeString];
+    nInds = Length[inds];
+    If[nInds == 2 && GetTensorAttribute[k, Symmetries] == {{2,1},1}, 
+        ttypeString = ttypeString <> "_sym"];
+    If[nInds == 3 && GetTensorAttribute[k, Symmetries] == {{1,3,2},1}, 
+        ttypeString = ttypeString <> "_sym"];
+    tags = {"tensortypealias" -> ttypeString, "tensorweight" -> GetTensorAttribute[k, TensorWeight]};
+    If[HasTensorAttribute[k, TensorSpecial],
+      tags = Append[tags, "tensorspecial" -> GetTensorAttribute[k, TensorSpecial]]];
+    If[HasTensorAttribute[k, TensorManualCartesianParities],
+      tags = Append[tags, "cartesianreflectionparities" -> 
+                          reflectionParityString[GetTensorAttribute[k, TensorManualCartesianParities]]]];
+    If[HasTensorAttribute[k, TensorParity],
+      tags = Append[tags, "tensorparity" -> GetTensorAttribute[k, TensorParity]]];
+
+    If[HasTensorAttribute[k, Checkpoint],
+      tags = Append[tags, "checkpoint" -> GetTensorAttribute[k, Checkpoint]]];
+
+    vars = If[nInds == 0, {k}, {Apply[Tensor, {k, Apply[Sequence,inds]}]}];
+    group = CreateGroup[ToString[k] <> "_group", vars, {Tags -> tags}];
+    Return[group]];
+
+CreateGroupFromTensor[x_] :=
+  If[IsTensor[x],
+    CreateGroupFromTensor[x, {}],
+    ThrowError["CreateGroupFromTensor: Not a tensor", x]];
+
+CheckEquationTensors[eq_] :=
+  Module[{},
+    CheckTensors[eq]];
+
+CheckCalculationTensors[calc_] :=
+  Module[{eqs},
+
+  If[mapContains[calc, Shorthands],
+  CheckTensors[lookup[calc, Shorthands]]];
+
+  eqs = lookup[calc, Equations];
+  Map[CheckEquationTensors, eqs]];
+
+End[];
+EndPackage[];
diff --git a/Tools/CodeGen/Thorn.m b/Tools/CodeGen/Thorn.m
index 06b3e8d..9381948 100644
--- a/Tools/CodeGen/Thorn.m
+++ b/Tools/CodeGen/Thorn.m
@@ -220,7 +220,8 @@ Options[CreateConfiguration] = ThornOptions;
 CreateConfiguration[opts:OptionsPattern[]] :=
   {whoWhen["CCL"],
    "REQUIRES GenericFD\n",
-   If[OptionValue[UseLoopControl], "REQUIRES LoopControl\n", {}]
+   If[OptionValue[UseLoopControl], "REQUIRES LoopControl\n", {}],
+   If[OptionValue[UseVectors], "REQUIRES Vectors\n", {}]
   };
 
 (* ------------------------------------------------------------------------ 
@@ -473,14 +474,19 @@ CreateSchedule[globalStorageGroups_, scheduledGroups_, scheduledFunctions_] :=
          optional LoopPreIncludes     -> {include file list},
 	          Equations           -> {{K11_rhs -> 2 A K11, ...}...}} *)
 
-calculationMacros[] :=
+calculationMacros[vectorise_] :=
   CommentedBlock["Define macros used in calculations",
-    Map[{"#define ", #, "\n"} &,
-       {"INITVALUE  (42)",
-        "INV(x) ((1.0) / (x))"   ,        
-        "SQR(x) ((x) * (x))"   ,        
-        "CUB(x) ((x) * (x) * (x))"   , 
-        "QAD(x) ((x) * (x) * (x) * (x))"}]];
+      Map[{"#define ", #, "\n"} &,
+         {"INITVALUE (42)",
+          "QAD(x) (SQR(SQR(x)))"} ~Join~
+          If[vectorise,
+           {"INV(x) (kdiv(ToReal(1.0),x))",
+            "SQR(x) (kmul(x,x))",
+            "CUB(x) (kmul(x,SQR(x)))"},
+           {"INV(x) ((1.0) / (x))",
+            "SQR(x) ((x) * (x))",
+            "CUB(x) ((x) * (x) * (x))"}]
+         ]];
 
 (* Given a list of Calculation structures as defined above, create a
    CodeGen representation of a source file that defines a function for
@@ -488,7 +494,7 @@ calculationMacros[] :=
 
 Options[CreateSetterSource] = ThornOptions;
 
-CreateSetterSource[calcs_, debug_, useCSE_, include_, imp_,
+CreateSetterSource[calcs_, debug_, include_, imp_,
   opts:OptionsPattern[]] :=
   Module[{},
 
@@ -503,21 +509,24 @@ CreateSetterSource[calcs_, debug_, useCSE_, include_, imp_,
          {IncludeSystemFile["assert.h"],
           IncludeSystemFile["math.h"],
           IncludeSystemFile["stdio.h"],
-          IncludeSystemFile["stdlib.h"]},
+          IncludeSystemFile["stdlib.h"],
+          IncludeSystemFile["string.h"]},
          {"\n"}
       ],
 
    Map[IncludeFile, Join[{"cctk.h", "cctk_Arguments.h", "cctk_Parameters.h",
-                     (*"precomputations.h",*) "GenericFD.h", "Differencing.h"}, include,
-                         If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}]]],
-   calculationMacros[],
+                         (*"precomputations.h",*) "GenericFD.h", "Differencing.h"},
+                         include,
+                         If[OptionValue[UseLoopControl], {"loopcontrol.h"}, {}],
+                         If[OptionValue[UseVectors], {"vectors.h"}, {}]]],
+   calculationMacros[OptionValue[UseVectors]],
 
    (* For each function structure passed, create the function and
       insert it *)
 
    CalculationBoundariesFunction[First[calcs], imp],
 
-   Map[CreateCalculationFunction[# , debug, useCSE, opts] &,
+   Map[CreateCalculationFunction[# , debug, imp, opts] &,
        calcs]}];
 
 
@@ -749,10 +758,10 @@ CreateMoLBoundariesSource[spec_] :=
      "if (CCTK_EQUALS(" <> boundpar <> ", \"none\"  ) ||\n",
      "    CCTK_EQUALS(" <> boundpar <> ", \"static\") ||\n",
      "    CCTK_EQUALS(" <> boundpar <> ", \"flat\"  ) ||\n",
-     "    CCTK_EQUALS(" <> boundpar <> ", \"zero\"  ) ) \n",
+     "    CCTK_EQUALS(" <> boundpar <> ", \"zero\"  ) )\n",
      "{\n",
 
-     "  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, \n",
+     "  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1,\n",
      "                    \"" <> fullgroupname <> "\", " <> boundpar <> ");\n",
 
      "  if (ierr < 0)\n",
@@ -771,10 +780,10 @@ CreateMoLBoundariesSource[spec_] :=
      "if (CCTK_EQUALS(" <> boundpar <> ", \"none\"  ) ||\n",
      "    CCTK_EQUALS(" <> boundpar <> ", \"static\") ||\n",
      "    CCTK_EQUALS(" <> boundpar <> ", \"flat\"  ) ||\n",
-     "    CCTK_EQUALS(" <> boundpar <> ", \"zero\"  ) ) \n",
+     "    CCTK_EQUALS(" <> boundpar <> ", \"zero\"  ) )\n",
      "{\n",
 
-     "  ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1, \n",
+     "  ierr = Boundary_SelectVarForBC(cctkGH, CCTK_ALL_FACES, 1, -1,\n",
      "                    \"" <> fullgfname <> "\", " <> boundpar <> ");\n",
 
      "  if (ierr < 0)\n",
@@ -1178,7 +1187,7 @@ headerComment2,
 "  assert (len);\n\n",
 
 "  for (d = 0; d < 3; ++d) {\n",
-"    assert (off[d] >= 0 && len[d] >= 0 && off[d] + len[d] <= cctk_lssh[CCTK_LSSH_IDX(0,d)]);\n",
+"    assert (off[d] >= 0 && len[d] >= 0 && off[d] + len[d] <= CCTK_LSSH(0,d));\n",
 "  }\n\n",
 
 "  assert (modes);\n",
@@ -1340,7 +1349,7 @@ CreateStartupFile[thornName_, bannerText_] :=
   tmp = {whoWhen["C"],
 
    IncludeFile["cctk.h"],
-   DefineFunction[thornName <> "_Startup", "int", "void",
+   DefineFunction[thornName <> "_Startup", "extern \"C\" int", "void",
      {DefineVariable["banner", "const char *", Quote[bannerText]],
       "CCTK_RegisterBanner(banner);\n",
       "return 0;\n"}]};
@@ -1354,7 +1363,7 @@ CreateStartupFile[thornName_, bannerText_] :=
    Thorn creation
    ------------------------------------------------------------------------ *)
 
-(* source = {Filename -> "MoLRegister.c", Contents -> "#include ..."} *)
+(* source = {Filename -> "MoLRegister.cc", Contents -> "#include ..."} *)
 
 (* thorn = {Name -> "ClassicADMMolEvolve", Directory -> "ClassicADM",
             Interface -> i, Schedule -> s, Param -> p, Makefile -> m, 
diff --git a/Tools/CodeGen/xTensorKranc.m b/Tools/CodeGen/xTensorKranc.m
new file mode 100644
index 0000000..f92072f
--- /dev/null
+++ b/Tools/CodeGen/xTensorKranc.m
@@ -0,0 +1,136 @@
+(* ::Package:: *)
+
+(*  Copyright 2010 Barry Wardell
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
+*)
+
+BeginPackage["xTensorKranc`", {"Differencing`", "Kranc`", "KrancGroups`", "xAct`xTensor`", "xAct`xCore`", "xAct`xCoba`"}];
+
+CreateGroupFromTensor::usage = "";
+ReflectionSymmetries::usage = "Produce a list of reflection symmetries of a tensor.";
+ExpandComponents::usage = "ExpandComponents[expr] converts an expression x containing abstract indices into one containing components instead."
+IncludeCharacter::usage = "IncludeCharacter is an option for makeExplicit which specifies whether the character should also be included in the generated variable names."
+TensorCharacterString::usage = "TensorCharacterString[tensor[inds]] returns a string consisting of a sequence of U's and D's representing the character of tensor."
+
+Begin["`Private`"];
+
+(* FIXME: Add support for ManualCartesian attribute *)
+TensorCharacterString[t_Symbol?xTensorQ[]] := "Scalar";
+TensorCharacterString[t_Symbol?xTensorQ[inds___]] := StringJoin[If[UpIndexQ[#],"U","D"]&/@{inds}];
+
+Options[makeExplicit] = {IncludeCharacter -> False};
+SetAttributes[makeExplicit, Listable];
+e : makeExplicit[_Plus, opts:OptionsPattern[]] := Distribute[Unevaluated[e]];
+makeExplicit[x_Times, opts:OptionsPattern[]] := Map[makeExplicit[#, opts]&, x];
+makeExplicit[Power[x_,p_], opts:OptionsPattern[]] := Power[makeExplicit[x, opts],p];
+makeExplicit[x_?NumericQ, opts:OptionsPattern[]] := x;
+makeExplicit[x_, {}, opts:OptionsPattern[]] := makeExplicit[x];
+makeExplicit[t_Symbol?xTensorQ[inds___], opts:OptionsPattern[]] := Module[{indexNumbers,character,indexString},
+  indexNumbers=First/@{inds};
+  
+  If[OptionValue[IncludeCharacter],
+    character = TensorCharacterString[t[inds]];
+    indexString = StringJoin[character,ToString/@indexNumbers],
+    indexString = StringJoin[ToString/@indexNumbers]
+  ];
+  SymbolJoin[PrintAs[t],Sequence@@indexString]
+];
+
+makeExplicit[x_, opts:OptionsPattern[]] := x;
+
+makeExplicit[(cd_?CovDQ)[ind_][expr_], opts:OptionsPattern[]] := Module[{indexNumbers},
+  indexNumbers=First/@{ind};
+
+  Global`PDstandard2nd[makeExplicit[expr, opts], Sequence@@indexNumbers]
+];
+
+Options[ExpandComponents] = Options[ExpandComponents];
+
+ExpandComponents[x_Rule, opts:OptionsPattern[makeExplicit]] := Thread[ExpandComponents[x[[1]], opts] -> ExpandComponents[x[[2]], opts]];
+ExpandComponents[dot[x_], opts:OptionsPattern[makeExplicit]] := dot/@ExpandComponents[x, opts];
+ExpandComponents[x_List, opts:OptionsPattern[makeExplicit]] := Flatten[Map[ExpandComponents[#, opts]&, x], 1];
+ExpandComponents[x_, opts:OptionsPattern[makeExplicit]] :=
+  Module[{eqs, options},
+
+  eqs = ComponentArray[TraceBasisDummy[x]];
+  options = Evaluate[FilterRules[{opts}, Options[makeExplicit]]];
+  If[Length[options]==0,
+    makeExplicit[eqs],
+    makeExplicit[eqs, options]
+  ]
+];
+
+
+(* Compute the reflection symmetries of a tensor *)
+ReflectionSymmetries[t_Symbol?xTensorQ[inds__]] :=
+  Module[{b=Global`Euclidean, cnums, components, componentIndices, counts},
+    (* Get the compoent indices of the basis *)
+    cnums = CNumbersOf[b, VBundleOfBasis[b]];
+
+    (* Get a list of components of the tensor t in the basis b *)
+    components = Flatten[ComponentArray[ToBasis[b][t[inds]]]];
+
+    (* Get the indices of each component *)
+    componentIndices = Map[IndicesOf[b], components];
+
+    (* Count the number of instances of each basis index. *)
+    countInds[expr_, basis_, cinds_] := Map[(Count[expr,{#,basis}]+Count[expr,{#,-basis}])&, cinds];
+    counts = Map[countInds[#, b, cnums]&, componentIndices];
+
+    (* For each instance, multiply by -1 *)
+    Thread[ExpandComponents[t[inds]] -> (-1)^counts]
+];
+
+ReflectionSymmetries[t_Symbol?xTensorQ[]] := t -> {1,1,1};
+ReflectionSymmetries[t_] := t -> {1, 1, 1};
+
+(* FIXME: Implement this fully *)
+GetTensorAttribute[t_Symbol?xTensorQ, TensorWeight] := WeightOfTensor[t];
+
+CreateGroupFromTensor[t_Symbol?xTensorQ[inds__]] := Module[{tCharString, nInds, tags, vars, group},
+  InfoMessage[InfoFull, "Creating group from tensor with kernel " <> SymbolName[t] <> " and indices " <> ToString[{inds}]];
+
+  (* Get a string representing the character of the tensor *)
+  tCharString = TensorCharacterString[t[inds]];
+  InfoMessage[InfoFull, "Tensor character string: ", tCharString];
+
+  (* Check if the tensor is symmetric *)
+  nInds = Length[SlotsOfTensor[t]];
+  If[SymmetryGroupOfTensor[t] == StrongGenSet[Range[nInds],GenSet[Cycles[Range[nInds]]]], 
+        tCharString = tCharString <> "_sym"];
+
+  (* FIXME: Add tensorspecial, cartesianreflectionparities  and tensorparity *)
+  tags = {"tensortypealias" -> tCharString, "tensorweight" -> GetTensorAttribute[t, TensorWeight]};
+
+  vars = If[nInds == 0, {t}, {t[inds]}];
+  group = CreateGroup[SymbolName[t] <> "_group", vars, {Tags -> tags}];
+  Return[group]
+];
+
+ReflectionSymmetries[x___]:= Throw["ReflectionSymmetries error: "<>ToString[x]];
+CreateGroupFromTensor[x___]:= Throw["CreateGroupFromTensor error: "<>ToString[x]];
+
+CheckTensors[expr_] := Validate[expr];
+
+End[];
+EndPackage[];
+
+
+
+
+
+
+