## Modelica_LinearSystems2.Math.Matrices

Additional functions for Modelica.Math

### Information

Extends from Modelica.Icons.Package (Icon for standard packages).

### Package Content

Name | Description
--- | ---
LAPACK |
Examples |
care | Solution of continuous-time algebraic Riccati equations
dare | Solution of discrete-time algebraic Riccati equations
sylvester | Solution of continuous-time Sylvester equation A*X + X*B = C
dsylvester | Solution of discrete-time Sylvester equation A*X*B + sgn*X = C
lyapunov | Solution of continuous-time Lyapunov equation X*A + A'*X = C
dlyapunov | Solution of continuous-time Lyapunov equation A'X*A - X = C
conditionNumber | Calculate the condition number norm(A)*norm(inv(A))
rcond | Reciprocal condition number
det | Determinant of a matrix (computed by LU decomposition)
norm | Returns the norm of a matrix
nullspace | Orthonormal nullspace of a matrix
eigenValues | Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix
equalityLeastSquares | Solve a linear equality constrained least squares problem
solve | Solve real system of linear equations A*x=b with a b vector (Gaussian elemination with partial pivoting)
solve2 | Solve real system of linear equations A*X=B with a B matrix (Gaussian elemination with partial pivoting)
solve2r | Solve real system of linear equations X*op(A)=B with a B matrix (Gaussian elemination with partial pivoting)
LU | LU decomposition of square or rectangular matrix
LU_solve | Solve real system of linear equations P*L*U*x=b with a b vector and an LU decomposition (from LU(..))
LU_solve2 | Solve real system of linear equations P*L*U*X=B with a B vector and an LU decomposition (from LU(..))
leastSquares | Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient)
leastSquares2 | Solve overdetermined or underdetermined real system of linear equations A*X=B in a least squares sense (A may be rank deficient)
QR | QR decomposition of a rectangular matrix without column pivoting (A = Q*R). Return the full square Q-matrix
LQ | LQ decomposition of a rectangular matrix without column pivoting (A = L*Q)
orthogonalQ | generates a real orthogonal matrix Q which is defined as the product of IHI-ILO elementary reflectors of order N, as returned by DGEHRD
cholesky | Compute the Cholesky factorization of a symmetric positive definte matrix
choleskyUpDate | Compute the cholesky factor Lu according to Au=Lu'*Lu=A + v*v' with A=L'*L
choleskyDownDate | Compute the cholesky factor Ld according to Ad=Ld'*Ld=A - v*v' with A=L'*L
rsf | Computes the real Schur form (RSF) of a square matrix
rsf2 | Computes the real Schur form (RSF) of a square matrix but uses lapack.dgees
hessenberg | Transform a matrix to upper Hessenberg form
toUpperHessenberg | transform a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H
householderReflexion | Reflect each of the vectors ai of matrix A=[a1, a2, ..., an] on a plane with orthogonal vector u
householderSimilarityTransformation | Calculate the similarity transformation S*A*S of matrix A with symmetric householder matrix S = I - 2u*u'
fliplr | flip the columns of a matrix in left/right direction
flipud | flip the columns of a matrix in up/down direction
triangle | Return the upper/lower triangular part of a square matrix
trace | tarce(A) is the sum of the diagonal elements of A
generalizedEigenvaluesTriangular | Compute invariant zeros of linear state space system with a generalized system matrix [A, B, C, D] which is of upper Hessenberg form
printMatrix | print matrix
fromFile | Read matrix from a matlab file

## Modelica_LinearSystems2.Math.Matrices.care

Solution of continuous-time algebraic Riccati equations

### Information

Function care computes the solution X of the continuous-time algebraic Riccati equation

>      Q + A'*X + X*A - X*G*X = 0

with

>            -1
>     G = B*R *B'

using the Schur vector approach proposed by Laub [1].

It is assumed that Q is symmetric and positve semidefinite and R is symmetric, nonsingular and positive definite, (A,B) is stabilizable and (A,Q) is detectable.

The assumptions are not checked in this function

The assumptions guarantee that Hamiltonian matrix

>     H = [A, -G; -Q, -A']

has no pure imaginary eigenvalue and can be put to an ordered real Schur form

>     U'*H*U = S = [S11, S12; 0, S22]

with orthogonal similarity transformation U. S is ordered in such a way, that S11 contains the n stable eigenvalues of the closed loop system with system matrix

>            -1
>     A - B*R *B'*X

If U is partitioned to

>     U = [U11, U12; U21, U22]

with dimenstions according to S, the solution X can be calculated by

>     X*U11 = U21.

The algorithm uses LAPACK routines dgehrd (to compute the upper Hessenberg matrix of H), dorghr (to calculate the orthogonal matrix from the elementary reflectors as returned from dgehrd), dhseqr (to put transformed H to Schur form and to calculate the eigenvalues of the closed loop system) and dtrsen (to compute the ordered real Schur form and matrix U).

References

      [1] Laub, A.J.
          A Schur Method for Solving Algebraic Riccati equations.
          IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] |
B[size(A, 1), :] |
R[size(B, 2), size(B, 2)] |
Q[size(A, 1), size(A, 1)] |
refine |

### Outputs

Name | Description
--- | ---
X[size(A, 1), size(A, 2)] | stabilizing solution of CARE
ev[size(A, 1)] | eigenvalues of the closed loop system

## Modelica_LinearSystems2.Math.Matrices.dare

Solution of discrete-time algebraic Riccati equations

### Information

Function dare computes the solution X of the discrete-time algebraic Riccati equation

