Subject: [FFC-dev] Consolidating dolfin-swig and dolfin

From: Johan Jansson <[email-masked]>
Date: Thu, 7 Sep 2006 03:23:07 +0200

FFC has a special output format for PyDOLFIN: "dolfin-swig", which is
very similar to the standard "dolfin" format, but generates a bit
simpler code since SWIG does not support some concepts in C++ (nested
classes primarily, and namespaces somewhat). This is bad for
maintainability, since a change in one format has to be made in the
other, possibly with some minor modification.

I've done some work in consolidating the two formats and I think I've
achieved an acceptable result. The new "dolfin" format is semantically
equivalent to the old format, so there is no visible change.

The two main changes are: moving nested class declarations out of the
enclosing class and generating a map (which is output as a SWIG file)
between DOLFIN class names and PyDOLFIN class names (the scope
resolution operator "::" doesn't exist in Python, so it has to be
removed or replaced).

The consequence is that SWIG can at least parse the nested classes,
though it (probably) ignores the fact that they are nested. This is a
hack (and an additional workaround is necessary to compensate for the
scoping), but it will have to do until there is better support in
SWIG. The SWIG developers have stated that nested classes won't be
implemented in the foreseeable future (i.e. until someone volunteers
responsibility of the problem), so we can't just wait.

I've attached the diff of the changes (against changeset
3a54d08c5669). Tell me if you have any objections, otherwise I'll go
ahead and make the corresponding changes to PyDOLFIN as well.

  Johan
diff -r 3a54d08c5669 src/ffc/format/dolfin.py
--- a/src/ffc/format/dolfin.py	Wed Jun 14 18:04:10 2006 +0200
+++ b/src/ffc/format/dolfin.py	Thu Sep 07 02:52:38 2006 +0200
@@ -37,6 +37,15 @@ format = { "sum": lambda l: " + ".join(l
            "tmp declaration": lambda j, k: "const real tmp%d_%d" % (j, k),
            "tmp access": lambda j, k: "tmp%d_%d" % (j, k) }
 
+def swig(swigmap):
+    output = ""
+
+    for p in swigmap.items():
+        output += "%%rename(%s) %s;\n" % (p[1], p[0])
+
+    return output
+
+
 def init(options):
     "Initialize code generation for DOLFIN format."
 
@@ -54,6 +63,9 @@ def write(forms, options):
 
     # Get name of form
     name = forms[0].name
+
+    # Initialize SWIG map
+    swigmap = {}
 
     # Write file header
     output = ""
@@ -80,8 +92,11 @@ be able to use it with DOLFIN.""", -1)
         else:
             xmlfile = "%s.xml" % forms[j].name
 
-        # Write form
-        output += __form(form, type, options, xmlfile)
+        # Write form prototype
+        output += __form(form, type, options, xmlfile, swigmap, True)
+
+        # Write form implementation
+        output += __form(form, type, options, xmlfile, swigmap, False)
 
     # Write file footer
     output += __file_footer()
@@ -92,6 +107,15 @@ be able to use it with DOLFIN.""", -1)
     file.write(output)
     file.close()
     debug("Output written to " + filename)
+
+    print swigmap
+
+    # Write SWIG file
+    swigfilename = name + ".i"
+    swigfile = open(swigfilename, "w")
+    swigfile.write(swig(swigmap))
+    swigfile.close()
+    debug("Swig output written to " + filename)
 
     # Write XML files if compiling for BLAS
     if options["blas"]:
@@ -196,26 +220,36 @@ def __file_footer_element():
 
 #endif\n"""
 
-def __elements(form):
+def __elements(form, subclass, swigmap, prototype = False):
     "Generate finite elements for DOLFIN."
 
     output = ""
     
     # Write test element (if any)
     if form.test:
-        output += __element(form.test, "TestElement")
+        output += __element(form.test, subclass,
+                            "TestElement", prototype)
+        swigmap["dolfin::" + form.name + "::" + subclass + "::" + "TestElement"] = \
+            form.name + subclass + "TestElement"
 
     # Write trial element (if any)
     if form.trial:
-        output += __element(form.trial, "TrialElement")
+        output += __element(form.trial, subclass,
+                            "TrialElement", prototype)
+        swigmap["dolfin::" + form.name + "::" + subclass + "::" + "TrialElement"] = \
+            form.name + subclass + "TrialElement"
 
     # Write function elements (if any)
     for j in range(len(form.elements)):
