[PSyAd] [LFRic] transpose intrinsic not correctly handled in sympy writer/reader.

[DrTVockerodtMO]

In MetOffice/lfric_apps#545 we found a failing adjoint test, where in our current software stack it properly features the line

jac_inv_transpose(:,:) = TRANSPOSE(jac_inv)

which is a direct copy from the forward. However, when running PSyAd using PSyclone 3.3 we found that the line is replaced with

jac_inv_transpose(:,:) = jac_inv

which is not expected. These matrices are passive variables so PSyAd ought to replicate them properly. This was already after the forward code was patched to include this line-by-line.

[arporter]

Thanks @DrTVockerodtMO, please see my comment on the original ticket.

[arporter]

Thanks to @sergisiso, we've discovered that if you roll back sympy to 1.13.3 this bug goes away. I'm going to be slightly cheeky and tag @hiker here as he is our resident sympy expert. Essentially Joerg, the Sympy writer/reader doesn't work properly with transpose which of course is a Sympy intrinsic itself.

[arporter]

(Note that, for the 3.3.0 release of PSyclone we will pin the version of sympy that is installed in order to workaround this issue in the short term.)

[hiker]

I could reproduce the problem:

/psyclone/tests$ git diff psyir/frontend/sympy_reader_test.py
diff --git a/src/psyclone/tests/psyir/frontend/sympy_reader_test.py b/src/psyclone/tests/psyir/frontend/sympy_reader_test.py
index 5259bb85fc..2961cfdb71 100644
--- a/src/psyclone/tests/psyir/frontend/sympy_reader_test.py
+++ b/src/psyclone/tests/psyir/frontend/sympy_reader_test.py
@@ -66,6 +66,7 @@ def test_sympy_reader_constructor():
                                          ("b(2:3:4)", "b(2:3:4)"),
                                          ("b(2:3:1)", "b(2:3)"),
                                          ("b", "b"),
+                                         ("transpose(b)", "TRANSPOSE(b)\n"),
                                          ("b(:)", "b(:)"),
                                          ("b(::)", "b(:)"),
                                          ("b(::1)", "b(:)"),

But the main issue is that sympy apparently 'resolves' the transmute call in 1.14 - it disappears and is not part of the tree anymore:

from sympy import Symbol, transpose
x_nonc=Symbol("X_nonc", commutative=False)
x_c=Symbol("X_c")
transpose(x_c)
X_c
transpose(x_nonc)
transpose(X_nonc)

Work arounds as shown to declare the symbol commutative (which causes many other failures), declare the symbol to be a matrix (which we might not know in general??).

There might be some options to either overwrite Sympy's transpose (to be not evaluated), or create a derived symbol class (which disables evaluation of transpose).

[arporter]

But the main issue is that sympy apparently 'resolves' the transmute call in 1.14

Aha, I wondered about that. I guess, since we're dealing with Fortran, if there's a TRANSPOSE then we know its argument must be a matrix. That might be complicated though?

[sergisiso]

But isn't it a good things that it resolve intrinsics? This makes is symbolic comparisons better.

The question is why this end up affecting the adjoint generated code?

[arporter]

I guess the problem is that Sympy thinks the argument to transpose is just a scalar symbol and therefore the transpose operation is redundant and can be optimised out. It doesn't know that it's actually a matrix.

[hiker]

I guess the problem is that Sympy thinks the argument to transpose is just a scalar symbol and therefore the transpose operation is redundant and can be optimised out. It doesn't know that it's actually a matrix.

Yes, indeed, at this stage the type information is not used, so all symbols are scalars. I am also not certain what would happen if we declare arrays as SymPy matrices, though now that we have better type information, that might be a great next step (I vaguely remember looking at matrices in SymPy, but its support was ... somehow lacking at that time and so I abandoned it. No idea what it was).

Note that it will be a few weeks before I can get back to that. One in-between solution would be to use create the SymPy Transpose node ourselves, but set it to not evaluate - that should also work, without introducing matrices.