>                                      -1
>      X = A'*X*A - A'*X*B*(R + B'*X*B)  *B'*X*A + Q

using the Schur vector approach proposed by Laub [1].

It is assumed that Q is symmetric and positve semidefinite and R is symmetric, nonsingular and positive definite, (A,B) is stabilizable and (A,Q) is detectable.

The assumptions are not checked in this function

The assumptions guarantee that the Hamiltonian matrix

>           -1  -1       -1        -1
>     H = [A, -A  *G; Q*A, A' + Q*A  *G ]

with

>            -1
>     G = B*R *B'

has no eigenvalue on the unit circle and can be put to an ordered real Schur form

>     U'*H*U = X = [S11, S12; 0, S22]

with orthogonal similarity transformation U. X is ordered in such a way, that S11 contains the n stable eigenvalues of the closed loop system with system matrix

>                       -1
>     A - B*(R + B'*X*B)  *B'*X*A

If U is partitioned to

>     U = [U11, U12; U21, U22]

according to X, the solution X can be calculated by

>     X*U11 = U21.

The algorithm uses LAPACK routines dgehrd (to compute the upper Hessenberg matrix of H), dorghr (to calculate the orthogonal matrix from the elementary reflectors as returned from dgehrd), dhseqr (to put transformed H to Schur form and to calculate the eigenvalues of the closed loop system) and dtrsen (to compute the ordered real Schur form and matrix U).

References

      [1] Laub, A.J.
          A Schur Method for Solving Algebraic Riccati equations.
          IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] |
B[size(A, 1), :] |
R[size(B, 2), size(B, 2)] |
Q[size(A, 1), size(A, 1)] |
refine |

### Outputs

Name | Description
--- | ---
X[size(A, 1), size(A, 2)] | orthogonal matrix of the Schur vectors associated to ordered rsf
ev[size(A, 1)] | eigenvalues of the closed loop system

## Modelica_LinearSystems2.Math.Matrices.sylvester

Solution of continuous-time Sylvester equation A*X + X*B = C

### Information

Function sylvester computes the solution X of the continuous-time Sylvester equation

>      A*X + X*B = C.

using the Schur method for Sylvester equations proposed by Bartels and Stewart [1].

References

      [1] Bartels, R.H. and Stewart G.W.
          Algorithm 432: Solution of the matrix equation AX + XB = C.
          Comm. ACM., Vol. 15, pp. 820-826, 1972.

### Inputs

Name | Description
--- | ---
A[:, :] | Matrix A
B[:, :] | Matrix B
C[size(A, 1), size(B, 2)] | Matrix C
aIsSchur | True if A has already real Schur form
bIsSchur | True if B has already real Schur form

### Outputs

Name | Description
--- | ---
X[size(A, 1), size(B, 2)] | Solution of Sylvester equation

## Modelica_LinearSystems2.Math.Matrices.dsylvester

Solution of discrete-time Sylvester equation A*X*B + sgn*X = C

### Information

#### Syntax

>              X = Matrices.dsylvester(A, B, C);
>              X = Matrices.dsylvester(A, B, C, AisHess, BTisSchur, sgn, eps);

#### Description

Function dsylvester computes the solution X of the discrete-time Sylvester equation

>      A*X*B + sgn*X = C.

where sgn = 1 or sgn = -1. The algorithm applies the Hessenberg-Schur method proposed by Golub et al [1]. For sgn = -1, the discrete Sylvester equation is also known as Stein equation:

>      A*X*B - X + Q = 0.

In a nutshell, the problem is reduced to the corresponding problem

>      H*Y*S' + sgn*Y = F.

with H=U'*A*U is the Hessenberg form of A and S=V'*B'*V is the real Schur form of B', F=U'*C*V and Y=U*X*V' are appropriate transformations of C and X. This problem is solved sequently by exploiting the specific forms of S and H. Finally the solution of the the original problem is recovered as X=U'*Y*V.
The boolean inputs "AisHess" and "BTisSchur" indicate to omit one or both of the transformation to Hessenberg form or Schur form repectively in the case that A and/or B have already Hessenberg form or Schur respectively.

#### References

      [1] Golub, G.H., Nash, S. and Van Loan, C.F.
          A Hessenberg-Schur method for the problem AX + XB = C.
          IEEE Transaction on Automatic Control, AC-24, no. 6, pp. 909-913, 1979.

#### Example

>       A = [1.0,   2.0,   3.0;
>            6.0,   7.0,   8.0;
>            9.0,   2.0,   3.0];

>       B = [7.0,   2.0,   3.0;
>            2.0,   1.0,   2.0;
>            3.0,   4.0,   1.0];

>       C = [271.0,   135.0,   147.0;
>            923.0,   494.0,   482.0;
>            578.0,   383.0,   287.0];

>       X = discreteSylvester(A, B, C);

>       results in:
>       X = [2.0,   3.0,   6.0;
>            4.0,   7.0,   1.0;
>            5.0,   3.0,   2.0];

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] | Square matrix A in A*X*B + sgn*X = C
B[:, size(B, 1)] | Square matrix B in A*X*B + sgn*X = C
C[size(A, 2), size(B, 1)] | Rectangular matrix C in A*X*B + sgn*X = C
AisHess | True if A has already Hessenberg form
BTisSchur | True if B' has already real Schur form
sgn | Specifies the sign in A*X*B + sgn*X = C
eps | Tolerance

### Outputs

Name | Description
--- | ---
X[size(A, 2), size(B, 1)] | solution of the discrete Sylvester equation A*X*B + sgn*X = C

## Modelica_LinearSystems2.Math.Matrices.lyapunov

Solution of continuous-time Lyapunov equation X*A + A'*X = C

### Information

Function laypunov computes the solution X of the continuous-time Lyapunov equation

>      XA + A'*X = C.

using the Schur method for Lyapunov equations proposed by Bartels and Stewart [1].

References

      [1] Bartels, R.H. and Stewart G.W.
          Algorithm 432: Solution of the matrix equation AX + XB = C.
          Comm. ACM., Vol. 15, pp. 820-826, 1972.

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] |
C[size(A, 1), size(A, 2)] |
eps |

### Outputs

Name | Description
--- | ---
X[size(A, 1), size(A, 2)] | solution of the Lyapunov equation

## Modelica_LinearSystems2.Math.Matrices.dlyapunov

Solution of continuous-time Lyapunov equation A'X*A - X = C

### Information

Function laypunov computes the solution X of the continuous-time Lyapunov equation

>      XA + A'*X = C.

using the Schur method for Lyapunov equations proposed by Bartels and Stewart [1].

References

      [1] Bartels, R.H. and Stewart G.W.
          Algorithm 432: Solution of the matrix equation AX + XB = C.
          Comm. ACM., Vol. 15, pp. 820-826, 1972.

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] |
C[size(A, 1), size(A, 2)] |
eps |

### Outputs

Name | Description
--- | ---
X[size(A, 1), size(A, 2)] | solution of the Lyapunov equation

## Modelica_LinearSystems2.Math.Matrices.conditionNumber

Calculate the condition number norm(A)*norm(inv(A))

### Information

#### Syntax

>     r = Matrices.conditionNumber(A);

#### Description

This function calculates the the condition number (norm(A) * norm(inv(A))) of a general real matrix A, in either the 1-norm, 2-norm or the infinity-norm. In the case of 2-norm the result is he ratio of the largest to the smallest singular value to A.

#### Example

>       A = [1, 2
>            2, 1];
>       r = conditionNumber(A);

>       results in:

>       r = 3.0

#### See also

Matrices.rcond

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, :] | Input matrix
p | Type of p-norm (only allowed: 1, 2 or Modelica.Constants.inf)

### Outputs

Name | Description
--- | ---
result | p-norm of matrix A

## Modelica_LinearSystems2.Math.Matrices.rcond

Reciprocal condition number

### Information

#### Syntax

>     r = Matrices.rcond(A);
>     (r, info) = Matrices.rcond(A, false);

#### Description

This function estimates the reciprocal of the condition number (norm(A) * norm(inv(A))) of a general real matrix A, in either the 1-norm or the infinity-norm, using the LAPACK function DGECON.

#### Example

>       A = [1, 2
>            2, 1];
>       r = rcond(A);

>       results in:

>       r = 0.3333

#### See also

Matrices.conditionNumber

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] | Square real matrix
inf | Is true if infinity norm is used and false for 1-norm

### Outputs

Name | Description
--- | ---
rcond | Reciprocal condition number of A
info | Information

## Modelica_LinearSystems2.Math.Matrices.det

Determinant of a matrix (computed by LU decomposition)

### Information

#### Syntax

>     Matrices.det(A);

#### Description

This function call returns the determinant of matrix A computed by a LU decomposition. Usally, this function should never be used, because there are nearly always better numerical algorithms as by computing the determinant. E.g., use function Matrices.rank to compute the rank of a matrix.

#### See also

Matrices.rank, Matrices.solve

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] |

### Outputs

Name | Description
--- | ---
result | Determinant of matrix A

## Modelica_LinearSystems2.Math.Matrices.norm

Returns the norm of a matrix

### Information

#### Syntax

>     Matrices.norm(A);
>     Matrices.norm(A, p=2);

#### Description

The function call "`Matrices.norm(A)`" returns the 2-norm of matrix A, i.e., the largest singular value of A.
The function call "`Matrices.norm(A, p)`" returns the p-norm of matrix A. The only allowed values for p are

- "p=1": the largest column sum of A
- "p=2": the largest singular value of A
- "p=Modelica.Constants.inf": the largest row sum of A

Note, for any matrices A1, A2 the following inequality holds:

>     Matrices.norm(A1+A2,p) ≤ Matrices.norm(A1,p) + Matrices.norm(A2,p)

Note, for any matrix A and vector v the following inequality holds:

>     Vectors.norm(A*v,p) ≤ Matrices.norm(A,p)*Vectors.norm(A,p)

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, :] | Input matrix
p | Type of p-norm (only allowed: 1, 2 or Modelica.Constants.inf)

### Outputs

Name | Description
--- | ---
result | p-norm of matrix A

## Modelica_LinearSystems2.Math.Matrices.nullspace

Orthonormal nullspace of a matrix

### Information

#### Syntax

>     Z = Matrices.nullspace(A);

>     (Z, nullity) = Matrices.nullspace(A);

#### Description

This function calculates a orthonormal basis Z=[z_1, z_2, ...] of the nullspace of a matrix A, i.e. A*z_i=0.

The nullspace is obtained by svd method. That is, matrix A is decomposed into the matrices S, U, V:

>             T
>      A = USV

with the orthonormal matrices U and V and the matrix S with

>      S = [S1, 0]
>      S1 = [diag(s); 0]

with the singular values s={s1, s2, ..., sr} of A and r=rank(A). Note, that S has the same size as A. Since, U and V are orthonormal, we may write

>       T
>      U*A*V = [S1, 0].

Matrix S1 obviously has full column rank and therefore, the left n-r rows (n is the number of columns of A or S) of matrix V span a nullspace of A.

