Complex matrix vector multiplication added in MSL trunk #1890

modelica-trac-importer · 2017-01-15T07:50:35Z

Reported by beutlich on 25 Jan 2016 09:04 UTC
The newly introduced matrix vector product of type Complex (e.g., in Modelica.Magnetic.QuasiStatic.FundamentalWave.Components.MultiPhaseElectroMagneticConverter.vSymmetricalComponent) is currently not supported by all tools.

We therefore propose to add a new function matrixVectorProduct to Modelica.ComplexMath in order to circumvent this problem.

The following test model contains three implementations of the complex matrix vector product.

model TestComplexMatrixVectorMultiplication
  parameter Complex m[2,2] = {{Complex(1,0),Complex(1,1)},{Complex(1,-1),Complex(0,1)}};
  parameter Complex v[2] = {Complex(1,0),Complex(0,1)};
  Complex mv1[2];
  Complex mv2[2];
  Complex mv3[2];

  // Modelica.ComplexMath.matrixVectorProduct
  function matrixVectorProduct "Matrix vector product c1*c2 of a complex matrix and a complex vector"
    import Complex;
    input Complex c1[:,:] "Vector of Complex numbers 1";
    input Complex c2[size(c1,2)] "Vector of Complex numbers 2";
    output Complex c3[size(c1,1)] "= c3[j]:= sum_k(c1[j,k]*c2[k])";
  algorithm
    for j in 1:size(c1,1) loop
      c3[j] := Complex(0);
      for k in 1:size(c1,2) loop
        c3[j] := c3[j] + c1[j,k]*c2[k];
      end for;
    end for;
  annotation(Documentation(info="<html>This function returns the matrix vector product of a given matrix of Complex and a given vector of Complex.</p></html>"));
  end matrixVectorProduct;

equation
  mv1 = m*v; // New in MSL trunk
  mv2 = {sum(m[j,k]*v[k] for k in 1:size(m,2)) for j in 1:size(m,1)}; // Illegal Modelica since arguments of sum must be of type Real or Integer
  mv3 = matrixVector(m, v); // Explicit use function matrixVector
  
annotation(uses(Complex));
end TestComplexMatrixVectorMultiplication;

Migrated-From: https://trac.modelica.org/Modelica/ticket/1890

modelica-trac-importer · 2017-01-15T07:50:44Z

Modified by beutlich on 2 Feb 2016 12:48 UTC

modelica-trac-importer · 2017-01-15T07:51:04Z

Comment by ahaumer on 3 Feb 2016 14:10 UTC
ok I see no problem implementing the function.
The problem is: How shall I find ALL occurences of Complex matrix x Complex vector replace the multiplication by the function call? Could you prove a list?

modelica-trac-importer · 2017-01-15T07:51:14Z

Comment by ahaumer on 3 Feb 2016 15:14 UTC
ok implemented the functon; I found 3 (all?) occurences:
1 in Modelica.Electrical.QuasiStationary.MulitPhase.Blocks.SymmetricalComponents
2 in Modelica.Magnetic.QuasiStatic.FundamentalWave.Components.MultiPhaseElectroMagneticConverter
Added your test case to ModelicaTest; added a line to the release notes.
Resolved by e24066a.
Closing the ticket; please reopen if you stumble upon more occurences.

modelica-trac-importer · 2017-01-15T07:51:44Z

Comment by ahaumer on 3 Feb 2016 15:46 UTC
Replacing symmetricTransformationMatrix(m)*v by matrixVectorProduct(symmetricTransformationMatrix(m), i) in Modelica.Magnetic.Quasistatic.FundamentalWave.Components.MultiPhaseElectroMagneticConverter breaks examples, e.g. Modelica.Magnetic.Quasistatic.FundamentalWave.Examples.Components.MultiPhaseInductance.
Searching for a solution ...
btw there's a more elegant formulaton of the matrixVectorProduct (maybe that's the solution):

function matrixVectorProduct 
  "Returns the product of a complex matrix and a complex vector"
  extends Modelica.Icons.Function;
  input Complex m[:,:] "Complex matrix";
  input Complex v[size(m,2)] "Complex vector";
  output Complex result[size(m,1)] "Complex result vector m*v";
algorithm 
  result:={Modelica.ComplexMath.'sum'({m[j,k]*v[k] for k in 1:size(m,2)}) for j in 1:size(m,1)};
end matrixVectorProduct;

modelica-trac-importer · 2017-01-15T07:52:14Z