-        output += __element(form.elements[j], "FunctionElement_%d" % j)
+        output += __element(form.elements[j], subclass,
+                            "FunctionElement_%d" % j, prototype)
+        swigmap["dolfin::" + form.name + "::" + subclass + "::" + \
+                "FunctionElement_%d" % j] = \
+                form.name + subclass + "FunctionElement_%d" % j
     
     return output
 
-def __element(element, name):
+def __element(element, subclass, name, prototype = False):
     "Generate finite element for DOLFIN."
 
     # Generate code for initialization of tensor dimensions
@@ -267,7 +301,8 @@ def __element(element, name):
     subelements = ""
     if isinstance(element, MixedElement):
         for i in range(len(element.elements)):
-            subelements += __element(element.elements[i], "SubElement_%d" % i)
+            subelements += __element(element.elements[i], "",
+                                     "SubElement_%d" % i)
 
     # Generate code for FiniteElementSpec
     if element.type_str == "mixed":
@@ -280,9 +315,23 @@ def __element(element, name):
                (element.type_str, element.shape_str, element.degree())
 
     # Generate output
-    output = """\
-
-class %s : public dolfin::FiniteElement
+    if subclass != "":
+        subclass += "::"
+
+    if prototype == True:
+        # Write element prototype
+
+        output = """\
+  class %s;
+
+""" % (name)
+        return output
+    else:
+        # Write element implementation
+
+        output = """\
+
+class %s%s : public dolfin::FiniteElement
 {
 public:
 
@@ -360,7 +409,7 @@ private:
   FiniteElement** subelements;
 
 };
-""" % (name, name,
+""" % (subclass, name, name,
        diminit,
        elementinit,
        name,
@@ -377,9 +426,9 @@ private:
        spec,
        subelements)
 
-    return indent(output, 2)
-
-def __form(form, type, options, xmlfile):
+        return output
+
+def __form(form, type, options, xmlfile, swigmap, prototype = False):
     "Generate form for DOLFIN."
     
     #ptr = "".join(['*' for i in range(form.rank)])
@@ -395,8 +444,10 @@ def __form(form, type, options, xmlfile)
     if constinit:
         constinit = ", " + constinit
     
-    # Write class header
-    output = """\
+    if prototype == True:
+        # Write form prototype
+
+        output = """\
 /// This class contains the form to be evaluated, including
 /// contributions from the interior and boundary of the domain.
 