The nullity of matrix A is the dimension of the nullspace of A. In view of the above, it becomes clear that nullity holds

>      nullity = n - r

with

>      n = number of columns of matrix A
>      r = rank(A)

#### Example

>       A = [1, 2,  3, 1;
>            3, 4,  5, 2;
>           -1, 2, -3, 3];
>       (Z, nullity) = nullspace(A);

>       results in:

>       Z=[0.1715;
>         -0.686;
>          0.1715;
>          0.686]

>       nullity = 1

#### See also

Matrices.singularValues

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, :] | input matrix

### Outputs

Name | Description
--- | ---
Z[size(A, 2), :] | orthonormal nullspace of matrix A
nullity | nullity, i.e. the dimension of the nullspace

## Modelica_LinearSystems2.Math.Matrices.eigenValues

Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix

### Information

#### Syntax

>                     eigenvalues = Matrices.eigenValues(A);
>     (eigenvalues, eigenvectors) = Matrices.eigenValues(A);

#### Description

This function call returns the eigenvalues and optionally the (right) eigenvectors of a square matrix A. The first column of "eigenvalues" contains the real and the second column contains the imaginary part of the eigenvalues. If the i-th eigenvalue has no imaginary part, then eigenvectors[:,i] is the corresponding real eigenvector. If the i-th eigenvalue has an imaginary part, then eigenvalues[i+1,:] is the conjugate complex eigenvalue and eigenvectors[:,i] is the real and eigenvectors[:,i+1] is the imaginary part of the eigenvector of the i-th eigenvalue. With function Matrices.eigenValueMatrix, a real block diagonal matrix is constructed from the eigenvalues such that

>     A = eigenvectors * eigenValueMatrix(eigenvalues) * inv(eigenvectors)

provided the eigenvector matrix "eigenvectors" can be inverted (an inversion is possible, if all eigenvalues are different and no eigenvalue is zero).

#### Example

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];
>       Real eval;
>     algorithm
>       eval := Matrices.eigenValues(A);  // eval = [-0.618, 0;
>                                         //          8.0  , 0;
>                                         //          1.618, 0];

i.e., matrix A has the 3 real eigenvalues -0.618, 8, 1.618.

#### See also

Matrices.eigenValueMatrix, Matrices.singularValues

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] | Matrix

### Outputs

Name | Description
--- | ---
eigenvalues[size(A, 1), 2] | Eigenvalues of matrix A (Re: first column, Im: second column)
leftEigenvectors[size(A, 1), size(A, 2)] | Real-valued eigenvector matrix
rightEigenvectors[size(A, 1), size(A, 2)] | Real-valued eigenvector matrix

## Modelica_LinearSystems2.Math.Matrices.equalityLeastSquares

Solve a linear equality constrained least squares problem

### Information

#### Syntax

>     x = Matrices.equalityLeastSquares(A,a,B,b);

#### Description

This function returns the solution x of the linear equality-constrained least squares problem:

> min|A*x - a|^2 over x, subject to B*x = b

It is required that the dimensions of A and B fulfill the following relationship:

> size(B,1) ≤ size(A,2) ≤ size(A,1) + size(B,1)

Note, the solution is computed with the LAPACK function "dgglse" using the generalized RQ factorization under the assumptions that B has full row rank (= size(B,1)) and the matrix [A;B] has full column rank (= size(A,2)). In this case, the problem has a unique solution.

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, :] | Minimize |A*x - a|^2
a[size(A, 1)] |
B[:, size(A, 2)] | subject to B*x=b
b[size(B, 1)] |

### Outputs

Name | Description
--- | ---
x[size(A, 2)] | solution vector

## Modelica_LinearSystems2.Math.Matrices.solve

Solve real system of linear equations A*x=b with a b vector (Gaussian elemination with partial pivoting)

### Information

#### Syntax

>     Matrices.solve(A,b);

#### Description

This function call returns the solution x of the linear system of equations

> A*x = b

If a unique solution x does not exist (since A is singular), an exception is raised.

Note, the solution is computed with the LAPACK function "dgesv", i.e., by Gaussian elemination with partial pivoting.

#### Example

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];
>       Real b[3] = {10,22,12};
>       Real x[3];
>     algorithm
>       x := Matrices.solve(A,b);  // x = {3,2,1}

#### See also

Matrices.LU, Matrices.LU_solve

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] | Matrix A of A*x = b
b[size(A, 1)] | Vector b of A*x = b

### Outputs

Name | Description
--- | ---
x[size(b, 1)] | Vector x such that A*x = b

## Modelica_LinearSystems2.Math.Matrices.solve2

Solve real system of linear equations A*X=B with a B matrix (Gaussian elemination with partial pivoting)

### Information

#### Syntax

>     Matrices.solve2(A,b);

#### Description

This function call returns the solution X of the linear system of equations

> A*X = B

If a unique solution X does not exist (since A is singular), an exception is raised.

Note, the solution is computed with the LAPACK function "dgesv", i.e., by Gaussian elemination with partial pivoting.

#### Example

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];
>       Real B[3,2] = [10, 20;
>                      22, 44;
>                      12, 24];
>       Real X[3,2];
>     algorithm
>       (LU, pivots) := Matrices.LU(A);
>       X := Matrices.solve2(A, B1);  /* X = [3, 6;
>                                             2, 4;
>                                             1, 2] */

#### See also

Matrices.LU, Matrices.LU_solve2

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] | Matrix A of A*X = B
B[size(A, 1), :] | Matrix B of A*X = B

### Outputs

Name | Description
--- | ---
X[size(B, 1), size(B, 2)] | Matrix X such that A*X = B

## Modelica_LinearSystems2.Math.Matrices.solve2r

Solve real system of linear equations X*op(A)=B with a B matrix (Gaussian elemination with partial pivoting)

### Information

#### Syntax

>       X = Matrices.solve2r(A,B);
>       X = Matrices.solve2r(A, B, transA=false, isTriangular=false);

#### Description

This function call returns the solution X of the linear system of equations

> X*op(A) = B

with

> op(A) = transpose((A)) if transA==true op(A) = (A) if transA==false

If matrix (A) is already lower triangular, the factorization is avoided if input "isTriangular" is set true. If a unique solution X does not exist (since A is singular), an exception is raised.

Note, the solution is computed with the LAPACK function "dgesv", i.e., by Gaussian elemination with partial pivoting.

#### Example

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];

>       Real B[2,3]  = [10, 22, 12;
>                       20, 44, 24];
>       Real X[2,3];
>     algorithm
>       X := Matrices.solve2r(A, B);  /* X = [-34.0, 17.2, 2.4;
>                                             -68.0, 34.4, 4.8] */

#### See also

Matrices.LU, Matrices.LU_solve2

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] | Matrix A of X*op(A) = B
B[:, size(A, 1)] | Matrix B of X*op(A) = B
transA | True if op(A)=A', false if op(A)=A
isTriangular | True if the A is already lower triangular

### Outputs

Name | Description
--- | ---
X[size(B, 1), size(B, 2)] | Matrix X such that X*op(A) = B

## Modelica_LinearSystems2.Math.Matrices.LU

LU decomposition of square or rectangular matrix

### Information

#### Syntax

>     (LU, pivots)       = Matrices.LU(A);
>     (LU, pivots, info) = Matrices.LU(A);

#### Description

This function call returns the LU decomposition of a "Real[m,n]" matrix A, i.e.,

> P*L*U = A

where P is a permutation matrix (implicitely defined by vector `pivots`), L is a lower triangular matrix with unit diagonal elements (lower trapezoidal if m > n), and U is an upper triangular matrix (upper trapezoidal if m < n). Matrices L and U are stored in the returned matrix `LU` (the diagonal of L is not stored). With the companion function Matrices.LU_solve, this decomposition can be used to solve linear systems (P*L*U)*x = b with different right hand side vectors b. If a linear system of equations with just one right hand side vector b shall be solved, it is more convenient to just use the function Matrices.solve.

The optional third (Integer) output argument has the following meaning:

info = 0:successful exit
info > 0: | if info = i, U[i,i] is exactly zero. The factorization has been completed, but the factor U is exactly singular, and division by zero will occur if it is used to solve a system of equations.

The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].