Comment by ahaumer on 3 Feb 2016 16:45 UTC
Reverted changes from e24066a but prepared a possible solution for this issue in b139add.
Usage of matrixVectorProduct seems to be problematic in Modelica.Magnetic.Quasistatic.FundamentalWave.Components.MultiPhaseElectroMagneticConverter.

modelica-trac-importer · 2017-01-15T07:52:24Z

Comment by ahaumer on 3 Feb 2016 18:32 UTC
I did withdraw the ComplexMath.Vector.matrixVectorProduct function, but I replaced the complex matrix times vector product by another formulation (equations) for the 3 occurences.
Resolved by ad8900e; if you still have issues, please reopen the ticket and describe the problems.

modelica-trac-importer · 2017-01-15T07:52:55Z

Comment by beutlich on 3 Feb 2016 21:49 UTC
Thanks for fixing it, Toni!

My only issue is that this is now the least elegant solution and I wonder why neither the proposed matrixVectorProduct nor your even more elegant matrixVectorProduct (using _ComplexMath.'sum' _) did not work out.

modelica-trac-importer · 2017-01-15T07:53:05Z

Comment by beutlich on 4 Feb 2016 09:59 UTC
I reopen for the above comment. The Complex matrix vector product used in Magnetic and Electrical must not know the underlying implementation of the Complex record. This is not the Modelica we really like. Sorry for the double work. I hope you understand my concerns.

modelica-trac-importer · 2017-01-15T07:53:35Z

Comment by ahaumer on 4 Feb 2016 13:10 UTC
I found two solutions that worked - no other solution worked (at least in Dymola):

Matrix * Vector multiplication by the tool
The formulation that is implemented now
I recommend not using function calls in models. The tool has no chance to simplify the equations (unless you propose for inlining). In the MultiPhaseElectroMagneticConverter the tool has to combine the currents to calculate MMF; this is hard if the "combining" is hidden in a function.
Ok you think the matrix vector product must not know the Complex implemenation.
How should one then access real and imaginary part?
On the other hand, the model MultiPhaseElectroMagneticConverter must not know the implementation of the product of a matrix and a vector.
I added comments "preferred implemenation", so we can change later when most tools have done thier homework and implemented the operator overloading from Complex numbers to Complex matrices.
For now: This is the only solution that works in Dymola AND gives other tools the chance to run the models without complex matrix multiplication. -> Close (works for me)

modelica-trac-importer · 2017-01-15T07:54:05Z

Comment by sjoelund.se on 4 Feb 2016 13:21 UTC
Replying to [comment:9 ahaumer]:

I recommend not using function calls in models. The tool has no chance to simplify the equations (unless you propose for inlining). In the MultiPhaseElectroMagneticConverter the tool has to combine the currents to calculate MMF; this is hard if the "combining" is hidden in a function.

Then you also propose not operator overloading since that is the same as function calls.

For now: This is the only solution that works in Dymola AND gives other tools the chance to run the models without complex matrix multiplication. -> Close (works for me)

Function calls did not work in Dymola? This is equivalent to operator overloading in Modelica (the only difference is that the function is manually selected).

modelica-trac-importer · 2017-01-15T07:54:15Z

Comment by beutlich on 4 Feb 2016 13:32 UTC
One-liner algorithms like your proposal from comment:4 or Modelica.ComplexMath.Vectors.length are inlined by default.

modelica-trac-importer · 2017-01-15T07:54:25Z

Modified by dietmarw on 4 Feb 2016 14:41 UTC

modelica-trac-importer · 2017-01-15T07:54:55Z

Comment by ahaumer on 5 Feb 2016 14:24 UTC
Replying to [comment:12 dietmarw]:
So why did you reopen the ticket without any comment?

It was desired to find an alternative for the original formulation.
I tested and found only one alternative that worked in Dymola and OpenModelica - and I'm pretty sure it is easier to handle for other tools.

This problem seems to be solved (at least for me). Maybe we can implement a more elegant solution in future - but then we sould shift the milestone (for sure not 3.2.2).

If you want to have a discussion regarding equations, operator overloading, functions, inlinig, etc. we should open a new ticket. I admit I do not know too much about compiler, but I'd like to know more abput that from a modeler's point of view.

modelica-trac-importer · 2017-01-15T07:55:06Z

Comment by ahaumer on 5 Feb 2016 14:33 UTC
Replying to [comment:11 beutlich]:
Using a one-liner direct in the component itself (not function call) Dymola complains about an algebraic loop involving constants - I don't know why, but I have to clarify that with the developers at DS.

modelica-trac-importer · 2017-01-15T07:55:16Z