@@ -405,98 +456,125 @@ public:
 public:
 """ % (subclass, baseclass)
 
-    # Write elements
-    output += __elements(form)
-    
-    # Write constructor
-    output += """\
-
-  %s(%s) : dolfin::%s(%d)%s
-  {
-""" % (subclass, arguments, baseclass, form.nfunctions, constinit)
-
-    # Initialize test and trial elements
-    if form.test:
-        output += """\
-    // Create finite element for test space
-    _test = new TestElement();
-"""
-    if form.trial:
-        output += """\
-
-    // Create finite element for trial space
-    _trial = new TrialElement();
-"""
+        # Write element prototypes
+        output += __elements(form, subclass, swigmap, True)
+
+        # Write constructor
+        output += """\
+
+  %s(%s);
+""" % (subclass, arguments)
+
+        # Interior contribution
+        output += """\
+
+  void eval(real block[], const AffineMap& map) const;
+"""
+
+        # Boundary contribution (if any)
+        output += """\
+
+  void eval(real block[], const AffineMap& map, unsigned int facet) const;
+"""
+
+        # Declare class members (if any)
+        if form.nconstants > 0:
+            output += """\
+
+private:
+
+"""
+
+        # Create declaration list for for constants (if any)
+        if form.nconstants > 0:
+            for j in range(form.nconstants):
+                output += """\
+  const real& c%d;""" % j
+            output += "\n"
+
+        # Class footer
+        output += """
+};
+
+"""
+
+        # Write elements
+        output += __elements(form, subclass, swigmap, False)
+
+    else:
+        # Write form implementation
         
-    # Add functions (if any)
-    if form.nfunctions > 0:
-        output += """\
-
-    // Add functions\n"""
-        for j in range(form.nfunctions):
-            output += "    add(w%d, new FunctionElement_%d());\n" % (j, j)
-
-    # Initialize BLAS array (if any)
-    if options["blas"]:
-        output += """\
-
-    // Initialize form data for BLAS
-    blas.init(\"%s\");\n""" % xmlfile
-
-    output += "  }\n"
-
-    # Interior contribution (if any)
-    if form.AKi.terms:
-        eval = __eval_interior(form, options)
-        output += """\
-
-  void eval(real block[], const AffineMap& map) const
-  {
-%s  }
-""" % eval
-    else:
-        output += """\
-
-  // No contribution from the interior
-  void eval(real block[], const AffineMap& map) const {}
-"""
-
-    # Boundary contribution (if any)
-    if form.AKb[0].terms:
-        eval = __eval_boundary(form, options)
-        output += """\
-
-  void eval(real block[], const AffineMap& map, unsigned int facet) const
-  {
-%s  }
-""" % eval
-    else:
-        output += """\
-
-  // No contribution from the boundary
-  void eval(real block[], const AffineMap& map, unsigned int facet) const {}   
-"""
-
-    # Declare class members (if any)
-    if form.nconstants > 0:
-        output += """\
-
-private:
-
-"""
-
-    # Create declaration list for for constants (if any)
-    if form.nconstants > 0:
-        for j in range(form.nconstants):
-            output += """\
-  const real& c%d;""" % j
-        output += "\n"
-
-    # Class footer
-    output += """
-};
-
-"""
+        # Write constructor
+        output = """\
+
+%s::%s(%s) : dolfin::%s(%d)%s
+{
+""" % (subclass, subclass, arguments, baseclass, form.nfunctions, constinit)
+
+        # Initialize test and trial elements
+        if form.test:
+            output += """\
+  // Create finite element for test space
+  _test = new TestElement();
+"""
+        if form.trial:
+            output += """\
+
+  // Create finite element for trial space
+  _trial = new TrialElement();
+"""
+        
+        # Add functions (if any)
+        if form.nfunctions > 0:
+            output += """\
+
+  // Add functions\n"""
+            for j in range(form.nfunctions):
+                output += "  add(w%d, new FunctionElement_%d());\n" % (j, j)
+
+        # Initialize BLAS array (if any)
+        if options["blas"]:
+            output += """\
+
+  // Initialize form data for BLAS
+  blas.init(\"%s\");\n""" % xmlfile
+
+        output += "}\n"
+
+        # Interior contribution (if any)
+        if form.AKi.terms:
+            eval = __eval_interior(form, options)
+            output += """\
+
+void %s::eval(real block[], const AffineMap& map) const
+{
+%s}
+""" % (subclass, eval)
+        else:
+            output += """\
+
+// No contribution from the interior
+void %s::eval(real block[], const AffineMap& map) const {}
+""" % subclass
+
+        # Boundary contribution (if any)
+        if form.AKb[0].terms:
+            eval = __eval_boundary(form, options)
+            output += """\
+
+void %s::eval(real block[], const AffineMap& map, unsigned int facet) const
+{
+%s}
+""" % (subclass, eval)
+        else:
+            output += """\
+
+// No contribution from the boundary
+void %s::eval(real block[], const AffineMap& map, unsigned int facet) const {}   
+""" % subclass
+
+    swigmap["dolfin::" + form.name + "::" + subclass] = \
+        form.name + subclass
 
     return output
 