#### Example

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];
>       Real b1[3] = {10,22,12};
>       Real b2[3] = { 7,13,10};
>       Real    LU[3,3];
>       Integer pivots[3];
>       Real    x1[3];
>       Real    x2[3];
>     algorithm
>       (LU, pivots) := Matrices.LU(A);
>       x1 := Matrices.LU_solve(LU, pivots, b1);  // x1 = {3,2,1}
>       x2 := Matrices.LU_solve(LU, pivots, b2);  // x2 = {1,0,2}

#### See also

Matrices.LU_solve, Matrices.solve,

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, :] | Square or rectangular matrix

### Outputs

Name | Description
--- | ---
LU[size(A, 1), size(A, 2)] | L,U factors (used with LU_solve(..))
pivots[min(size(A, 1), size(A, 2))] | pivot indices (used with LU_solve(..))
info | Information

## Modelica_LinearSystems2.Math.Matrices.LU_solve

Solve real system of linear equations P*L*U*x=b with a b vector and an LU decomposition (from LU(..))

### Information

#### Syntax

>     Matrices.LU_solve(LU, pivots, b);

#### Description

This function call returns the solution x of the linear systems of equations

> P*L*U*x = b;

where P is a permutation matrix (implicitely defined by vector `pivots`), L is a lower triangular matrix with unit diagonal elements (lower trapezoidal if m > n), and U is an upper triangular matrix (upper trapezoidal if m < n). The matrices of this decomposition are computed with function Matrices.LU that returns arguments `LU` and `pivots` used as input arguments of `Matrices.LU_solve`. With `Matrices.LU` and `Matrices.LU_solve` it is possible to efficiently solve linear systems with different right hand side vectors. If a linear system of equations with just one right hand side vector shall be solved, it is more convenient to just use the function Matrices.solve.

If a unique solution x does not exist (since the LU decomposition is singular), an exception is raised.

The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].

#### Example

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];
>       Real b1[3] = {10,22,12};
>       Real b2[3] = { 7,13,10};
>       Real    LU[3,3];
>       Integer pivots[3];
>       Real    x1[3];
>       Real    x2[3];
>     algorithm
>       (LU, pivots) := Matrices.LU(A);
>       x1 := Matrices.LU_solve(LU, pivots, b1);  // x1 = {3,2,1}
>       x2 := Matrices.LU_solve(LU, pivots, b2);  // x2 = {1,0,2}

#### See also

Matrices.LU, Matrices.solve,

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
LU[:, size(LU, 1)] | L,U factors of Matrices.LU(..) for a square matrix
pivots[size(LU, 1)] | Pivots indices of Matrices.LU(..)
b[size(LU, 1)] | Right hand side vector of P*L*U*x=b

### Outputs

Name | Description
--- | ---
x[size(b, 1)] | Solution vector such that P*L*U*x = b

## Modelica_LinearSystems2.Math.Matrices.LU_solve2

Solve real system of linear equations P*L*U*X=B with a B vector and an LU decomposition (from LU(..))

### Information

#### Syntax

>     Matrices.LU_solve(LU, pivots, B);

#### Description

This function call returns the solution x of the linear systems of equations

> P*L*U*X = B;

where P is a permutation matrix (implicitely defined by vector `pivots`), L is a lower triangular matrix with unit diagonal elements (lower trapezoidal if m > n), and U is an upper triangular matrix (upper trapezoidal if m < n). The matrices of this decomposition are computed with function Matrices.LU that returns arguments `LU` and `pivots` used as input arguments of `Matrices.LU_solve`. With `Matrices.LU` and `Matrices.LU_solve` it is possible to efficiently solve linear systems with different right hand side matrices. If a linear system of equations with just one right hand side matrix shall be solved, it is more convenient to just use the function Matrices.solve2.

If a unique solution X does not exist (since the LU decomposition is singular), an exception is raised.

The LU factorization is computed with the LAPACK function "dgetrf", i.e., by Gaussian elemination using partial pivoting with row interchanges. Vector "pivots" are the pivot indices, i.e., for 1 ≤ i ≤ min(m,n), row i of matrix A was interchanged with row pivots[i].

#### Example

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];
>       Real B1[3] = [10, 20;
>                     22, 44;
>                     12, 24];
>       Real B2[3] = [ 7, 14;
>                     13, 26;
>                     10, 20];
>       Real    LU[3,3];
>       Integer pivots[3];
>       Real    X1[3,2];
>       Real    X2[3,2];
>     algorithm
>       (LU, pivots) := Matrices.LU(A);
>       X1 := Matrices.LU_solve2(LU, pivots, B1);  /* X1 = [3, 6;
>                                                           2, 4;
>                                                           1, 2] */
>       X2 := Matrices.LU_solve2(LU, pivots, B2);  /* X2 = [1, 2;
>                                                           0, 0;
>                                                           2, 4] */

#### See also

Matrices.LU, Matrices.solve2,

### Inputs

Name | Description
--- | ---
LU[:, size(LU, 1)] | L,U factors of Matrices.LU(..) for a square matrix
pivots[size(LU, 1)] | Pivots indices of Matrices.LU(..)
B[size(LU, 1), :] | Right hand side matrix of P*L*U*X=B

### Outputs

Name | Description
--- | ---
X[size(B, 1), size(B, 2)] | Solution matrix such that P*L*U*X = B

## Modelica_LinearSystems2.Math.Matrices.leastSquares

Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient)

### Information

#### Syntax

>     x = Matrices.leastSquares(A,b);

#### Description

A linear system of equations A*x = b has no solutions or infinitely many solutions if A is not square. Function "leastSquares" returns a solution in a least squarse sense:

      size(A,1) > size(A,2):  returns x such that |A*x - b|^2 is a minimum
      size(A,1) = size(A,2):  returns x such that A*x = b
      size(A,1) < size(A,2):  returns x such that |x|^2 is a minimum for all
                              vectors x that fulfill A*x = b

Note, the solution is computed with the LAPACK function "dgelsx", i.e., QR or LQ factorization of A with column pivoting. If A does not have full rank, the solution is not unique and from the infinitely many solutions the one is selected that minimizes both |x|^2 and |A*x - b|^2.

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, :] | Matrix A
b[size(A, 1)] | Vector b

### Outputs

Name | Description
--- | ---
x[size(A, 2)] | Vector x such that min|A*x-b|^2 if size(A,1) >= size(A,2) or min|x|^2 and A*x=b, if size(A,1) < size(A,2)

## Modelica_LinearSystems2.Math.Matrices.leastSquares2

Solve overdetermined or underdetermined real system of linear equations A*X=B in a least squares sense (A may be rank deficient)

### Information

#### Syntax

>     X = Matrices.leastSquares2(A,B);

#### Description

Returns a solution of equation A*X = B in a least square sense (A may be rank deficient):

      minimize | A*X - B |

Several different cases can be distinguished (note, rank is an output argument of this function):

size(A,1) = size(A,2)

A solution is returned for a regular, as well as a singular matrix A:

- rank = size(A,1):
  A is regular and the returned solution X fulfills the equation A*X = B uniquely.
- rank < size(A,1):
  A is singular and no unique solution for equation A*X = B exists.
  - If an infinite number of solutions exists, the one is selected that fulfills the equation and at the same time has the minimum norm |x| for all solution vectors that fulfill the equation.
  - If no solution exists, X is selected such that |A*X - B| is as small as possible (but A*X - B is not zero).

size(A,1) > size(A,2):

The equation A*X = B has no unique solution. The solution X is selected such that |A*X - B| is as small as possible. If rank = size(A,2), this minimum norm solution is unique. If rank < size(A,2), there are an infinite number of solutions leading to the same minimum value of |A*X - B|. From these infinite number of solutions, the one with the minimum norm |X| is selected. This gives a unique solution that minimizes both |A*X - B| and |X|.

size(A,1) < size(A,2):

- rank = size(A,1):
  There are an infinite number of solutions that fulfill the equation A*X = B. From this infinite number, the unique solution is selected that minimizes |X|.
- rank < size(A,1):
  There is either no solution of equation A*X = B, or there are again an infinite number of solutions. The unique solution X is returned that minimizes both |A*X - B| and |X|.

Note, the solution is computed with the LAPACK function "dgelsx", i.e., QR or LQ factorization of A with column pivoting.

#### Algorithmic details

The function first computes a QR factorization with column pivoting:

          A * P = Q * [ R11 R12 ]
                      [  0  R22 ]

