Unrolling array expressions

[Christian Schubert]

Hi,

Bruno has written a great library for simulating flexible bodies. Unfortunately, OMC is not able to handle it because of the extensive use of multi dimensional matrices in the code.

I created a small example to illustrate the problem:

model FemPackage_Test
  constant Integer N = 10;
  constant Integer M = 30;
  parameter Real A[N,N] = identity(N);
  parameter Real B[N,N,M] = ones(N,N,M);
  Real q[M](each start = 1, each fixed=true);
  Real C[N,N];
equation
  C = A + sum(B[:,:,i]*q[i] for i in 1:M);
  der(q) = -q;
end FemPackage_Test;

Since every expression is rolled out in the frontend, the expressions become very large and thus memory consumption as well as flattening time are enormous. Life sized models have many such expressions with much larger dimensions.

My question is how we can circumvent this? All the multidimensional matrices are known constants, so it is not required to unroll them for matching in the backend.
Would it be possible to not unroll expressions above a 'certain size'?

It would be great if that library could be used with OpenModelica!

PS: This model does not compile with OpenModelica:

expected 'const struct real_array_t *' but argument is of type 'modelica_real *'

[Martin Sjölund]

This ticket was not closed for 1.9.1, which has now been released. It was batch modified for milestone 1.9.2 (but maybe an empty milestone was more appropriate; feel free to change it).

[Martin Sjölund]

Milestone changed to 1.9.3 since 1.9.2 was released.

[Martin Sjölund]

Moved to new milestone 1.9.4

[Martin Sjölund]

Milestone pushed to 1.9.5

[Martin Sjölund]

Milestone renamed

[Martin Sjölund]

Ticket retargeted after milestone closed

[Martin Sjölund]

Milestone moved to 1.12.0 due to 1.11.0 already being released.

[Francesco Casella]

The milestone of this ticket has been reassigned to "Future".

If you think the issue is still valid and relevant for you, please select milestone 1.13.0 for back-end, code generation and run-time issues, or 2.0.0 for front-end issues.

If you are aware that the problem is no longer present, please select the milestone corresponding to the version of OMC you used to check that, and set the status to "worksforme".

In both cases, a short informative comment would be welcome.

[Francesco Casella]

I tried this example in Dymola, it compiles and runs in a few seconds. OMC fails with this message:

[1] 00:47:34 Translation Error
[FemPackage_Test: 9:3-9:42]: Variable B[1,:,i,1] not found in scope
FemPackage_Test.$foriter loop scope$.

It would be good if this worked nicely with the new FE.

[Martin Sjölund]

This model works fine in the old frontend with -v=1; it seems when expanding the matrix something goes wrong.

[Martin Sjölund]

​https://github.com/OpenModelica/OMCompiler/pull/2018 makes this work slightly better. Code generation fails because we cannot generate code for B[:,:,i] in simulation context. -v=1 works fine though.

[Francesco Casella]

I tried a smaller version (M = 4, N = 4) with the NF and the latest nightly build. I get this error message:

[1] 01:53:02 Translation Error 
[FemPackage_Test: 9:3-9:42]:
Cannot resolve type of expression A + sum(array(B[:, :, i] .* (q[i]) for i in 1:4)).
The operands have types Real[4, 4], Real in component <NO_COMPONENT>.

• it would be nice if the NF could handle this model
• it would also be nice if the multiplications by constant zeros and ones were carried out in order to avoid the explosion of memory allocation and code generation times when N and M become larger. This should probably be carried out in the front-end, but it could also be handled in the back-end as long as the front-end does not expand the expressions.

[Per Östlund]

The type matching issue is because reductions are implemented wrong in the NF, and it substitutes sum(... for ...) with sum(array(... for ...)). This isn't a valid substitution since sum of an array will sum all the scalars in the array, always giving a scalar result. Hence why the NF thinks the sum expression has type Real instead of Real[4, 4] like it should be. See also #4753.

[Francesco Casella]

Still valid as of v1.13.0-dev-798-g1c8bb86de

[Francesco Casella]

In general, the theme of avoiding expansion in the front-end and having the back-end and code generation handle arrays is a hot topic, and there will be developments in this direction in OMC in the next few years. Some early attempts are reported here: #5110, #5144.

Regarding this specific model, the new frontend (-d=newInst) can now handle this model a lot more efficiently. It is necessary to add the annotation Evaluate = true to all parameters (which I guess is something Dymola does automatically), so that the back-end can symbolically eliminate all the products by zero, and also to hide the A and B matrices, to avoid cluttering the result file with a lot of zeros.

model FemPackage_Test
  parameter Integer N = 10 annotation(Evaluate = true);
  parameter Integer M = 30 annotation(Evaluate = true);
  parameter Real A[N,N] = identity(N) annotation(Evaluate = true, HideResult = true);
  parameter Real B[N,N,M] = ones(N,N,M) annotation(Evaluate = true, HideResult = true);
  Real q[M](each start = 1, each fixed=true);
  Real C[N,N];
equation
  C = A + sum(B[:,:,i]*q[i] for i in 1:M);
  der(q) = -q;
annotation(__OpenModelica_commandLineOptions = "-d=newInst");
end FemPackage_Test;

Currently, the code generation process still goes through full expansion, so it won't scale well for really large models. With the current nightly build v1.14.0-dev-105-g2de3e7582, the above model is compiled in about 12 seconds and runs in 20 ms, which is not bad. If the -d=newInst,dumpindxdae compiler flags are used, the reported optimized DAEs after symbolic manipulation

1/1 (1): C[1,1] = 1.0 + q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
2/2 (1): C[1,2] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
3/3 (1): C[1,3] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
4/4 (1): C[1,4] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
5/5 (1): C[2,1] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
6/6 (1): C[2,2] = 1.0 + q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
7/7 (1): C[2,3] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
8/8 (1): C[2,4] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
9/9 (1): C[3,1] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
10/10 (1): C[3,2] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
11/11 (1): C[3,3] = 1.0 + q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
12/12 (1): C[3,4] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
13/13 (1): C[4,1] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
14/14 (1): C[4,2] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
15/15 (1): C[4,3] = q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
16/16 (1): C[4,4] = 1.0 + q[1] + q[2] + q[3] + q[4]   [dynamic |0|0|0|0|] 
17/17 (1): der(q[1]) = -q[1]   [dynamic |0|0|0|0|] 
18/18 (1): der(q[2]) = -q[2]   [dynamic |0|0|0|0|] 
19/19 (1): der(q[3]) = -q[3]   [dynamic |0|0|0|0|] 
20/20 (1): der(q[4]) = -q[4]   [dynamic |0|0|0|0|] 

show that the computation of C is as efficient as possible, since the B array has been optimized away.

I would say that what we now have in OMC is good for this size of the model, which is probably enough for most applications of flexible multibody systems, which account for a relatively low number of modes. Once we have a fully vectorized back-end, we could completely avoid expanding the matrices in the front end, and handle even larger models easily.