@@ -514,22 +592,22 @@ def __eval_interior_default(form, option
     if not options["debug-no-geometry-tensor"]:
         if len(form.cKi) > 0:
             output += """\
-    // Compute coefficients
-%s
-""" % "".join(["    const real %s = %s;\n" % (cKi.name, cKi.value) for cKi in form.cKi if cKi.used])
-        output += """\
-    // Compute geometry tensors
-%s"""  % "".join(["    const real %s = %s;\n" % (gK.name, gK.value) for gK in form.AKi.gK if gK.used])
-    else:
-        output += """\
-    // Compute geometry tensors
-%s""" % "".join(["    const real %s = 0.0;\n" % gK.name for gK in form.AKi.gK if gK.used])
+  // Compute coefficients
+%s
+""" % "".join(["  const real %s = %s;\n" % (cKi.name, cKi.value) for cKi in form.cKi if cKi.used])
+        output += """\
+  // Compute geometry tensors
+%s"""  % "".join(["  const real %s = %s;\n" % (gK.name, gK.value) for gK in form.AKi.gK if gK.used])
+    else:
+        output += """\
+  // Compute geometry tensors
+%s""" % "".join(["  const real %s = 0.0;\n" % gK.name for gK in form.AKi.gK if gK.used])
 
     if not options["debug-no-element-tensor"]:
         output += """\
 
-    // Compute element tensor
-%s""" % "".join(["    %s = %s;\n" % (aK.name, aK.value) for aK in form.AKi.aK])
+  // Compute element tensor
+%s""" % "".join(["  %s = %s;\n" % (aK.name, aK.value) for aK in form.AKi.aK])
 
     return output
 
@@ -541,24 +619,24 @@ def __eval_interior_blas(form, options):
     if not options["debug-no-geometry-tensor"]:
         if len(form.cKi) > 0:
             output += """\
-    // Compute coefficients
-%s
-""" % "".join(["    const real %s = %s;\n" % (cKi.name, cKi.value) for cKi in form.cKi if cKi.used])
-        output += """\
-    // Reset geometry tensors
-    for (unsigned int i = 0; i < blas.ni; i++)
-      blas.Gi[i] = 0.0;
-
-    // Compute entries of G multiplied by nonzero entries of A
-%s
-""" % "".join(["    blas.Gi[%d] = %s;\n" % (j, form.AKi.gK[j].value)
+  // Compute coefficients
+%s
+""" % "".join(["  const real %s = %s;\n" % (cKi.name, cKi.value) for cKi in form.cKi if cKi.used])
+        output += """\
+  // Reset geometry tensors
+  for (unsigned int i = 0; i < blas.ni; i++)
+    blas.Gi[i] = 0.0;
+
+  // Compute entries of G multiplied by nonzero entries of A
+%s
+""" % "".join(["  blas.Gi[%d] = %s;\n" % (j, form.AKi.gK[j].value)
                for j in range(len(form.AKi.gK)) if form.AKi.gK[j].used])
 
     # Compute element tensor
     if not options["debug-no-element-tensor"]:
         output += """\
-    // Compute element tensor using level 2 BLAS
-    cblas_dgemv(CblasRowMajor, CblasNoTrans, blas.mi, blas.ni, 1.0, blas.Ai, blas.ni, blas.Gi, 1, 0.0, block, 1);
+  // Compute element tensor using level 2 BLAS
+  cblas_dgemv(CblasRowMajor, CblasNoTrans, blas.mi, blas.ni, 1.0, blas.Ai, blas.ni, blas.Gi, 1, 0.0, block, 1);
 """
 
     return output
@@ -577,31 +655,31 @@ def __eval_boundary_default(form, option
     if not options["debug-no-geometry-tensor"]:
         if len(form.cKb) > 0:
             output += """\
-    // Compute coefficients
-%s
-""" % "".join(["    const real %s = %s;\n" % (cKb.name, cKb.value) for cKb in form.cKb if cKb.used])
-        output += """\
-    // Compute geometry tensors
-%s""" % "".join(["    const real %s = %s;\n" % (gK.name, gK.value) for gK in form.AKb[-1].gK if gK.used])
-    else:
-        output += """\
-    // Compute geometry tensors
-%s""" % "".join(["    const real %s = 0.0;\n" % gK.name for gK in form.AKb[-1].gK if gK.used])
+  // Compute coefficients
+%s
+""" % "".join(["  const real %s = %s;\n" % (cKb.name, cKb.value) for cKb in form.cKb if cKb.used])
+        output += """\
+  // Compute geometry tensors
+%s""" % "".join(["  const real %s = %s;\n" % (gK.name, gK.value) for gK in form.AKb[-1].gK if gK.used])
+    else:
+        output += """\
+  // Compute geometry tensors
+%s""" % "".join(["  const real %s = 0.0;\n" % gK.name for gK in form.AKb[-1].gK if gK.used])
 
     if not options["debug-no-element-tensor"]:
         output += """\
 
-    // Compute element tensor
-    switch ( facet )
-    { """
+  // Compute element tensor
+  switch ( facet )
+  { """
         for akb in form.AKb:
           output += """ 
-    case %s:"""  % akb.facet   
+  case %s:"""  % akb.facet   
           output += """ 
-%s      break; \n""" % "".join(["      %s = %s;\n" % (aK.name, aK.value) for aK in akb.aK])
-
-        output += """\
-    } \n"""
+%s      break; \n""" % "".join(["    %s = %s;\n" % (aK.name, aK.value) for aK in akb.aK])
+
+        output += """\
+  } \n"""
     return output
 
 def __eval_boundary_blas(form, options):
@@ -612,24 +690,24 @@ def __eval_boundary_blas(form, options):
     if not options["debug-no-geometry-tensor"]:
         if len(form.cKb) > 0:
             output += """\
-    // Compute coefficients
-%s
-""" % "".join(["    const real %s = %s;\n" % (cKb.name, cKb.value) for cKb in form.cKb if cKb.used])        
-        output += """\
-    // Reset geometry tensors
-    for (unsigned int i = 0; i < blas.nb; i++)
-      blas.Gb[i] = 0.0;
-
-    // Compute entries of G multiplied by nonzero entries of A
-%s
-""" % "".join(["    blas.Gb[%d] = %s;\n" % (j, form.AKb.gK[j].value)
+  // Compute coefficients
+%s
+""" % "".join(["  const real %s = %s;\n" % (cKb.name, cKb.value) for cKb in form.cKb if cKb.used])        
+        output += """\
+  // Reset geometry tensors
+  for (unsigned int i = 0; i < blas.nb; i++)
+    blas.Gb[i] = 0.0;
+
+  // Compute entries of G multiplied by nonzero entries of A
+%s
+""" % "".join(["  blas.Gb[%d] = %s;\n" % (j, form.AKb.gK[j].value)
                 for j in range(len(form.AKb.gK)) if form.AKb.gK[j].used])
 
     # Compute element tensor
     if not options["debug-no-element-tensor"]:
         output += """\
-    // Compute element tensor using level 2 BLAS
-    cblas_dgemv(CblasRowMajor, CblasNoTrans, blas.mb, blas.nb, 1.0, blas.Ab, blas.nb, blas.Gb, 1, 0.0, block, 1);
+  // Compute element tensor using level 2 BLAS
+  cblas_dgemv(CblasRowMajor, CblasNoTrans, blas.mb, blas.nb, 1.0, blas.Ab, blas.nb, blas.Gb, 1, 0.0, block, 1);
 """
 
     return output

---

From: Anders Logg <[email-masked]>
Date: Thu, 7 Sep 2006 09:28:06 +0200

I'll take a look.

There will be a fix to this problem in the not so distant future. I'm
working with Ola Skavhaug and Kent Andre Mardal on a new common
interface specification for the code generation. This will be named
UFC - Unified Form-assembly Code and be the default (and perhaps only)
output format from FFC. The idea is to have a common output format for
FFC and SyFi which is a similar tool developed by Kent.

There will also be a common input language: UFL - Unified Form
Language for both compilers.

The UFC will be chosen to work well with SWIG out of the box, so we
would not need any special solutions for SWIG.

I hope we'll have a first draft ready for the UFC within a month or
two.

/Anders

> _______________________________________________
> FFC-dev mailing list
> [email-masked]
> http://www.fenics.org/mailman/listinfo/ffc-dev

---

From: Johan Jansson <[email-masked]>
Date: Thu, 7 Sep 2006 12:50:04 +0200

> I hope we'll have a first draft ready for the UFC within a month or

Ok, great news.

  Johan

---

From: Anders Logg <[email-masked]>
Date: Fri, 8 Sep 2006 08:17:56 +0200

The new format breaks for some forms, in particular P5tet.form in the
demo directory of FFC:

logg@gwaihir:~/work/src/fenics/ffc/ffc-dev/src/demo$ ffc -d1 P5tet.py 
This is FFC, the FEniCS Form Compiler, version 0.3.3.
For further information, go to http://www/fenics.org/ffc/.
Parsing P5tet.py
Output written to P5tet.py.py
Compiling finite element: Lagrange finite element of degree 5 on a
tetrahedron

Generating output for DOLFIN
Traceback (most recent call last):
  File "/usr/bin/ffc", line 128, in ?
    main(sys.argv[1:])
  File "/usr/bin/ffc", line 82, in main
    execfile(outname, ns)
  File "P5tet.py.py", line 38, in ?
    writeFiniteElement(element, name, "dolfin", {'blas': False,
    'no-gpl': False, 'optimize': False, 'debug-no-geometry-tensor':
    False, 'debug-no-element-tensor': False})
  File "/usr/lib/python2.4/site-packages/ffc/compiler/compiler.py",
  line 191, in writeFiniteElement
    format.writeFiniteElement(element, name, options)
  File "/usr/lib/python2.4/site-packages/ffc/format/dolfin.py", line
  136, in writeFiniteElement
    output += __element(element, name)
TypeError: __element() takes at least 3 arguments (2 given)

Otherwise, the new format looks ok. If you can send a new version that
can be verified to work for all FFC demos and all DOLFIN demos, I'll
commit the new version right away.

/Anders

---

From: Johan Jansson <[email-masked]>
Date: Wed, 13 Sep 2006 16:48:30 +0200

> The new format breaks for some forms, in particular P5tet.form in the
...

Ok, I didn't consider compiling just the finite element. This is fixed
now.

> commit the new version right away.

I've tested all the FFC forms and all the DOLFIN demo forms. I've also
run some of the demos but not all (should be put in as unit tests in
the future).

I'm attaching the changesets as a bundle. Do "hg unbundle" to "pull"
them from the bundle. This seems to be a much better concept than
diffs.

  Johan

---

From: Anders Logg <[email-masked]>
Date: Wed, 13 Sep 2006 16:51:04 +0200

Very good. I'll run some more tests tomorrow and then commit the
update.

/Anders

---

From: Anders Logg <[email-masked]>
Date: Thu, 14 Sep 2006 16:35:01 +0200

Fixed. I've recompiled all the forms in DOLFIN (srcipts/compileforms)
and will run some more tests before I make the update to DOLFIN.

Can we remove the DOLFIN SWIG format?

/Anders

---

From: Anders Logg <[email-masked]>
Date: Thu, 14 Sep 2006 17:36:53 +0200

Seems to work but it breaks PyDOLFIN (multiply defined funtions during
linking). Do you have a fix waiting for PyDOLFIN?

/Anders

---

From: Johan Jansson <[email-masked]>
Date: Thu, 14 Sep 2006 19:20:42 +0200

> linking). Do you have a fix waiting for PyDOLFIN?

Yes, PyDOLFIN also needs some updates which are not checked in
yet. Perhaps you could wait until after that has been tested with
removing the DOLFIN-SWIG format.

  Johan

---

From: Anders Logg <[email-masked]>
Date: Thu, 14 Sep 2006 19:43:02 +0200

> removing the DOLFIN-SWIG format.

ok.

/Anders

---

From: Anders Logg <[email-masked]>
Date: Thu, 21 Sep 2006 10:19:48 +0200

> removing the DOLFIN-SWIG format.

Has this been fixed or are FFC and DOLFIN still out of sync?

/Anders

---

From: Johan Jansson <[email-masked]>
Date: Thu, 21 Sep 2006 12:38:55 +0200

> Has this been fixed or are FFC and DOLFIN still out of sync?

I've checked it in now, I've tested some of the demos and they work
ok.

  Johan

---

From: Johan Jansson <[email-masked]>
Date: Thu, 21 Sep 2006 12:42:06 +0200

On Thu, Sep 21, 2006 at 12:38:55PM +0200, Johan Jansson wrote:

...

> I've checked it in now, I've tested some of the demos and they work

This also means the old "dolfin-swig" format is obsolete and can be
removed.

  Johan

---

From: Anders Logg <[email-masked]>
Date: Thu, 21 Sep 2006 20:02:54 +0200

> I've checked it in now, I've tested some of the demos and they work

Great. I'll remove the swig format and recompile all forms in DOLFIN.

/Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Wed, 27 Sep 2006 14:24:12 +0200

> I hope we'll have a first draft ready for the UFC within a month or

In the light of this, is it worthwhile fixing the blas output option for
boundary integrals or should I wait?

Garth

---

From: Anders Logg <[email-masked]>
Date: Wed, 27 Sep 2006 19:25:10 +0200

> boundary integrals or should I wait?

It depends. It will be another couple of months until the new
framework is in place, so if you want it now, feel free to fix it.

Having something that works now in FFC will simplify reimplementation
but I will still need to (re-)think about all details, including the
BLAS support (which should also be extended to BLAS 3 operations).

/Anders