with R11 defined as the largest leading submatrix whose estimated condition number is less than 1/rcond. The order of R11, rank, is the effective rank of A.

Then, R22 is considered to be negligible, and R12 is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization:

         A * P = Q * [ T11 0 ] * Z
                     [  0  0 ]

The minimum-norm solution is then

         X = P * Z' [ inv(T11)*Q1'*B ]
                    [        0       ]

where Q1 consists of the first "rank" columns of Q.

#### See also

Matrices.leastSquares (same as leastSquares2, but with a right hand side vector),
Matrices.solve2 (for square, regular matrices A)

Extends from Modelica.Icons.Function (Icon for functions).

### Inputs

Name | Description
--- | ---
A[:, :] | Matrix A
B[size(A, 1), :] | Matrix B
rcond | Reciprocal condition number to estimate rank of A

### Outputs

Name | Description
--- | ---
X[size(A, 2), size(B, 2)] | Matrix X such that min|A*X-B|^2 if size(A,1) >= size(A,2) or min|X|^2 and A*X=B, if size(A,1) < size(A,2)
rank | Rank of A

## Modelica_LinearSystems2.Math.Matrices.QR

QR decomposition of a rectangular matrix without column pivoting (A = Q*R). Return the full square Q-matrix

### Inputs

Name | Description
--- | ---
A[:, :] | Rectangular matrix with size(A,1) >= size(A,2)

### Outputs

Name | Description
--- | ---
Q[size(A, 1), size(A, 2)] | Rectangular matrix with orthonormal columns such that Q*R=A[:,p]
R[min(size(A, 1), size(A, 2)), size(A, 2)] | Square upper triangular matrix
tau[min(size(A, 1), size(A, 2))] |
Q2[size(A, 1), size(A, 1)] |

## Modelica_LinearSystems2.Math.Matrices.LQ

LQ decomposition of a rectangular matrix without column pivoting (A = L*Q)

### Information

#### Syntax

>     (Q,R,p) = Matrices.QR(A);

#### Description

This function returns the QR decomposition of a rectangular matrix A (the number of columns of A must be less than or equal to the number of rows):

> Q*R = A[:,p]

where Q is a rectangular matrix that has orthonormal columns and has the same size as A (QTQ=I), R is a square, upper triangular matrix and p is a permutation vector. Matrix R has the following important properties:

- The absolute value of a diagonal element of R is the largest value in this row, i.e., abs(R[i,i]) ≥ abs(R[i,j]).
- The diagonal elements of R are sorted according to size, such that the largest absolute value is abs(R[1,1]) and abs(R[i,i]) ≥ abs(R[j,j]) with i < j.

This means that if abs(R[i,i]) ≤ ε then abs(R[j,k]) ≤ ε for j ≥ i, i.e., the i-th row up to the last row of R have small elements and can be treated as being zero. This allows to, e.g., estimate the row-rank of R (which is the same row-rank as A). Furthermore, R can be partitioned in two parts

>        A[:,p] = Q * [R1, R2;
>                      0,  0]

where R1 is a regular, upper triangular matrix.

Note, the solution is computed with the LAPACK functions "dgeqp3" and "dorgqr", i.e., by Housholder transformations with column pivoting. If Q is not needed, the function may be called as: `(,R,p) = QR(A)`.

#### Example

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];
>       Real R[3,3];
>     algorithm
>       (,R) := Matrices.QR(A);  // R = [-7.07.., -4.24.., -3.67..;
>                                         0     , -1.73.., -0.23..;
>                                         0     ,  0     ,  0.65..];

### Inputs

Name | Description
--- | ---
A[:, :] | Rectangular matrix

### Outputs

Name | Description
--- | ---
L[size(A, 1), size(A, 1)] | Rectangular matrix containing the lower triangular matrix
Q[size(A, 1), size(A, 2)] | Rectangular matrix with orthonormal columns such that L*Q=A

## Modelica_LinearSystems2.Math.Matrices.orthogonalQ

generates a real orthogonal matrix Q which is defined as the product of IHI-ILO elementary reflectors of order N, as returned by DGEHRD

### Information

       Purpose
    /*  ======= */

    /*  DORGHR generates a real orthogonal matrix Q which is defined as the */
    /*  product of IHI-ILO elementary reflectors of order N, as returned by */
    /*  DGEHRD: */

    /*  Q = H(ilo) H(ilo+1) . . . H(ihi-1). */

    /*  Arguments */
    /*  ========= */

    /*  N       (input) INTEGER */
    /*          The order of the matrix Q. N >= 0. */

    /*  ILO     (input) INTEGER */
    /*  IHI     (input) INTEGER */
    /*          ILO and IHI must have the same values as in the previous call */
    /*          of DGEHRD. Q is equal to the unit matrix except in the */
    /*          submatrix Q(ilo+1:ihi,ilo+1:ihi). */
    /*          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. */

    /*  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N) */
    /*          On entry, the vectors which define the elementary reflectors, */
    /*          as returned by DGEHRD. */
    /*          On exit, the N-by-N orthogonal matrix Q. */

    /*  LDA     (input) INTEGER */
    /*          The leading dimension of the array A. LDA >= max(1,N). */

    /*  TAU     (input) DOUBLE PRECISION array, dimension (N-1) */
    /*          TAU(i) must contain the scalar factor of the elementary */
    /*          reflector H(i), as returned by DGEHRD. */

    /*  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK) */
    /*          On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */

    /*  LWORK   (input) INTEGER */
    /*          The dimension of the array WORK. LWORK >= IHI-ILO. */
    /*          For optimum performance LWORK >= (IHI-ILO)*NB, where NB is */
    /*          the optimal blocksize. */

    /*          If LWORK = -1, then a workspace query is assumed; the routine */
    /*          only calculates the optimal size of the WORK array, returns */
    /*          this value as the first entry of the WORK array, and no error */
    /*          message related to LWORK is issued by XERBLA. */

    /*  INFO    (output) INTEGER */
    /*          = 0:  successful exit */
    /*          < 0:  if INFO = -i, the i-th argument had an illegal value */

    /*  ===================================================================== */

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] |
tau[size(A, 1) - 1] | scalar factors of the elementary reflectors
ilo | lowest index where the original matrix had been Hessenbergform - ilo must have the same value as in the previous call of DGEHRD
ihi | highest index where the original matrix had been Hessenbergform - ihi must have the same value as in the previous call of DGEHRD

### Outputs

Name | Description
--- | ---
Q[size(A, 1), size(A, 2)] | Orthogonal matrix as a result of elementary reflectors
info |

## Modelica_LinearSystems2.Math.Matrices.cholesky

Compute the Cholesky factorization of a symmetric positive definte matrix

### Information

### Inputs

Name | Description
--- | ---
A[:, size(A, 1)] |
upper | True if the upper triangle of A is provided

### Outputs

Name | Description
--- | ---
H[size(A, 1), size(A, 2)] | Cholesky factor U or L for A = U'*U or A = L*L'

## Modelica_LinearSystems2.Math.Matrices.choleskyUpDate

Compute the cholesky factor Lu according to Au=Lu'*Lu=A + v*v' with A=L'*L

### Information

#### Syntax

>       Lud = Matrices.Utilities.choleskyUpDate(L, v);
>       Lud = Matrices.Utilities.choleskyUpDate(L, v, true);

#### Description

Function choleskyUpDate(L, v) computes the rank-1-updated Cholesky factorization Lud, with

>                  T          T     T      T
>     Aud = Lud*Lud = A +  v*v = L*L +  v*v

from the input L, i.e. the left (lower) Cholesky factor of the original matrix A.
The approach is a transformation H*[v, L]' = [0, Lud]' with orthonormal matrix H such, that

>                           T          T         T               T     T
>        [0, Lud] * [0, Lud] = [v, L]*H *H*[v, L] = [v, L]*[v, L] = v*v + A

and matrix Lud is lower (upper) triangular. The transformation is performed by n (order of A) Givens rotations. The following sequence illustrates the principle of stepwise transformation of matrix [v, L]'. "*" are arbitrary elements. For each step the changed elements are bold.