Comment by beutlich on 5 Feb 2016 14:42 UTC
OK. Because of all these open questions it does not hurt to have the ticket open until we got the answers.

modelica-trac-importer · 2017-01-15T07:55:26Z

Comment by hansolsson on 18 Feb 2016 13:00 UTC
Replying to [comment:14 ahaumer]:

Replying to [comment:11 beutlich]:
Using a one-liner direct in the component itself (not function call) Dymola complains about an algebraic loop involving constants - I don't know why, but I have to clarify that with the developers at DS.

I was going to ask about the exact problem - and then I realized it:

Even if matrixVectorProduct is a one-line it will not be fully inlined unless
Modelica.ComplexMath.'sum' also is a one-liner, i.e.:

function 'sum' "Return sum of complex vector"
  extends Modelica.Icons.Function;
  input Complex v[:] "Vector";
  output Complex result "Complex sum of vector elements";
algorithm 
  result:=Complex(sum(Modelica.ComplexMath.real(v)), sum(Modelica.ComplexMath.imag(v)));
end sum;

(well, or using v.re if you want shorter code).

modelica-trac-importer · 2017-01-15T07:55:36Z

Comment by beutlich on 18 Feb 2016 13:17 UTC
Note: This also holds for function 'product'.

modelica-trac-importer · 2017-01-15T07:55:46Z

Comment by hansolsson on 18 Feb 2016 15:55 UTC
Replying to [comment:17 beutlich]:

Note: This also holds for function 'product'.

I see that as less important and also more complicated - unless you 'cheat' by using the logarithm to compute the product.

However, for matrixVectorProduct I would still prefer if we just wrote m*v.

modelica-trac-importer · 2017-01-15T07:55:56Z

Comment by beutlich on 18 Feb 2016 16:05 UTC
Agreed, m*v is superior, matrixVectorProduct is nice and current implementation is not elegant.

However, m*v requires a tool to first consider the matrix-vector-product and second the overloaded times operator. Is this according to the MLS?

Can you please give the reasons why we have function scalarProduct in Complex?

modelica-trac-importer · 2017-01-15T07:56:06Z

Comment by otter on 18 Feb 2016 16:19 UTC
Replying to [comment:16 hansolsson]:

Replying to [comment:14 ahaumer]:

Replying to [comment:11 beutlich]:
Using a one-liner direct in the component itself (not function call) Dymola complains about an algebraic loop involving constants - I don't know why, but I have to clarify that with the developers at DS.

I was going to ask about the exact problem - and then I realized it:

Even if matrixVectorProduct is a one-line it will not be fully inlined unless
Modelica.ComplexMath.'sum' also is a one-liner, i.e.:
function 'sum' "Return sum of complex vector"
  extends Modelica.Icons.Function;
  input Complex v[:] "Vector";
  output Complex result "Complex sum of vector elements";
algorithm 
  result:=Complex(sum(Modelica.ComplexMath.real(v)), sum(Modelica.ComplexMath.imag(v)));
end sum;
(well, or using v.re if you want shorter code).

Improved implementation of 'sum' so that it can be inlined (and introducing Inline=true annotation) in 7d7f4ee
(This is fine by itself, but I do not fully understand the explanation because sum(..) is used in model Modelica.Magnetic.QuasiStatic.FundamentalWave.Components.MultiPhaseElectroMagneticConverter and not ComplexMath.'sum').

modelica-trac-importer · 2017-01-15T07:56:36Z

Comment by hansolsson on 18 Feb 2016 16:20 UTC
Replying to [comment:19 beutlich]:

Agreed, m*v is superior, matrixVectorProduct is nice and current implementation is not elegant.

However, m*v requires a tool to first consider the matrix-vector-product and second the overloaded times operator. Is this according to the MLS?

Yes. The rules in M Spec 32 R2 in section 14.4 are that vector*vector and vector*matrix are undefined, but matrix*vector and matrix*matrix are defined - using the same definition as in 10.6.4 (and relying on Complex.'0' for zero inner size).

Can you please give the reasons why we have function scalarProduct in Complex?

Because sum(a[k]*b[k] for k in 1:n) as is used for reals would give the wrong result. We need a good solution for conjugating one of the operands - when appropriate.

(Note: Wikipedia gives a definition that conjugates 'b' for scalar product, but the bra-ket notation conjugates 'a'.)

E.g. Matlab solves it by writing a*b' - where ' is the conjugate transpose.

modelica-trac-importer · 2017-01-15T07:56:46Z

Comment by hansolsson on 18 Feb 2016 16:23 UTC
Replying to [comment:20 otter]:

(This is fine by itself, but I do not fully understand the explanation because sum(..) is used in model Modelica.Magnetic.QuasiStatic.FundamentalWave.Components.MultiPhaseElectroMagneticConverter and not ComplexMath.'sum').

Because matrixVectorProduct uses ComplexMath.'sum', and we needed to inline both functions to get to the scalar level.

modelica-trac-importer · 2017-01-15T07:56:56Z

Comment by otter on 18 Feb 2016 16:49 UTC
Introduced MatrixVectorProduct as suggested and used it in 563e24c. Simulates in Dymola. Used it in

1 in Modelica.Electrical.QuasiStationary.MulitPhase.Blocks.SymmetricalComponents
2 in Modelica.Magnetic.QuasiStatic.FundamentalWave.Components.MultiPhaseElectroMagneticConverter

modelica-trac-importer · 2017-01-15T07:57:07Z

Comment by beutlich on 18 Feb 2016 20:49 UTC
Thank you for fixing. I'd even remove the commented sections though.

modelica-trac-importer · 2017-01-15T07:57:17Z

Comment by otter on 26 Feb 2016 01:05 UTC
Reverted 563e24c in bbff65b. Reason:
The change in 563e24c broke OpenModelica and the reason was not obvious (a simple test model was able to inline the introduced MatrixVectorProduction function).
The reverted changed introduced in bbff65b uses the previous implementation with array comprehensions. It has been tested (with example Modelica.Magnetic.QuasiStatic.FundamentalWave.Examples.Components.MultiPhaseInductance) to work in Dymola 2017 Alpha and in OpenModelica 1.9.4 Beta.2 (previously, this example failed in OpenModelica).

modelica-trac-importer · 2017-01-15T07:57:27Z

Comment by beutlich on 26 Feb 2016 07:32 UTC
For the record: om:#3713 is the corresponding OpenModelica ticket.

modelica-trac-importer · 2017-01-15T07:57:37Z

Comment by beutlich on 29 Feb 2016 22:14 UTC
For me, it rather looks like a workaround to have the array comprehensions three times within the models. We really should rethink what is the most standard-like option.

modelica-trac-importer · 2017-01-15T07:58:07Z

Comment by otter on 1 Mar 2016 07:23 UTC
Replying to [comment:27 beutlich]:

For me, it rather looks like a workaround to have the array comprehensions three times within the models. We really should rethink what is the most standard-like option.

I totally disagree. We have a solution now that works in several tools. The only good alternative is that all tools support Modelica 3.2 fully. It is not worth to force tools to support now overloaded matrix operations on operator records (which is the only fully "clean" solution).

Fill a new ticket, if you think this should be fixed in the future and we change the implementation in a future release once better tool support is available.

git-svn-id: https://svn.modelica.org/projects/Modelica/trunk@9033 7ce873d0-865f-4ce7-a662-4bb36ea78beb

…delica#1890 git-svn-id: https://svn.modelica.org/projects/Modelica/trunk@9034 7ce873d0-865f-4ce7-a662-4bb36ea78beb

git-svn-id: https://svn.modelica.org/projects/Modelica/trunk@9035 7ce873d0-865f-4ce7-a662-4bb36ea78beb

…added in MSL trunk) As suggested, implemented 'sum' so that it can be inlined (and introduced Inline=true annotation) git-svn-id: https://svn.modelica.org/projects/Modelica/trunk@9050 7ce873d0-865f-4ce7-a662-4bb36ea78beb

…added in MSL trunk) As suggested, introduced function MatrixVectorProduct and used it. This works in Dymola. git-svn-id: https://svn.modelica.org/projects/Modelica/trunk@9051 7ce873d0-865f-4ce7-a662-4bb36ea78beb

…added in MSL trunk) Reverted changes from r9051, so the matrix-vector product function is removed and the matrix-vector product of complex arrays is replaced by the previoius array comprehension implementation. The reason is that r9051 did not work in OpenModelica. The reverted change has been tested to work in Dymola and in OpenModelica. git-svn-id: https://svn.modelica.org/projects/Modelica/trunk@9073 7ce873d0-865f-4ce7-a662-4bb36ea78beb

modelica-trac-importer assigned AHaumer Jan 15, 2017

modelica-trac-importer added bug Critical/severe issue L: Complex* Issue addresses Complex, Modelica.ComplexBlocks or Modelica.ComplexMath labels Jan 15, 2017

modelica-trac-importer added this to the MSL3.2.2 milestone Jan 15, 2017