Skip to content
Snippets Groups Projects
Commit 05c64326 authored by Marcel Koch's avatar Marcel Koch
Browse files

use fma for inverse 

3d -> 9 fma + 18 mul (without -1 *)
2 fma more than C++, 4 mul/add more than C++
parent dd00c10b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
python/dune/codegen/pdelab/tensors.py
+ 21
− 31
View file @ 05c64326
......@@ -24,8 +24,7 @@ def define_determinant(name, matrix, shape, visitor):
matrix_entry = [[prim.Subscript(prim.Variable(matrix), (i, j)) for j in range(dim)] for i in range(dim)]
if dim == 2:
expr_determinant = FMA(matrix_entry[0][0], matrix_entry[1][1],
-1 * prim.Product((matrix_entry[1][0], matrix_entry[0][1])))
expr_determinant = FMA(matrix_entry[0][0], matrix_entry[1][1], -1 * matrix_entry[1][0] * matrix_entry[0][1])
elif dim == 3:
fma_A = FMA(matrix_entry[1][1], matrix_entry[2][2], -1 * matrix_entry[1][2] * matrix_entry[2][1])
......@@ -73,35 +72,26 @@ def define_matrix_inverse(name, name_inv, shape, visitor):
sign = 1. if i == j else -1.
exprs[i][j] = prim.Product((sign, prim.Variable(det_inv), matrix_entry[1 - i][1 - j]))
elif dim == 3:
exprs[0][0] = prim.Product((1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[1][1], matrix_entry[2][2])),
-1 * prim.Product((matrix_entry[1][2], matrix_entry[2][1]))))))
exprs[1][0] = prim.Product((-1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[0][1], matrix_entry[2][2])),
-1 * prim.Product((matrix_entry[0][2], matrix_entry[2][1]))))))
exprs[2][0] = prim.Product((1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[0][1], matrix_entry[1][2])),
-1 * prim.Product((matrix_entry[0][2], matrix_entry[1][1]))))))
exprs[0][1] = prim.Product((-1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[1][0], matrix_entry[2][2])),
-1 * prim.Product((matrix_entry[1][2], matrix_entry[2][0]))))))
exprs[1][1] = prim.Product((1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[2][2])),
-1 * prim.Product((matrix_entry[0][2], matrix_entry[2][0]))))))
exprs[2][1] = prim.Product((-1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][2])),
-1 * prim.Product((matrix_entry[0][2], matrix_entry[1][0]))))))
exprs[0][2] = prim.Product((1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[1][0], matrix_entry[2][1])),
-1 * prim.Product((matrix_entry[1][1], matrix_entry[2][0]))))))
exprs[1][2] = prim.Product((-1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[2][1])),
-1 * prim.Product((matrix_entry[0][1], matrix_entry[2][0]))))))
exprs[2][2] = prim.Product((1., prim.Variable(det_inv),
prim.Sum((prim.Product((matrix_entry[0][0], matrix_entry[1][1])),
-1 * prim.Product((matrix_entry[0][1], matrix_entry[1][0]))))))
exprs[0][0] = prim.Variable(det_inv) * FMA(matrix_entry[1][1], matrix_entry[2][2],
-1 * matrix_entry[1][2] * matrix_entry[2][1])
exprs[1][0] = prim.Variable(det_inv) * FMA(matrix_entry[0][1], matrix_entry[2][2],
-1 * matrix_entry[0][2] * matrix_entry[2][1]) * -1
exprs[2][0] = prim.Variable(det_inv) * FMA(matrix_entry[0][1], matrix_entry[1][2],
-1 * matrix_entry[0][2] * matrix_entry[1][1])
exprs[0][1] = prim.Variable(det_inv) * FMA(matrix_entry[1][0], matrix_entry[2][2],
-1 * matrix_entry[1][2] * matrix_entry[2][0]) * -1
exprs[1][1] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[2][2],
-1 * matrix_entry[0][2] * matrix_entry[2][0])
exprs[2][1] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[1][2],
-1 * matrix_entry[0][2] * matrix_entry[1][0]) * -1
exprs[0][2] = prim.Variable(det_inv) * FMA(matrix_entry[1][0], matrix_entry[2][1],
-1 * matrix_entry[1][1] * matrix_entry[2][0])
exprs[1][2] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[2][1],
-1 * matrix_entry[0][1] * matrix_entry[2][0]) * -1
exprs[2][2] = prim.Variable(det_inv) * FMA(matrix_entry[0][0], matrix_entry[1][1],
-1 * matrix_entry[0][1] * matrix_entry[1][0])
else:
raise NotImplementedError
for j in range(dim):
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Datenschutz und Impressum