>     | v' |    | * * * * |       | 0 * * * |       | 0 0 * * |       | 0 0 0 * |       | 0 0 0 0 |
>     |    |    | * * * * |       | * * * * |       | * * * * |       | * * * * |       | * * * * |
>     | L' | =  | 0 * * * |  ->   | 0 * * * |  ->   | 0 * * * |  ->   | 0 * * * |  ->   | 0 * * * |
>     |    |    | 0 0 * * |       | 0 0 * * |       | 0 0 * * |       | 0 0 * * |       | 0 0 * * |
>     |    |    | 0 0 0 * |       | 0 0 0 * |       | 0 0 0 * |       | 0 0 0 * |       | 0 0 0 * |

With the boolean input "upper" the user specifies wether the matrix L is lower or upper triangular matrix (left or right Cholesky factor). If "upper==true", the output Lud is also upper triangular. Default is "upper==false".

#### Example

> ###

> ###

> ##

> ###

> ####

> ####

>     Inputs

>     NameDescription
>     L[:, size(L, 1)]Cholesky factor
>     v[size(L, 1)]Real vector A+v*v'
>     upperTrue if the upper triangle of A is provided and the modified upper triangle will be returned

>     Outputs

>     NameDescription
>     Lud[size(L, 1), size(L, 2)]Updated Cholesky factor

>     Modelica_LinearSystems2.Math.Matrices.choleskyDownDate
>     Compute the cholesky factor Ld according to Ad=Ld'*Ld=A - v*v' with A=L'*L
>     Information

>     Syntax

>       Ldd = Matrices.Utilities.choleskyUpDate(L, v);
>       Ldd = Matrices.Utilities.choleskyUpDate(L, v, true);

>     Description

>     Function choleskyDownDate(L, v) computes the rank-1-downdated Cholesky factorization Ldd, with

>                  T          T     T      T
>     Add = Ldd*Ldd = A -  v*v = L*L -  v*v

>     from the input L, i.e. the left (lower) Cholesky factor of the original matrix A.

>     The algortihm is taken from [1].

>     Matrix Ldd is calculated by

>               T           T
>       [v, Ldd] = H *[0, L]

> >                           T          T         T               T     T
> >        [v, Ldd] * [v, Ldd] = [0, L]*H *H*[0, L] = [0, L]*[0, L] = L*L + A,

> i.e., by orthogonal transformation

> >      H = H_1*...*H_n.

> The matrices H_i are Givens matrices computed such, that

> >                          T T                T
> >     H_1*H_2*...*H_n*[z, a ] = [1, 0, ..., 0],

> with a is the solution of

> >     L*a = x

> and

> >     z = ||a||.

> The following sequence illustrate the principle of calculating the H_i, starting with H_n

> >     |z|      |z|      |z|      |z|
> >     |a| H_3  |a| H_2  |a| H_1  |0|
> >     |a| -->  |a| -->  |0| -->  |0|
> >     |a|      |0|      |0|      |0|

> Note, that the z and a are different in each column.
> It is shown in [1] that this algorithms results in the modified Cholesky factor Ldd.
> With the boolean input "upper" the user specifies wether the matrix L is lower or upper triangular matrix (left or right Cholesky factor). If "upper==true", the output Ldd is also upper triangular. Default is "upper==false".

> #### References

> [1] | Dongarra, J. J., Bunch, J. R., Moler, G. B., Stewart, G.W. | "LINPACK Users' Guide " | Society for Industrial Mathematics, 1987

#### Example

> ###

> ###

> ##

> ###

> ####

> ####

> ####

> ###

> ###

> ##

> ###

> ####

> ####

> ####

> ###

> ###

> ##

> ###

> ####

> ####

> ####

> ####

> ###

> ###

> ##

> ###

> ####

> ####

> ####

> ###

> ###

> ##

> ###

> ####

> ####

> ####

> ####

> ###

> ###

> ##

> ###

> ####

> ####

> ####

> ####

> ###

> ###

> ##

> ###

> ###

> ##

> ###

> ###

> ##

> ###

> ###

> ##

> ###

> ###

> ###

> ##

> ###

> ###

> ###

> ##

> ###

> ###

> ##

> ###

> ###

>     Inputs

>     NameDescription
>     L[:, size(L, 1)]Cholesky factor
>     v[size(L, 1)]Real vector A+v*v'
>     upperTrue if the upper triangle of A is provided

>     Outputs

>     NameDescription
>     Ldd[size(L, 1), size(L, 2)]Updated Cholesky factor

>     Modelica_LinearSystems2.Math.Matrices.rsf
>     Computes the real Schur form (RSF) of a square matrix
>     Information

>     Syntax

>        (T, Z, alphaReal, alphaImag)  =     Modelica_LinearSystems2.Math.Matrices.(A)

>     Description

>     Function rsf (real Schur form) calculates the real Schur form af a real square matrix A, i.e.

>              T
>     A = Z*T*Z

>     with the real nxn matrices T and Z. Z is an orthogonal matrix.  T is an block upper triangular matrix with 1x1 and 2x2 blocks in the diagonal.The 1x1 blocks contains the real eigenvalues of a. The 2x2 blocks are matrices with the conjugated complex pairs of eigenvalues, whereas the real parts of the eigenvalues are the elements of the diagonal.

>     The calculation is performed stepwise using several lapack routines. First, lapack.dgehrd reduces matrix A is to upper Hessenberg form H=Q'AQ, whereas Q is computed by lapack.dodrghr.Finally, lapack.dhseqr transforms H to T. The eigenvalues of A are calculated straightforward from T.

>     Function rsf does not apply lapack.dgees, a routine to directly compute the real Schur from. See also
>     Math.Matrices.rsf2

>     Example

>        Real A[3,3] = [1, 2, 3; 4, 5, 6; 7, 8, 9];
>        Real T[3,3];
>        Real Z[3,3];
>        Real alphaReal[3];
>        Real alphaImag[3];

>     algorithm
>       (T, Z, alphaReal, alphaImag):=Modelica_LinearSystems2.Math.Matrices.rsf(A);
>     //   T = [16.12, 4.9,   1.59E-015;
>     //        0,    -1.12, -1.12E-015;
>     //        0,     0,    -1.30E-015]
>     //   Z = [-0.23,  -0.88,   0.41;
>     //        -0.52,  -0.24,  -0.82;
>     //        -0.82,   0.4,    0.41]
>     //alphaReal = {16.12, -1.12, -1.32E-015}
>     //alphaImag = {0, 0, 0}

>     Inputs

>     NameDescription
>     A[:, size(A, 1)]

>     Outputs

