root/branches/maintenance/3.0/Modelica/Math/package.mo

Revision 1092, 196.3 kB (checked in by otter, 9 months ago)

Math.Matrices.eigenValues/rank/singularValues lead to a crash, for matrices with zero dimensions.
This has been fixed
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
Line 
1within Modelica;
2package Math "Library of mathematical functions (e.g., sin, cos) and of functions operating on vectors and matrices"
3  import SI = Modelica.SIunits;
4
5
6extends Modelica.Icons.Library2;
7
8
9annotation (
10  Invisible=true,
11  Icon(coordinateSystem(preserveAspectRatio=true, extent={{-100,-100},{100,100}}), 
12      graphics={Text(
13        extent={{-59,-9},{42,-56}},
14        lineColor={0,0,0},
15        textString="f(x)")}),
16  Documentation(info="<HTML>
17<p>
18This package contains <b>basic mathematical functions</b> (such as sin(..)),
19as well as functions operating on <b>vectors</b> and <b>matrices</b>.
20</p>
21 
22<dl>
23<dt><b>Main Author:</b>
24<dd><a href=\"http://www.robotic.dlr.de/Martin.Otter/\">Martin Otter</a><br>
25    Deutsches Zentrum f&uuml;r Luft und Raumfahrt e.V. (DLR)<br>
26    Institut f&uuml;r Robotik und Mechatronik<br>
27    Postfach 1116<br>
28    D-82230 Wessling<br>
29    Germany<br>
30    email: <A HREF=\"mailto:Martin.Otter@dlr.de\">Martin.Otter@dlr.de</A><br>
31</dl>
32 
33<p>
34Copyright &copy; 1998-2008, Modelica Association and DLR.
35</p>
36<p>
37<i>This Modelica package is <b>free</b> software; it can be redistributed and/or modified
38under the terms of the <b>Modelica license</b>, see the license conditions
39and the accompanying <b>disclaimer</b> 
40<a href=\"Modelica://Modelica.UsersGuide.ModelicaLicense\">here</a>.</i>
41</p><br>
42</HTML>
43", revisions="<html>
44<ul>
45<li><i>October 21, 2002</i>
46       by <a href=\"http://www.robotic.dlr.de/Martin.Otter/\">Martin Otter</a>
47       and <a href=\"http://www.robotic.dlr.de/Christian.Schweiger/\">Christian Schweiger</a>:<br>
48       Function tempInterpol2 added.</li>
49<li><i>Oct. 24, 1999</i>
50       by <a href=\"http://www.robotic.dlr.de/Martin.Otter/\">Martin Otter</a>:<br>
51       Icons for icon and diagram level introduced.</li>
52<li><i>June 30, 1999</i>
53       by <a href=\"http://www.robotic.dlr.de/Martin.Otter/\">Martin Otter</a>:<br>
54       Realized.</li>
55</ul>
56 
57</html>"));
58
59
60package Vectors "Library of functions operating on vectors"
61  extends Modelica.Icons.Library;
62
63  annotation (
64    preferedView = "info",
65    Documentation(info="<HTML>
66<h4>Library content</h4>
67<p>
68This library provides functions operating on vectors:
69</p>
70<table border=1 cellspacing=0 cellpadding=2>
71  <tr><th><i>Function</i></th>
72      <th><i>Description</i></th>
73  </tr>
74  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Vectors.isEqual\">isEqual</a>(v1, v2)</td>
75      <td valign=\"top\">Determines whether two vectors have the same size and elements</td>
76  </tr>
77  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Vectors.norm\">norm</a>(v,p)</td>
78      <td valign=\"top\">p-norm of vector v</td>
79  </tr>
80  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Vectors.length\">length</a>(v)</td>
81      <td valign=\"top\">Length of vector v (= norm(v,2), but inlined and therefore usable in
82          symbolic manipulations) </td>
83  </tr>
84  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Vectors.normalize\">normalize</a>(v)</td>
85      <td valign=\"top\">Return normalized vector such that length = 1 and prevent
86          zero-division for zero vector</td>
87  </tr>
88  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Vectors.reverse\">reverse</a>(v)</td>
89      <td valign=\"top\">Reverse vector elements</td>
90  </tr>
91  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Vectors.sort\">sort</a>(v)</td>
92      <td valign=\"top\">Sort elements of vector in ascending or descending order</td>
93  </tr>
94</table>
95<h4>See also</h4>
96<a href=\"Modelica://Modelica.Math.Matrices\">Matrices</a>
97</HTML>"));
98
99  function isEqual "Determine if two Real vectors are numerically identical"
100    extends Modelica.Icons.Function;
101    input Real v1[:] "First vector";
102    input Real v2[:] "Second vector (may have different length as v1";
103    input Real eps(min=0) = 0
104      "Two elements e1 and e2 of the two vectors are identical if abs(e1-e2) <= eps";
105    output Boolean result
106      "= true, if vectors have the same length and the same elements";
107
108    annotation ( Documentation(info="<HTML>
109<h4>Syntax</h4>
110<blockquote><pre>
111Vectors.<b>isEqual</b>(v1, v2);
112Vectors.<b>isEqual</b>(v1, v2, eps=0);
113</pre></blockquote>
114<h4>Description</h4>
115<p>
116The function call \"<code>Vectors.isEqual(v1, v2)</code>\" returns <b>true</b>,
117if the two Real vectors v1 and v2 have the same dimensions and
118the same elements. Otherwise the function
119returns <b>false</b>. Two elements e1 and e2 of the two vectors
120are checked on equality by the test \"abs(e1-e2) &le; eps\", where \"eps\"
121can be provided as third argument of the function. Default is \"eps = 0\".
122</p>Modelica.Utilities.Strings.isEqual
123<h4>Example</h4>
124<blockquote><pre>
125  Real v1[3] = {1, 2, 3};
126  Real v2[3] = {1, 2, 3, 4};
127  Real v3[3] = {1, 2, 3.0001};
128  Boolean result;
129<b>algorithm</b>
130  result := Vectors.isEqual(v1,v2);     // = <b>false</b>
131  result := Vectors.isEqual(v1,v3);     // = <b>false</b>
132  result := Vectors.isEqual(v1,v1);     // = <b>true</b>
133  result := Vectors.isEqual(v1,v3,0.1); // = <b>true</b>
134</pre></blockquote>
135<h4>See also</h4>
136<a href=\"Modelica://Modelica.Math.Matrices.isEqual\">Matrices.isEqual</a>,
137<a href=\"Modelica://Modelica.Utilities.Strings.isEqual\">Strings.isEqual</a>
138</HTML>"));
139  protected
140    Integer n=size(v1, 1) "Dimension of vector v1";
141    Integer i=1;
142  algorithm
143    result := false;
144    if size(v2, 1) == n then
145      result := true;
146      while i <= n loop
147        if abs(v1[i] - v2[i]) > eps then
148          result := false;
149          i := n;
150        end if;
151        i := i + 1;
152      end while;
153    end if;
154  end isEqual;
155
156  function norm "Return the p-norm of a vector"
157    extends Modelica.Icons.Function;
158    input Real v[:] "Vector";
159    input Real p(min=1) = 2
160      "Type of p-norm (often used: 1, 2, or Modelica.Constants.inf)";
161    output Real result "p-norm of vector v";
162
163    annotation ( Documentation(info="<HTML>
164<h4>Syntax</h4>
165<blockquote><pre>
166Vectors.<b>norm</b>(v);
167Vectors.<b>norm</b>(v,p=2);   // 1 &le; p &le; &#8734;
168</pre></blockquote>
169<h4>Description</h4>
170<p>
171The function call \"<code>Vectors.<b>norm</b>(v)</code>\" returns the
172<b>Euclidean norm</b> \"<code>sqrt(v*v)</code>\" of vector v.
173With the optional
174second argument \"p\", any other p-norm can be computed:
175</p>
176<center>
177<IMG SRC=\"../Images/Math/vectorNorm.png\" ALT=\"function Vectors.norm\">
178</center>
179<p>
180Besides the Euclidean norm (p=2), also the 1-norm and the
181infinity-norm are sometimes used:
182</p>
183<table border=1 cellspacing=0 cellpadding=2>
184  <tr><td valign=\"top\"><b>1-norm</b></td>
185      <td valign=\"top\">= sum(abs(v))</td>
186      <td valign=\"top\"><b>norm</b>(v,1)</td>
187  </tr>
188  <tr><td valign=\"top\"><b>2-norm</b></td>
189      <td valign=\"top\">= sqrt(v*v)</td>
190      <td valign=\"top\"><b>norm</b>(v) or <b>norm</b>(v,2)</td>
191  </tr>
192  <tr><td valign=\"top\"><b>infinity-norm</b></td>
193      <td valign=\"top\">= max(abs(v))</td>
194      <td valign=\"top\"><b>norm</b>(v,Modelica.Constants.<b>inf</b>)</td>
195  </tr>
196</table>
197<p>
198Note, for any vector norm the following inequality holds:
199</p>
200<blockquote><pre>
201<b>norm</b>(v1+v2,p) &le; <b>norm</b>(v1,p) + <b>norm</b>(v2,p)
202</pre></blockquote>
203<h4>Example</h4>
204<blockquote><pre>
205  v = {2, -4, -2, -1};
206  <b>norm</b>(v,1);    // = 9
207  <b>norm</b>(v,2);    // = 5
208  <b>norm</b>(v);      // = 5
209  <b>norm</b>(v,10.5); // = 4.00052597412635
210  <b>norm</b>(v,Modelica.Constants.inf);  // = 4
211</pre></blockquote>
212<h4>See also</h4>
213<a href=\"Modelica://Modelica.Math.Matrices.norm\">Matrices.norm</a>
214</HTML>"));
215  algorithm
216    if p == 2 then
217      result:=sqrt(v*v);
218    elseif p == Modelica.Constants.inf then
219      result:=max(abs(v));
220    elseif p == 1 then
221      result:=sum(abs(v));
222    else
223      result:=(sum(abs(v[i])^p for i in 1:size(v, 1)))^(1/p);
224    end if;
225  end norm;
226
227  function length
228    "Return length of a vectorReturn length of a vector (better as norm(), if further symbolic processing is performed)"
229    extends Modelica.Icons.Function;
230    input Real v[:] "Vector";
231    output Real result "Length of vector v";
232    annotation (Inline=true, Documentation(info="<html>
233<h4>Syntax</h4>
234<blockquote><pre>
235Vectors.<b>length</b>(v);
236</pre></blockquote>
237<h4>Description</h4>
238<p>
239The function call \"<code>Vectors.<b>length</b>(v)</code>\" returns the
240<b>Euclidean length</b> \"<code>sqrt(v*v)</code>\" of vector v.
241The function call is equivalent to Vectors.norm(v). The advantage of
242length(v) over norm(v)\"is that function length(..) is implemented
243in one statement and therefore the function is usually automatically
244inlined. Further symbolic processing is therefore possible, which is
245not the case with function norm(..).
246</p>
247<h4>Example</h4>
248<blockquote><pre>
249  v = {2, -4, -2, -1};
250  <b>length</b>(v);  // = 5
251</pre></blockquote>
252<h4>See also</h4>
253<a href=\"Modelica://Modelica.Math.Vectors.norm\">Vectors.norm</a>
254</html>"));
255  algorithm
256    result := sqrt(v*v);
257  end length;
258
259  function normalize
260    "Return normalized vector such that length = 1Return normalized vector such that length = 1 and prevent zero-division for zero vector"
261    extends Modelica.Icons.Function;
262    input Real v[:] "Vector";
263    input Real eps = 100*Modelica.Constants.eps
264      "if |v| < eps then result = v/eps";
265    output Real result[size(v, 1)] "Input vector v normalized to length=1";
266
267    annotation (Inline=true, Documentation(info="<html>
268<h4>Syntax</h4>
269<blockquote><pre>
270Vectors.<b>normalize</b>(v);
271Vectors.<b>normalize</b>(v,eps=100*Modelica.Constants.eps);
272</pre></blockquote>
273<h4>Description</h4>
274<p>
275The function call \"<code>Vectors.<b>normalize</b>(v)</code>\" returns the
276<b>unit vector</b> \"<code>v/length(v)</code>\" of vector v.
277If length(v) is close to zero (more precisely, if length(v) &lt; eps),
278v/eps is returned in order to avoid
279a division by zero. For many applications this is useful, because
280often the unit vector <b>e</b> = <b>v</b>/length(<b>v</b>) is used to compute
281a vector x*<b>e</b>, where the scalar x is in the order of length(<b>v</b>),
282i.e., x*<b>e</b> is small, when length(<b>v</b>) is small and then
283it is fine to replace <b>e</b> by <b>v</b> to avoid a division by zero.
284</p>
285<p>
286Since the function is implemented in one statement,
287it is usually inlined and therefore symbolic processing is
288possible.
289</p>
290<h4>Example</h4>
291<blockquote><pre>
292  <b>normalize</b>({1,2,3});  // = {0.267, 0.534, 0.802}
293  <b>normalize</b>({0,0,0});  // = {0,0,0}
294</pre></blockquote>
295<h4>See also</h4>
296<a href=\"Modelica://Modelica.Math.Vectors.length\">Vectors.length</a>
297</html>"));
298  algorithm
299    result := smooth(0,if length(v) >= eps then v/length(v) else v/eps);
300  end normalize;
301
302  function reverse "Reverse vector elements (e.g. v[1] becomes last element)"
303    extends Modelica.Icons.Function;
304    input Real v[:] "Vector";
305    output Real result[size(v, 1)] "Elements of vector v in reversed order";
306
307    annotation (Inline=true, Documentation(info="<html>
308<h4>Syntax</h4>
309<blockquote><pre>
310Vectors.<b>reverse</b>(v);
311</pre></blockquote>
312<h4>Description</h4>
313<p>
314The function call \"<code>Vectors.<b>reverse</b>(v)</code>\" returns the
315vector elements in reverse order.
316</p>
317<h4>Example</h4>
318<blockquote><pre>
319  <b>reverse</b>({1,2,3,4});  // = {4,3,2,1}
320</pre></blockquote>
321</html>"));
322  algorithm
323    result := {v[end-i+1] for i in 1:size(v,1)};
324  end reverse;
325
326  function sort "Sort elements of vector in ascending or descending order"
327    extends Modelica.Icons.Function;
328    input Real v[:] "Vector to be sorted";
329    input Boolean ascending = true
330      "= true if ascending order, otherwise descending order";
331    output Real sorted_v[size(v,1)] = v "Sorted vector";
332    output Integer indices[size(v,1)] = 1:size(v,1) "sorted_v = v[indices]";
333
334    annotation (Documentation(info="<HTML>
335<h4>Syntax</h4>
336<blockquote><pre>
337           sorted_v = Vectors.<b>sort</b>(v);
338(sorted_v, indices) = Vectors.<b>sort</b>(v, ascending=true);
339</pre></blockquote>
340<h4>Description</h4>
341<p>
342Function <b>sort</b>(..) sorts a Real vector v
343in ascending order and returns the result in sorted_v.
344If the optional argument \"ascending\" is <b>false</b>, the vector
345is sorted in descending order. In the optional second
346output argument the indices of the sorted vector with respect
347to the original vector are given, such that sorted_v = v[indices].
348</p>
349<h4>Example</h4>
350<blockquote><pre>
351  (v2, i2) := Vectors.sort({-1, 8, 3, 6, 2});
352       -> v2 = {-1, 2, 3, 6, 8}
353          i2 = {1, 5, 3, 4, 2}
354</pre></blockquote>
355</HTML>"));
356    /* shellsort algorithm; should be improved later */
357  protected
358    Integer gap;
359    Integer i;
360    Integer j;
361    Real wv;
362    Integer wi;
363    Integer nv = size(v,1);
364    Boolean swap;
365  algorithm
366    gap := div(nv,2);
367
368    while gap > 0 loop
369       i := gap;
370       while i < nv loop
371          j := i-gap;
372          if j>=0 then
373             if ascending then
374                swap := sorted_v[j+1] > sorted_v[j + gap + 1];
375             else
376                swap := sorted_v[j+1] < sorted_v[j + gap + 1];
377             end if;
378          else
379             swap := false;
380          end if;
381
382          while swap loop
383             wv := sorted_v[j+1];
384             wi := indices[j+1];
385             sorted_v[j+1] := sorted_v[j+gap+1];
386             sorted_v[j+gap+1] := wv;
387             indices[j+1] := indices[j+gap+1];
388             indices[j+gap+1] := wi;
389             j := j - gap;
390             if j >= 0 then
391                if ascending then
392                   swap := sorted_v[j+1] > sorted_v[j + gap + 1];
393                else
394                   swap := sorted_v[j+1] < sorted_v[j + gap + 1];
395                end if;
396             else
397                swap := false;
398             end if;
399          end while;
400          i := i + 1;
401       end while;
402       gap := div(gap,2);
403    end while;
404  end sort;
405end Vectors;
406
407
408package Matrices "Library of functions operating on matrices"
409
410  extends Modelica.Icons.Library;
411
412  annotation (
413    version="0.8.1",
414    versionDate="2004-08-21",
415    Documentation(info="<HTML>
416<h4>Library content</h4>
417<p>
418This library provides functions operating on matrices:
419</p>
420<table border=1 cellspacing=0 cellpadding=2>
421  <tr><th><i>Function</i></th>
422      <th><i>Description</i></th>
423  </tr>
424  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.isEqual\">isEqual</a>(M1, M2)</td>
425      <td valign=\"top\">Determines whether two matrices have the same size and elements</td>
426  </tr>
427  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.norm\">norm</a>(A)</td>
428      <td valign=\"top\">1-, 2- and infinity-norm of matrix A</td>
429  </tr>
430  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.sort\">sort</a>(M)</td>
431      <td valign=\"top\">Sort rows or columns of matrix in ascending or descending order</td>
432  </tr>
433  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.solve\">solve</a>(A,b)</td>
434      <td valign=\"top\">Solve real system of linear equations A*x=b with a b vector</td>
435  </tr>
436  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.solve2\">solve2</a>(A,B)</td>
437      <td valign=\"top\">Solve real system of linear equations A*X=B with a B matrix</td>
438  </tr>
439  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.leastSquares\">leastSquares</a>(A,b)</td>
440      <td valign=\"top\">Solve overdetermined or underdetermined real system of <br>
441          linear equations A*x=b in a least squares sense</td>
442  </tr>
443  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.equalityLeastSquares\">equalityLeastSquares</a>(A,a,B,b)</td>
444      <td valign=\"top\">Solve a linear equality constrained least squares problem:<br>
445          min|A*x-a|^2 subject to B*x=b</td>
446  </tr>
447  <tr><td valign=\"top\">(LU,p,info) = <a href=\"Modelica://Modelica.Math.Matrices.LU\">LU</a>(A)</td>
448      <td valign=\"top\">LU decomposition of square or rectangular matrix</td>
449  </tr>
450  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.LU_solve\">LU_solve</a>(LU,p,b)</td>
451      <td valign=\"top\">Solve real system of linear equations P*L*U*x=b with a<br>
452          b vector and an LU decomposition from \"LU(..)\"</td>
453  </tr>
454  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.LU_solve2\">LU_solve2</a>(LU,p,B)</td>
455      <td valign=\"top\">Solve real system of linear equations P*L*U*X=B with a<br>
456          B matrix and an LU decomposition from \"LU(..)\"</td>
457  </tr>
458  <tr><td valign=\"top\">(Q,R,p) = <a href=\"Modelica://Modelica.Math.Matrices.QR\">QR</a>(A)</td>
459      <td valign=\"top\"> QR decomposition with column pivoting of rectangular matrix (Q*R = A[:,p]) </td>
460  </tr>
461  <tr><td valign=\"top\">eval = <a href=\"Modelica://Modelica.Math.Matrices.eigenValues\">eigenValues</a>(A)<br>
462          (eval,evec) = <a href=\"Modelica://Modelica.Math.Matrices.eigenValues\">eigenValues</a>(A)</td>
463      <td valign=\"top\"> Compute eigenvalues and optionally eigenvectors<br>
464           for a real, nonsymmetric matrix </td>
465  </tr>
466  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.eigenValueMatrix\">eigenValueMatrix</a>(eigen)</td>
467      <td valign=\"top\"> Return real valued block diagonal matrix J of eigenvalues of
468            matrix A (A=V*J*Vinv) </td>
469  </tr>
470  <tr><td valign=\"top\">sigma = <a href=\"Modelica://Modelica.Math.Matrices.singularValues\">singularValues</a>(A)<br>
471      (sigma,U,VT) = <a href=\"Modelica://Modelica.Math.Matrices.singularValues\">singularValues</a>(A)</td>
472      <td valign=\"top\"> Compute singular values and optionally left and right singular vectors </td>
473  </tr>
474  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.det\">det</a>(A)</td>
475      <td valign=\"top\"> Determinant of a matrix (do <b>not</b> use; use rank(..))</td>
476  </tr>
477  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.inv\">inv</a>(A)</td>
478      <td valign=\"top\"> Inverse of a matrix </td>
479  </tr>
480  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.rank\">rank</a>(A)</td>
481      <td valign=\"top\"> Rank of a matrix </td>
482  </tr>
483  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.balance\">balance</a>(A)</td>
484      <td valign=\"top\">Balance a square matrix to improve the condition</td>
485  </tr>
486  <tr><td valign=\"top\"><a href=\"Modelica://Modelica.Math.Matrices.exp\">exp</a>(A)</td>
487      <td valign=\"top\"> Compute the exponential of a matrix by adaptive Taylor series<br> 
488           expansion with scaling and balancing</td>
489  </tr>
490  <tr><td valign=\"top\">(P, G) = <a href=\"Modelica://Modelica.Math.Matrices.integralExp\">integralExp</a>(A,B)</td>
491      <td valign=\"top\"> Compute the exponential of a matrix and its integral</td>
492  </tr>
493  <tr><td valign=\"top\">(P, G, GT) = <a href=\"Modelica://Modelica.Math.Matrices.integralExpT\">integralExpT</a>(A,B)</td>
494      <td valign=\"top\"> Compute the exponential of a matrix and two integrals</td>
495  </tr>
496</table>
497 
498<p>
499Most functions are solely an interface to the external LAPACK library
500(<a href=\"http://www.netlib.org/lapack\">http://www.netlib.org/lapack</a>).
501The details of this library are described in:
502</p>
503 
504<dl>
505<dt>Anderson E., Bai Z., Bischof C., Blackford S., Demmel J., Dongarra J.,
506    Du Croz J., Greenbaum A., Hammarling S., McKenney A., and Sorensen D.:</dt>
507<dd> <b>Lapack Users' Guide</b>.
508     Third Edition, SIAM, 1999.</dd>
509</dl>
510 
511<h4>See also</h4>
512<a href=\"Modelica://Modelica.Math.Vectors\">Vectors</a>
513
514</HTML>
515"));
516
517  function isEqual "Compare whether two Real matrices are identical"
518    extends Modelica.Icons.Function;
519    input Real M1[:, :] "First matrix";
520    input Real M2[:, :] "Second matrix (may have different size as M1";
521    input Real eps(min=0) = 0
522      "Two elements e1 and e2 of the two matrices are identical if abs(e1-e2) <= eps";
523    output Boolean result
524      "= true, if matrices have the same size and the same elements";
525
526    annotation ( Documentation(info="<HTML>
527<h4>Syntax</h4>
528<blockquote><pre>
529Matrices.<b>isEqual</b>(M1, M2);
530Matrices.<b>isEqual</b>(M1, M2, eps=0);
531</pre></blockquote>
532<h4>Description</h4>
533<p>
534The function call \"<code>Matrices.isEqual(M1, M2)</code>\" returns <b>true</b>,
535if the two Real matrices M1 and M2 have the same dimensions and
536the same elements. Otherwise the function
537returns <b>false</b>. Two elements e1 and e2 of the two matrices
538are checked on equality by the test \"abs(e1-e2) &le; eps\", where \"eps\"
539can be provided as third argument of the function. Default is \"eps = 0\".
540</p>
541<h4>Example</h4>
542<blockquote><pre>
543  Real A1[2,2] = [1,2; 3,4];
544  Real A2[3,2] = [1,2; 3,4; 5,6];
545  Real A3[2,2] = [1,2, 3,4.0001];
546  Boolean result;
547<b>algorithm</b>
548  result := Matrices.isEqual(M1,M2);     // = <b>false</b>
549  result := Matrices.isEqual(M1,M3);     // = <b>false</b>
550  result := Matrices.isEqual(M1,M1);     // = <b>true</b>
551  result := Matrices.isEqual(M1,M3,0.1); // = <b>true</b>
552</pre></blockquote>
553<h4>See also</h4>
554<a href=\"Modelica://Modelica.Math.Vectors.isEqual\">Vectors.isEqual</a>,
555<a href=\"Modelica://Modelica.Utilities.Strings.isEqual\">Strings.isEqual</a>
556</HTML>"));
557  protected
558    Integer nrow=size(M1, 1) "Number of rows of matrix M1";
559    Integer ncol=size(M1, 2) "Number of columns of matrix M1";
560    Integer i=1;
561    Integer j;
562  algorithm
563    result := false;
564    if size(M2, 1) == nrow and size(M2, 2) == ncol then
565      result := true;
566      while i <= nrow loop
567        j := 1;
568        while j <= ncol loop
569          if abs(M1[i, j] - M2[i, j]) > eps then
570            result := false;
571            i := nrow;
572            j := ncol;
573          end if;
574          j := j + 1;
575        end while;
576        i := i + 1;
577      end while;
578    end if;
579
580  end isEqual;
581
582  function norm "Returns the norm of a matrix"
583    extends Modelica.Icons.Function;
584    input Real A[:, :] "Input matrix";
585    input Real p(min=1) = 2
586      "Type of p-norm (only allowed: 1, 2 or Modelica.Constants.inf)";
587    output Real result=0.0 "p-norm of matrix A";
588
589    annotation ( Documentation(info="<HTML>
590<h4>Syntax</h4>
591<blockquote><pre>
592Matrices.<b>norm</b>(A);
593Matrices.<b>norm</b>(A, p=2);
594</pre></blockquote>
595<h4>Description</h4>
596<p>
597The function call \"<code>Matrices.norm(A)</code>\" returns the
5982-norm of matrix A, i.e., the largest singular value of A.<br>
599The function call \"<code>Matrices.norm(A, p)</code>\" returns the
600p-norm of matrix A. The only allowed values for p are</p>
601<ul>
602<li> \"p=1\": the largest column sum of A</li>
603<li> \"p=2\": the largest singular value of A</li> 
604<li> \"p=Modelica.Constants.inf\": the largest row sum of A</li>
605</ul>
606<p>
607Note, for any matrices A1, A2 the following inequality holds:
608</p>
609<blockquote><pre>
610Matrices.<b>norm</b>(A1+A2,p) &le; Matrices.<b>norm</b>(A1,p) + Matrices.<b>norm</b>(A2,p)
611</pre></blockquote>
612<p>
613Note, for any matrix A and vector v the following inequality holds:
614</p>
615<blockquote><pre>
616Vectors.<b>norm</b>(A*v,p) &le; Matrices.<b>norm</b>(A,p)*Vectors.<b>norm</b>(A,p)
617</pre></blockquote>
618</HTML>"));
619  algorithm
620    if p == 1 then
621      // column sum norm
622      for i in 1:size(A, 2) loop
623        result := max(result, sum(abs(A[:, i])));
624      end for;
625    elseif p == 2 then
626      // largest singular value
627      result := max(singularValues(A));
628    elseif p == Modelica.Constants.inf then
629      // row sum norm
630      for i in 1:size(A, 1) loop
631        result := max(result, sum(abs(A[i, :])));
632      end for;
633    else
634      assert(false, "Optional argument \"p\" of function \"norm\" must be
6351, 2 or Modelica.Constants.inf");
636    end if;
637  end norm;
638
639  function sort
640    "Sort rows or columns of matrix in ascending or descending order"
641    extends Modelica.Icons.Function;
642    input Real M[:,:] "Matrix to be sorted";
643    input Boolean sortRows = true
644      "= true if rows are sorted, otherwise columns";
645    input Boolean ascending = true
646      "= true if ascending order, otherwise descending order";
647    output Real sorted_M[size(M,1), size(M,2)] = M "Sorted matrix";
648    output Integer indices[if sortRows then size(M,1) else size(M,2)]
649      "sorted_M = if sortRows then M[indices,:] else M[:,indices]";
650
651    annotation (Documentation(info="<HTML>
652<h4>Syntax</h4>
653<blockquote><pre>
654           sorted_M = Matrices.<b>sort</b>(M);
655(sorted_M, indices) = Matrices.<b>sort</b>(M, sortRows=true, ascending=true);
656</pre></blockquote>
657<h4>Description</h4>
658<p>
659Function <b>sort</b>(..) sorts the rows of a Real matrix M
660in ascending order and returns the result in sorted_M.
661If the optional argument \"sortRows\" is <b>false</b>, the columns
662of the matrix are sorted.
663If the optional argument \"ascending\" is <b>false</b>, the rows or
664columns are sorted in descending order. In the optional second
665output argument, the indices of the sorted rows or columns with respect
666to the original matrix are given, such that
667</p>
668<pre>
669   sorted_M = <b>if</b> sortedRow <b>then</b> M[indices,:] <b>else</b> M[:,indices];
670</pre>
671<h4>Example</h4>
672<blockquote><pre>
673  (M2, i2) := Matrices.sort([2, 1,  0;
674                             2, 0, -1]);
675       -> M2 = [2, 0, -1;
676                2, 1, 0 ];
677          i2 = {2,1};
678</pre></blockquote>
679</HTML>"));
680    /* shellsort algorithm; should be improved later */
681  protected
682    Integer gap;
683    Integer i;
684    Integer j;
685    Real wM2[size(M,2)];
686    Integer wi;
687    Integer nM1 = size(M,1);
688    Boolean swap;
689    Real sorted_MT[size(M,2), size(M,1)];
690
691  encapsulated function greater "Compare whether vector v1 > v2"
692      import Modelica;
693    extends Modelica.Icons.Function;
694      import Modelica.Utilities.Types.Compare;
695    input Real v1[:];
696    input Real v2[size(v1,1)];
697    output Boolean result;
698    protected
699    Integer n = size(v1,1);
700    Integer