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

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