>     NameDescription
>     T[size(A, 1), size(A, 2)]
>     Z[size(A, 1), size(A, 2)]
>     alphaReal[size(A, 1)]Real part of eigenvalue=alphaReal+i*alphaImag
>     alphaImag[size(A, 1)]Imaginary part of eigenvalue=(alphaReal+i*alphaImag

>     Modelica_LinearSystems2.Math.Matrices.rsf2
>     Computes the real Schur form (RSF) of a square matrix but uses lapack.dgees
>     Information

>     Syntax

>        (T, Z, alphaReal, alphaImag)  =     Modelica_LinearSystems2.Math.Matrices.(A)

>     Description

>     Function rsf2 (real Schur form) calculates the real Schur form af a real square matrix A, i.e.

>              T
>     A = Z*T*Z

>     with the real nxn matrices T and Z. Z is an orthogonal matrix.  T is an block upper triangular matrix with 1x1 and 2x2 blocks in the diagonal.The 1x1 blocks contains the real eigenvalues of a. The 2x2 blocks are matrices with the conjugated complex pairs of eigenvalues, whereas the real parts of the eigenvalues are the elements of the diagonal.

>     The calculation is performed stepwise using lapack.dgees, i.e. using the internal mehtods of balacing and scaling of dgees.

>     See also Math.Matrices.rsf

>     Example

>        Real A[3,3] = [1, 2, 3; 4, 5, 6; 7, 8, 9];
>        Real T[3,3];
>        Real Z[3,3];
>        Real alphaReal[3];
>        Real alphaImag[3];

>     algorithm
>       (T, Z, alphaReal, alphaImag):=Modelica_LinearSystems2.Math.Matrices.rsf2(A);
>     //   T = [16.12, 4.9,   1.59E-015;
>     //        0,    -1.12, -1.12E-015;
>     //        0,     0,    -1.30E-015]
>     //   Z = [-0.23,  -0.88,   0.41;
>     //        -0.52,  -0.24,  -0.82;
>     //        -0.82,   0.4,    0.41]
>     //alphaReal = {16.12, -1.12, -1.32E-015}
>     //alphaImag = {0, 0, 0}

>     Inputs

>     NameDescription
>     A[:, size(A, 1)]

>     Outputs

>     NameDescription
>     S[size(A, 1), size(A, 2)]
>     QZ[size(A, 1), size(A, 2)]
>     alphaReal[size(A, 1)]Real part of eigenvalue=alphaReal+i*alphaImag
>     alphaImag[size(A, 1)]Imaginary part of eigenvalue=(alphaReal+i*alphaImag

>     Modelica_LinearSystems2.Math.Matrices.hessenberg
>     Transform a matrix to upper Hessenberg form
>     Information

>     Syntax

>              H = Matrices.hessenberg(A);
>         (H, U) = Matrices.hessenberg(A);

>     Description
>     Function hessenberg computes the Hessenberg matrix H of matrix A as well as the orthogonal transformation matrix U that holds H = U'*A*U.
>     The Hessenberg form of a matrix is computed by repeated Householder similarity transformation. The elementary reflectors and the corresponding scalar factors are provided
>     by function "Utilities.toUpperHessenberg()". The transformation matrix U is then computed by
>     LAPACK.dorghr.

>     Example

>      A  = [1, 2,  3;
>            6, 5,  4;
>            1, 0,  0];

>      (H, U) = hessenberg(A);

>       results in:

>      H = [1.0,  -2.466,  2.630;
>          -6.083, 5.514, -3.081;
>           0.0,   0.919, -0.514]

>      U = [1.0,    0.0,      0.0;
>           0.0,   -0.9864,  -0.1644;
>           0.0,   -0.1644,   0.9864]

>       and therefore,

>      u*H*transpose(U) = [1.0, 2.0, 3.0;
>                          6.0, 5.0, 4.0;
>                          1.0, 0.0, 0.0]

>     See also
>     Matrices.Utilities.toUpperHessenberg

>     Inputs

>     NameDescription
>     A[:, :]Square matrix A

>     Outputs

>     NameDescription
>     H[size(A, 1), size(A, 2)]Hessenberg form of A
>     U[size(A, 1), size(A, 2)]Transformation matrix

>     Modelica_LinearSystems2.Math.Matrices.toUpperHessenberg
>     transform a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation:  Q' * A * Q = H
>     Information

>        Syntax

>              H = Matrices.UtilitiestoUpperHessenberg(A);
>              (H, V, tau, info) = Matrices.UtilitiestoUpperHessenberg(A,ilo, ihi);

>     Description
>     Function toUpperHessenberg computes a upper Hessenberg form H of a matrix A by orthogonal similarity transformation:  Q' * A * Q = H.
>     It calls LAPACK function DGEHRD. See Matrices.Lapack.dgehrd for more information about the additional outputs V, tau, info and inputs ilo, ihi for more information.

>     Example

>      A  = [1, 2,  3;
>            6, 5,  4;
>            1, 0,  0];

>      H = toUpperHessenberg(A);

>       results in:

>      H = [1.0,  -2.466,  2.630;
>          -6.083, 5.514, -3.081;
>           0.0,   0.919, -0.514]

>     Inputs

>     NameDescription
>     A[:, size(A, 1)]Square matrix A
>     ilolowest index where the original matrix had been Hessenbergform
>     ihihighest index where the original matrix had been Hessenbergform

>     Outputs

>     NameDescription
>     H[size(A, 1), size(A, 2)]Upper Hessenberg form
>     V[size(A, 1), size(A, 2)]V=[v1,v2,..vn-1,0] with vi are vectors which define the elementary reflectors
>     tau[max(0, size(A, 1) - 1)]Scalar factors of the elementary reflectors
>     info

>     Modelica_LinearSystems2.Math.Matrices.householderReflexion
>     Reflect each of the vectors ai of matrix  A=[a1, a2, ..., an] on a plane with orthogonal vector u
>     Information

>     Syntax

>     Matrices.householderReflection(A,u);

>     Description

>     This function computes the Housholder reflection (transformation)

>      Ar = Q*A

>     with

>      Q = I -2*u*u'/(u'*u)

>     where u*u is housholder vector, i.e. the normal vector of the reflection plane.

>     Householder reflection is widely used in numerical linear algebra, e.g. to perform QR decompositions.

>     Example

>     // First step of QR decomposition
>       import   Modelica.Math.Vectors.Utilities;

>       Real A[3,3] = [1,2,3;
>                      3,4,5;
>                      2,1,4];
>       Real Ar[3,3];
>       Real u[:];

>       u=Utilities.householderVector(A[:,1],{1,0,0});
>       // u={ = {0.763, 0.646, 0}
>       Ar=householderReflexion(A,u);
>      // Ar = [-6.0828,   -5.2608,   -4.4388;
>      //        0.0,      -1.1508,   -2.3016;
>      //        0.0,       2.0,       0.0]

>     See also

>     Matrices.housholderSimilarityTransformation,

>     Inputs

>     NameDescription
>     A[:, :]Rectangular matrix
>     u[size(A, 1)]Householder vector

>     Outputs

>     NameDescription
>     RA[size(A, 1), size(A, 2)]Reflexion of A

>     Modelica_LinearSystems2.Math.Matrices.householderSimilarityTransformation
>     Calculate the similarity transformation S*A*S of matrix A with symmetric householder matrix S = I - 2u*u'
>     Information

>     Syntax

>     Matrices.householderSimilarityTransformation(A,u);

>     Description

>     This function computes the Housholder similarity transformation

>      As = S*A*S

>     with

>      S = I -2*u*u'/(u'*u).

>     This transformation is widely used for transforming non-symmetric matrices to a Hessenberg form.

>     Example

>     // First step of Hessenberg decomposition
>       import   Modelica.Math.Vectors.Utilities;

>       Real A[4,4] = [1,2,3,4;
>                      3,4,5,6;
>                      9,8,7,6;
>                      1,2,0,0];
>       Real Ar[4,4];
>       Real u[4]={0,0,0,0};

>       u[2:4]=Utilities.householderVector(A[2:4,1],{1,0,0});
>       // u= = {0, 0.8107, 0.5819, 0.0647}
>       Ar=householderSimilarityTransformation(A,u);
>      //  Ar = [1.0,     -3.8787,    -1.2193,    3.531;
>               -9.5394, 11.3407,      6.4336,   -5.9243;
>                0.0,     3.1307,      0.7525,   -3.3670;
>                0.0,     0.8021,     -1.1656,   -1.0932]

>     See also

>     Matrices.housholderReflection,

>     Inputs

>     NameDescription
>     A[:, size(A, 1)]Square matrix A
>     u[size(A, 1)]Householder vector

>     Outputs

>     NameDescription
>     SAS[size(A, 1), size(A, 1)]

>     Modelica_LinearSystems2.Math.Matrices.fliplr
>     flip the columns of a matrix in left/right direction
>     Inputs

>     NameDescription
>     A[:, :]Matrix to be fliped

>     Outputs

>     NameDescription
>     Aflip[size(A, 1), size(A, 2)]fliped matrix

>     Modelica_LinearSystems2.Math.Matrices.flipud
>     flip the columns of a matrix in up/down direction
>     Inputs

>     NameDescription
>     A[:, :]Matrix to be flipped

>     Outputs

>     NameDescription
>     Aflip[size(A, 1), size(A, 2)]flipped matrix

>     Modelica_LinearSystems2.Math.Matrices.triangle
>     Return the upper/lower triangular part of a square matrix
>     Inputs

>     NameDescription
>     A[:, size(A, 1)]Square matrix
>     upperTrue for upper triangle to return

>     Outputs

>     NameDescription
>     Tri[size(A, 1), size(A, 2)]Triangular matrix

>     Modelica_LinearSystems2.Math.Matrices.trace
>     tarce(A) is the sum of the diagonal elements of A
>     Information
>     Extends from Modelica.Icons.Function (Icon for functions).
>     Inputs

>     NameDescription
>     A[:, size(A, 1)]

>     Outputs

>     NameDescription
>     result

>     Modelica_LinearSystems2.Math.Matrices.generalizedEigenvaluesTriangular
>     Compute invariant zeros of linear state space system with a generalized system matrix [A, B, C, D] which is of upper Hessenberg form
>     Information
>     This function is an interface to LAPACK routine DHGEQZ to calculate invariant
>     zeros of systems with generalized system matrices of upper Hessenberg form.
>     DHGEQZ is described below:

>          Purpose
>        ==========================================================

>        DHGEQZ implements a single-/double-shift version of the QZ method for
>        finding the generalized eigenvalues

>        w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j)   of the equation

>             det( A - w(i) B ) = 0

>        In addition, the pair A,B may be reduced to generalized Schur form:
>        B is upper triangular, and A is block upper triangular, where the
>        diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having
>        complex generalized eigenvalues (see the description of the argument
>        JOB.)

>        If JOB='S', then the pair (A,B) is simultaneously reduced to Schur
>        form by applying one orthogonal tranformation (usually called Q) on
>        the left and another (usually called Z) on the right.  The 2-by-2
>        upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks
>        of A will be reduced to positive diagonal matrices.  (I.e.,
>        if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and
>        B(j+1,j+1) will be positive.)

>        If JOB='E', then at each iteration, the same transformations
>        are computed, but they are only applied to those parts of A and B
>        which are needed to compute ALPHAR, ALPHAI, and BETAR.

>        If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal
>        transformations used to reduce (A,B) are accumulated into the arrays
>        Q and Z s.t.:

>             Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
>             Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*

>        Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
>             Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
>             pp. 241--256.

>        Arguments
>        =========

>        JOB     (input) CHARACTER*1
>                = 'E': compute only ALPHAR, ALPHAI, and BETA.  A and B will
>                       not necessarily be put into generalized Schur form.
>                = 'S': put A and B into generalized Schur form, as well
>                       as computing ALPHAR, ALPHAI, and BETA.

>        COMPQ   (input) CHARACTER*1
>                = 'N': do not modify Q.
>                = 'V': multiply the array Q on the right by the transpose of
>                       the orthogonal tranformation that is applied to the
>                       left side of A and B to reduce them to Schur form.
>                = 'I': like COMPQ='V', except that Q will be initialized to
>                       the identity first.

>        COMPZ   (input) CHARACTER*1
>                = 'N': do not modify Z.
>                = 'V': multiply the array Z on the right by the orthogonal
>                       tranformation that is applied to the right side of
>                       A and B to reduce them to Schur form.
>                = 'I': like COMPZ='V', except that Z will be initialized to
>                       the identity first.

>        N       (input) INTEGER
>                The order of the matrices A, B, Q, and Z.  N >= 0.

>        ILO     (input) INTEGER
>        IHI     (input) INTEGER
>                It is assumed that A is already upper triangular in rows and
>                columns 1:ILO-1 and IHI+1:N.
>                1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.

>        A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
>                On entry, the N-by-N upper Hessenberg matrix A.  Elements
>                below the subdiagonal must be zero.
>                If JOB='S', then on exit A and B will have been
>                   simultaneously reduced to generalized Schur form.
>                If JOB='E', then on exit A will have been destroyed.
>                   The diagonal blocks will be correct, but the off-diagonal
>                   portion will be meaningless.

>        LDA     (input) INTEGER
>                The leading dimension of the array A.  LDA >= max( 1, N ).

>        B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
>                On entry, the N-by-N upper triangular matrix B.  Elements
>                below the diagonal must be zero.  2-by-2 blocks in B
>                corresponding to 2-by-2 blocks in A will be reduced to
>                positive diagonal form.  (I.e., if A(j+1,j) is non-zero,
>                then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be
>                positive.)
>                If JOB='S', then on exit A and B will have been
>                   simultaneously reduced to Schur form.
>                If JOB='E', then on exit B will have been destroyed.
>                   Elements corresponding to diagonal blocks of A will be
>                   correct, but the off-diagonal portion will be meaningless.

>        LDB     (input) INTEGER
>                The leading dimension of the array B.  LDB >= max( 1, N ).

>        ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
>                ALPHAR(1:N) will be set to real parts of the diagonal
>                elements of A that would result from reducing A and B to
>                Schur form and then further reducing them both to triangular
>                form using unitary transformations s.t. the diagonal of B
>                was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block
>                (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j).
>                Note that the (real or complex) values
>                (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
>                generalized eigenvalues of the matrix pencil A - wB.

>        ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
>                ALPHAI(1:N) will be set to imaginary parts of the diagonal
>                elements of A that would result from reducing A and B to
>                Schur form and then further reducing them both to triangular
>                form using unitary transformations s.t. the diagonal of B
>                was non-negative real.  Thus, if A(j,j) is in a 1-by-1 block
>                (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0.
>                Note that the (real or complex) values
>                (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
>                generalized eigenvalues of the matrix pencil A - wB.

>        BETA    (output) DOUBLE PRECISION array, dimension (N)
>                BETA(1:N) will be set to the (real) diagonal elements of B
>                that would result from reducing A and B to Schur form and
>                then further reducing them both to triangular form using
>                unitary transformations s.t. the diagonal of B was
>                non-negative real.  Thus, if A(j,j) is in a 1-by-1 block
>                (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j).
>                Note that the (real or complex) values
>                (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
>                generalized eigenvalues of the matrix pencil A - wB.
>                (Note that BETA(1:N) will always be non-negative, and no
>                BETAI is necessary.)

>        Q       (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
>                If COMPQ='N', then Q will not be referenced.
>                If COMPQ='V' or 'I', then the transpose of the orthogonal
>                   transformations which are applied to A and B on the left
>                   will be applied to the array Q on the right.

>        LDQ     (input) INTEGER
>                The leading dimension of the array Q.  LDQ >= 1.
>                If COMPQ='V' or 'I', then LDQ >= N.

>        Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
>                If COMPZ='N', then Z will not be referenced.
>                If COMPZ='V' or 'I', then the orthogonal transformations
>                   which are applied to A and B on the right will be applied
>                   to the array Z on the right.

>        LDZ     (input) INTEGER
>                The leading dimension of the array Z.  LDZ >= 1.
>                If COMPZ='V' or 'I', then LDZ >= N.

>        WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
>                On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.

>        LWORK   (input) INTEGER
>                The dimension of the array WORK.  LWORK >= max(1,N).

>                If LWORK = -1, then a workspace query is assumed; the routine
>                only calculates the optimal size of the WORK array, returns
>                this value as the first entry of the WORK array, and no error
>                message related to LWORK is issued by XERBLA.

>        INFO    (output) INTEGER
>                = 0: successful exit
>                < 0: if INFO = -i, the i-th argument had an illegal value
>                = 1,...,N: the QZ iteration did not converge.  (A,B) is not
>                           in Schur form, but ALPHAR(i), ALPHAI(i), and
>                           BETA(i), i=INFO+1,...,N should be correct.
>                = N+1,...,2*N: the shift calculation failed.  (A,B) is not
>                           in Schur form, but ALPHAR(i), ALPHAI(i), and
>                           BETA(i), i=INFO-N+1,...,N should be correct.
>                > 2*N:     various "impossible" errors.

>        Further Details
>        ===============

>        Iteration counters:

>        JITER  -- counts iterations.
>        IITER  -- counts iterations run since ILAST was last
>                  changed.  This is therefore reset only when a 1-by-1 or
>                  2-by-2 block deflates off the bottom.

>        =====================================================================

>     Inputs

>     NameDescription
>     A[:, size(A, 1)]
>     B[size(A, 1), size(A, 1)]

>     Outputs

>     NameDescription
>     alphaReal[size(A, 1)]Real part of alpha (eigenvalue=(alphaReal+i*alphaImag)/beta)
>     alphaImag[size(A, 1)]Imaginary part of alpha
>     beta[size(A, 1)]Denominator of eigenvalue
>     info

>     Modelica_LinearSystems2.Math.Matrices.printMatrix
>     print matrix
>     Inputs

>     NameDescription
>     M[:, :]
>     significantDigitsNumber of significant digits that are shown
>     nameIndependent variable name used for printing

>     Outputs

>     NameDescription
>     s

>     Modelica_LinearSystems2.Math.Matrices.fromFile
>     Read matrix from a matlab file
>     Inputs

>     NameDescription
>     fileName
>     matrixNameName of the matrix

>     Outputs

>     NameDescription
>     A[n, m]