## Modelica_LinearSystems2.Math.Matrices

Additional functions for Modelica.Math

### Information

Extends from Modelica.Icons.Library (Icon for library).

### Package Content

Name | Description
--- | ---
LAPACK |
Examples |
care | Solution of continuous-time algebraic Riccati equations
dare | Solution of discrete-time algebraic Riccati equations
conditionNumber | Calculate the condition number norm(A)*norm(inv(A))
det | Determinant of a matrix (computed by LU decomposition)
eigenValues | Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix
equalityLeastSquares | Solve a linear equality constrained least squares problem
fliplr | flip the columns of a matrix in left/right direction
flipud | flip the columns of a matrix in up/down direction
fromFile | Read matrix from a matlab file
generalizedEigenvaluesTriangular | Compute invariant zeros of linear state space system with a genralized system matrix [A, B, C, D] which is of upper Hessenberg form
hessenberg | Compute an upper Hessenberg matrix by repeatedly applicated householder similarity transformation
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 SAS of matrix A with householder matrix S = I - 2u*u'
leastSquares | Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient)
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(..))
lyapunov | Solution of continuous-time Lyapunov equation X*A + A'*X = C
norm | Returns the norm of a matrix
nullspace | Orthonormal nullspace of a matrix
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
printMatrix | print matrix
QR | QR decomposition of a square matrix without column pivoting (A = Q*R)
rcond | Reciprocal condition number
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
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)
sylvester | Solution of continuous-time Sylvester equation A*X + X*B = C
toUpperHessenberg | transform a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H
trace | tarce(A) is the sum of the diagonal elements of A
Internal |

## 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

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, size(A, 1)] |  |
Real | B[size(A, 1), :] |  |
Real | R[size(B, 2), size(B, 2)] | identity(size(B, 2)) |
Real | Q[size(A, 1), size(A, 1)] | identity(size(A, 1)) |
Boolean | refine | false |

### Outputs

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

### Modelica definition

    function care
      "Solution of continuous-time algebraic Riccati equations"
      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Complex;

      input Real A[:,size(A, 1)];
      input Real B[size(A, 1),:];
      input Real R[size(B, 2),size(B, 2)]=identity(size(B, 2));
      input Real Q[size(A, 1),size(A, 1)]=identity(size(A, 1));
      input Boolean refine = false;
    protected
      Integer n=size(A, 1);
      Real G[size(A, 1),size(A, 1)]=B*Modelica.Math.Matrices.solve2(R, transpose(B));
      Real H[:,:]=[A,-G; -Q,-transpose(A)];
      Real H_RSF[:,:]=H;
      Real Z[size(H, 1),size(H, 2)];
      Real Z11[size(A, 1),size(A, 2)];
      Real Z21[size(A, 1),size(A, 2)];
      Real alphaReal[size(H, 1)] "Real part of eigenvalue=alphaReal+i*alphaImag";
      Real alphaImag[size(H, 1)]
        "Imaginary part of eigenvalue=(alphaReal+i*alphaImag";
      Integer info;
      Integer evSize;
      Complex xc[2];
    public
      output Real X[size(A, 1),size(A, 2)] "stabilizing solution of CARE";
      output Complex ev[size(A, 1)] "eigenvalues of the closed loop system";

    algorithm
      if n > 1 then
        (H_RSF,Z,alphaReal,alphaImag) := Matrices.rsf(H);
        (H_RSF,Z,alphaReal,alphaImag) := Matrices.Internal.reorderRSF(
          true,
          H_RSF,
          Z,
          alphaReal,
          alphaImag);

        evSize := size(ev, 1);
        for i in 1:evSize loop
          ev[i] := Complex(alphaReal[i], alphaImag[i]);
        end for;

        Z11 := Z[1:n, 1:n];
        Z21 := Z[n + 1:2*n, 1:n];
        if size(Z11, 1) > 0 then
    //  X := transpose(Matrices.solve2(transpose(Z11), transpose(Z21)));
          (X,info) := Matrices.LAPACK.dgesvx(Z11, transpose(Z21));
          //this function does not need to transpose Z11 as solve2 does
          X := transpose(X);
          assert(info == 0,
            "Solving a linear system of equations with function \"Matrices.LAPACK.dgesvx\" is not possible, because the system has either no or infinitely many solutions (input A is singular).");
          if refine then
            X := Matrices.Internal.carenls(A, B, R, Q, X);
          end if;
        else
          X := fill(
            0,
            size(Z21, 1),
            size(Z11, 1));
        end if;

      elseif n == 1 then
        xc := Polynomial.roots(Polynomial({-G[1, 1],2*A[1, 1],Q[1, 1]}));
        X := matrix(-abs(xc[1].re));
      else
        X := fill(0, 0, 0);
      end if;

    end care;

## 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 = 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*B)  *B'*X*A

If U is partitioned to

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

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

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, size(A, 1)] |  |
Real | B[size(A, 1), :] |  |
Real | R[size(B, 2), size(B, 2)] | identity(size(B, 2)) |
Real | Q[size(A, 1), size(A, 1)] | identity(size(A, 1)) |

### Outputs

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

### Modelica definition

    function dare "Solution of discrete-time algebraic Riccati equations"
      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Complex;

      input Real A[:,size(A, 1)];
      input Real B[size(A, 1),:];
      input Real R[size(B, 2),size(B, 2)]=identity(size(B, 2));
      input Real Q[size(A, 1),size(A, 1)]=identity(size(A, 1));
    protected
      Integer n=size(A, 1);
      Real G[size(A, 1),size(A, 1)]=B*Modelica.Math.Matrices.solve2(R, transpose(B));
      Real invAT[:,:]=transpose(Modelica.Math.Matrices.inv(A));
      Real H[:,:]=[A + G*invAT*Q,-G*invAT; -invAT*Q,invAT];
      Real H_RSF[:,:]=H;
      Real Z[size(H, 1),size(H, 2)];
      Real Z11[size(A, 1),size(A, 2)];
      Real Z21[size(A, 1),size(A, 2)];
      Real alphaReal[size(H, 1)] "Real part of eigenvalue=alphaReal+i*alphaImag";
      Real alphaImag[size(H, 1)]
        "Imaginary part of eigenvalue=(alphaReal+i*alphaImag";

      Integer info;
      Integer evSize;
      Complex j = Modelica_LinearSystems2.Math.Complex.j();
    public
      output Real S[size(A, 1),size(A, 2)]
        "orthogonal matrix of the Schur vectors associated to ordered rsf";
      output Complex ev[size(A, 1)] "eigenvalues of the closed loop system";

    algorithm
      (H_RSF,Z,alphaReal,alphaImag) := Matrices.rsf(H);// put H to Schur form
      (H_RSF,Z,alphaReal,alphaImag) := Matrices.Internal.reorderRSF(
        false,
        H_RSF,
        Z,
        alphaReal,
        alphaImag);// ordered Schur form
        evSize := size(ev, 1);
      for i in 1:evSize loop
        ev[i] := alphaReal[i] + j*alphaImag[i];
      end for;

      Z11 := Z[1:n, 1:n];
      Z21 := Z[n + 1:2*n, 1:n];
      if size(Z11, 1) > 0 then
    //  S := transpose(Matrices.solve2(transpose(Z11), transpose(Z21)));
        (S,info) := Matrices.LAPACK.dgesvx(Z11, transpose(Z21));//function does not need to transpose Z11 as solve2 does
        S := transpose(S);
        assert(info == 0, "Solving a linear system of equations with function
    \"Matrices.LAPACK.dgesvx\" is not possible, because the system has either
    no or infinitely many solutions (input A is singular).");
      else
        S := fill(
          0,
          size(Z21, 1),
          size(Z11, 1));
      end if;

    end dare;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

    function conditionNumber
      "Calculate the condition number norm(A)*norm(inv(A))"
      extends Modelica.Icons.Function;

      input Real A[:,:] "Input matrix";
      input Real p(min=1) = 2
        "Type of p-norm (only allowed: 1, 2 or Modelica.Constants.inf)";
      output Real result=0.0 "p-norm of matrix A";

    protected
      Real eps=1e-25;
      Real s[size(A, 1)] "singular values";

    algorithm
      if p == 2 then
        s := Modelica.Math.Matrices.singularValues(A);
        if min(s) < eps then
    result := Modelica.Constants.inf;
        else
    result := max(s)/min(s);
        end if;
      else
        result := Modelica.Math.Matrices.norm(A, p)*Modelica.Math.Matrices.norm(
          Modelica.Math.Matrices.inv(A), p);
      end if;

    end conditionNumber;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

    function det "Determinant of a matrix (computed by LU decomposition)"

      extends Modelica.Icons.Function;
      input Real A[:,size(A, 1)];
      output Real result "Determinant of matrix A";
    protected
      Real LU[size(A, 1),size(A, 1)];
      Integer pivots[size(A, 1)];

    algorithm
      if size(LU, 1) > 0 then
        (LU,pivots) := Modelica_LinearSystems2.Math.Matrices.LU(A);
        result := product(LU[i, i] for i in 1:size(A, 1))*product(if pivots[i] == i then
                1 else -1 for i in 1:size(pivots, 1));
      else
        result := -1e100;
      end if;
    end det;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

    function eigenValues
      "Compute eigenvalues and eigenvectors for a real, nonsymmetric matrix"

      extends Modelica.Icons.Function;

      input Real A[:, size(A, 1)] "Matrix";
      output Real eigenvalues[size(A, 1), 2]
        "Eigenvalues of matrix A (Re: first column, Im: second column)";

      output Real leftEigenvectors[size(A,1), size(A,2)]
        "Real-valued eigenvector matrix";
     output Real rightEigenvectors[size(A,1), size(A,2)]
        "Real-valued eigenvector matrix";

    protected
      Integer info;
      Boolean onlyEigenvalues = false;
    algorithm
    if size(A,1) > 0 then
      if onlyEigenvalues then
          (eigenvalues[:, 1],eigenvalues[:, 2],info) :=
            Modelica.Math.Matrices.LAPACK.dgeev_eigenValues(A);
         rightEigenvectors :=zeros(size(A, 1), size(A, 1));
         leftEigenvectors :=zeros(size(A, 1), size(A, 1));
      else
          (eigenvalues[:, 1],eigenvalues[:, 2],leftEigenvectors,rightEigenvectors,info) := Modelica_LinearSystems2.Math.Matrices.LAPACK.dgeev(A);
      end if;
      assert(info == 0, "Calculating the eigenvalues with function
    \"Matrices.eigenvalues\" is not possible, since the
    numerical algorithm does not converge.");
    end if;
    end eigenValues;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

    function equalityLeastSquares
      "Solve a linear equality constrained least squares problem"
      extends Modelica.Icons.Function;
      input Real A[:,:] "Minimize |A*x - a|^2";
      input Real a[size(A, 1)];
      input Real B[:,size(A, 2)] "subject to B*x=b";
      input Real b[size(B, 1)];
      output Real x[size(A, 2)] "solution vector";

    protected
      Integer info;
    algorithm
      if Modelica.Math.Matrices.rank(cat(
          1,
          A,
          B)) == size(A, 2) and min(size(A)) > 0 then

        (x,info) := Modelica.Math.Matrices.LAPACK.dgglse_vec(
          A,
          a,
          B,
          b);

        assert(info == 0, "Solving a linear equality-constrained least squares problem
    with function \"Matrices.equalityLeastSquares\" failed.");
      else
        if size(A, 2) == 0 then
          x := fill(0, 0);
        else
          x := Modelica_LinearSystems2.Math.Matrices.leastSquares(cat(
            1,
            A,
            B), cat(
            1,
            a,
            b));
        end if;
      end if;
    end equalityLeastSquares;

## Modelica_LinearSystems2.Math.Matrices.fliplr

flip the columns of a matrix in left/right direction

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, :] |  | Matrix to be fliped

### Outputs

Type | Name | Description
--- | --- | ---
Real | Aflip[size(A, 1), size(A, 2)] | fliped matrix

### Modelica definition

    function fliplr
      "flip the columns of a matrix in left/right direction"
      import Modelica_LinearSystems2.Math.Matrices;

      input Real A[:,:] "Matrix to be fliped";
      output Real Aflip[size(A, 1),size(A, 2)] "fliped matrix";
    algorithm
      for i in 1:size(A, 2) loop
        Aflip[:, i] := A[:, size(A, 2) + 1 - i];
      end for;
    end fliplr;

## Modelica_LinearSystems2.Math.Matrices.flipud

flip the columns of a matrix in up/down direction

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, :] |  | Matrix to be fliped

### Outputs

Type | Name | Description
--- | --- | ---
Real | Aflip[size(A, 1), size(A, 2)] | fliped matrix

### Modelica definition

    function flipud "flip the columns of a matrix in up/down direction"
      import Modelica_LinearSystems2.Math.Matrices;

      input Real A[:,:] "Matrix to be fliped";
      output Real Aflip[size(A, 1),size(A, 2)] "fliped matrix";
    algorithm
      for i in 1:size(A, 2) loop
        Aflip[i,:] := A[size(A, 2) + 1 - i,:];
      end for;
    end flipud;

## Modelica_LinearSystems2.Math.Matrices.fromFile

Read matrix from a matlab file

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | DataDir + "m.mat" |
String | matrixName | "m" | Name of the matrix

### Outputs

Type | Name | Description
--- | --- | ---
Real | A[n, m] |

### Modelica definition

    function fromFile "Read matrix from a matlab file"
      input String fileName=DataDir + "m.mat";
      input String matrixName="m" "Name of the matrix";

    protected
      input Integer Msizes[2]=readMatrixSize(fileName, matrixName);
      input Integer n=Msizes[1];
      input Integer m=Msizes[2];
      Real M[n,m]=readMatrix(
          fileName,
          matrixName,
          n,
          m);

    public
      output Real A[n,m]=M;
    algorithm

    end fromFile;

## Modelica_LinearSystems2.Math.Matrices.generalizedEigenvaluesTriangular

Compute invariant zeros of linear state space system with a genralized 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

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, size(A, 1)] |  |
Real | B[size(A, 1), size(A, 1)] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Real | alphaReal[size(A, 1)] | Real part of alpha (eigenvalue=(alphaReal+i*alphaImag)/beta)
Real | alphaImag[size(A, 1)] | Imaginary part of alpha
Real | beta[size(A, 1)] | Denominator of eigenvalue
Integer | info |

### Modelica definition

    function generalizedEigenvaluesTriangular
      "Compute invariant zeros of linear state space system with a genralized system matrix [A, B, C, D] which is of upper Hessenberg form"

      import Modelica_LinearSystems2.Math.Matrices.LAPACK;
      input Real A[:,size(A, 1)];
      input Real B[size(A, 1),size(A, 1)];
      output Real alphaReal[size(A, 1)]
        "Real part of alpha (eigenvalue=(alphaReal+i*alphaImag)/beta)";
      output Real alphaImag[size(A, 1)] "Imaginary part of alpha";
      output Real beta[size(A, 1)] "Denominator of eigenvalue";

      output Integer info;

    algorithm
      (alphaReal,alphaImag,beta,info) := LAPACK.dhgeqz(A, B);
      assert(info == 0, "Failed to compute eigenvalues with function dhgeqz(..)");
    end generalizedEigenvaluesTriangular;

## Modelica_LinearSystems2.Math.Matrices.hessenberg

Compute an upper Hessenberg matrix by repeatedly applicated householder similarity transformation

### Information

This function computes the Hessenberg matrix of matrix A by repetitive application of Householder similarity transformation

        Ai+1 = (I-2*u_i*u_i')*Ai*(I-2*u_i*u_i')

with Householder vector u_i.

The elementary transformations can be subsumed under

     A -> Q*A*Q

and Q*A*Q is Hessenberg matrix.

Modelica_LinearSystems2.Math.Matrices.hess uses LAPACK routine dgehrd. In contrast to this function Modelica_LinearSystems2.Math.Matrices.Internal.hessenberg does not use any LAPACK routine.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | H[size(H, 1), size(H, 2)] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Real | Ht[size(H, 1), size(H, 2)] |

### Modelica definition

    function hessenberg
      "Compute an upper Hessenberg matrix by repeatedly applicated householder similarity transformation"
      import Modelica;
      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Vectors;

      input Real H[size(H, 1),size(H, 2)];

      output Real Ht[size(H, 1),size(H, 2)];

    protected
      Integer q=size(H, 1);
      Real u[q] "householder vector";
      Integer ll;

    algorithm
      Ht := Modelica_LinearSystems2.Math.Matrices.toUpperHessenberg(
        H,
        1,
        size(H, 1));

    end hessenberg;

## Modelica_LinearSystems2.Math.Matrices.householderReflexion

reflect each of the vectors ai of matrix A=[a1, a2, ..., an] on a plane with orthogonal vector u

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, :] |  |
Real | u[size(A, 1)] |  | householder vector

### Outputs

Type | Name | Description
--- | --- | ---
Real | RA[size(A, 1), size(A, 2)] | reflexion of A

### Modelica definition

    function householderReflexion
      "reflect each of the vectors ai of matrix  A=[a1, a2, ..., an] on a plane with orthogonal vector u"
      import Modelica_LinearSystems2.Math.Vectors;

      input Real A[:,:];
      input Real u[size(A, 1)] "householder vector";

      output Real RA[size(A, 1),size(A, 2)] "reflexion of A";

    protected
      Integer n=size(A, 2);
      Real h;
      Real lu=Vectors.length(u)*Vectors.length(u);

    algorithm
      for i in 1:n loop
        h := scalar(2*transpose(matrix(u))*A[:, i]/lu);
        RA[:, i] := A[:, i] - h*u;
      end for;
    end householderReflexion;

## Modelica_LinearSystems2.Math.Matrices.householderSimilarityTransformation

Calculate the similarity transformation SAS of matrix A with householder matrix S = I - 2u*u'

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, size(A, 1)] |  |
Real | u[size(A, 1)] |  | householder vector

### Outputs

Type | Name | Description
--- | --- | ---
Real | SAS[size(A, 1), size(A, 1)] |

### Modelica definition

    function householderSimilarityTransformation
      "Calculate the similarity transformation SAS of matrix A with householder matrix S = I - 2u*u'"

      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Vectors;

      input Real A[:,size(A, 1)];
      input Real u[size(A, 1)] "householder vector";
      output Real SAS[size(A, 1),size(A, 1)];

    protected
      Integer na=size(A, 1);
      Real S[:,:]=-2*matrix(u)*transpose(matrix(u))/(Vectors.length(u)*
          Vectors.length(u));
      Integer i;
    algorithm
      for i in 1:na loop
                       //S=I-2u*u'
        S[i, i] := 1.0 + S[i, i];
      end for;
      SAS := S*A*S;

    end householderSimilarityTransformation;

## 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 a function).

### Inputs

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

### Outputs

Type | Name | Description
--- | --- | ---
Real | 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 definition

    function leastSquares
      "Solve overdetermined or underdetermined real system of linear equations A*x=b in a least squares sense (A may be rank deficient)"
      extends Modelica.Icons.Function;
      input Real A[:,:] "Matrix A";
      input Real b[size(A, 1)] "Vector b";
      output Real 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)";

    protected
      Integer info;
      Integer rank;
      Real xx[max(size(A, 1), size(A, 2))];
    algorithm
      if min(size(A)) > 0 then
        (xx,info,rank) := Modelica.Math.Matrices.LAPACK.dgelsx_vec(
          A,
          b,
          100*Modelica.Constants.eps);
        x := xx[1:size(A, 2)];
        assert(info == 0, "Solving an overdetermined or underdetermined linear system of
    equations with function \"Matrices.leastSquares\" failed.");
      else
        x := fill(0, size(A, 2));
      end if;
    end leastSquares;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

    function LU "LU decomposition of square or rectangular matrix"
      extends Modelica.Icons.Function;
      input Real A[:,:] "Square or rectangular matrix";
      output Real LU[size(A, 1),size(A, 2)]=A
        "L,U factors (used with LU_solve(..))";
      output Integer pivots[min(size(A, 1), size(A, 2))]
        "pivot indices (used with LU_solve(..))";
      output Integer info "Information";

    algorithm
      if min(size(A)) > 0 then
        (LU,pivots,info) := Modelica.Math.Matrices.LAPACK.dgetrf(A);
      else
        LU := fill(
          0,
          size(A, 1),
          size(A, 2));
        pivots := fill(0, min(size(A, 1), size(A, 2)));
        info := 0;
      end if;

    end LU;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

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

      extends Modelica.Icons.Function;
      input Real LU[:,size(LU, 1)]
        "L,U factors of Matrices.LU(..) for a square matrix";
      input Integer pivots[size(LU, 1)] "Pivots indices of Matrices.LU(..)";
      input Real b[size(LU, 1)] "Right hand side vector of P*L*U*x=b";
      output Real x[size(b, 1)] "Solution vector such that P*L*U*x = b";

    algorithm
      for i in 1:size(LU, 1) loop
        assert(LU[i, i] <> 0, "Solving a linear system of equations with function
    \"Matrices.LU_solve\" is not possible, since the LU decomposition
    is singular, i.e., no unique solution exists.");
      end for;
      if size(LU, 1) > 0 then
        x := Modelica.Math.Matrices.LAPACK.dgetrs_vec(
          LU,
          pivots,
          b);
      else
        x := fill(0, 0);
      end if;
    end LU_solve;

## 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

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

### Outputs

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

### Modelica definition

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

      input Real LU[:,size(LU, 1)]
        "L,U factors of Matrices.LU(..) for a square matrix";
      input Integer pivots[size(LU, 1)] "Pivots indices of Matrices.LU(..)";
      input Real B[size(LU, 1),:] "Right hand side matrix of P*L*U*X=B";
      output Real X[size(B, 1),size(B, 2)] "Solution matrix such that P*L*U*X = B";

    algorithm
      for i in 1:size(LU, 1) loop
        assert(LU[i, i] <> 0, "Solving a linear system of equations with function
    \"Matrices.LU_solve\" is not possible, since the LU decomposition
    is singular, i.e., no unique solution exists.");
      end for;
      if size(LU, 1) > 0 then
        X := LAPACK.dgetrs(
          LU,
          pivots,
          B);
      else
        X := fill(
          0,
          size(B, 1),
          size(B, 2));
      end if;
    end LU_solve2;

## 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

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, size(A, 1)] |  |
Real | C[size(A, 1), size(A, 2)] |  |
Real | eps | Modelica.Math.Matrices.norm(... |

### Outputs

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

### Modelica definition

    function lyapunov
      "Solution of continuous-time Lyapunov equation X*A + A'*X = C"
      import Modelica_LinearSystems2.Math.Matrices;

      input Real A[:,size(A, 1)];
      input Real C[size(A, 1),size(A, 2)];
      input Real eps=Modelica.Math.Matrices.norm(A,1)*10*Modelica.Constants.eps;

    protected
      Integer n=size(A, 1);
      Real R[size(A, 1),size(A, 2)] "rsf of A', i.e. R=U'A'U";
      Real U[size(A, 1),size(A, 2)] "transformation matrix U for R=U'A'U";
      Real C2[size(A, 1),size(A, 2)];
      Real R11[size(A, 1),size(A, 2)];
      Real R22[size(A, 1),size(A, 2)];
      Real R12[size(A, 1),size(A, 2)];
      Real R21[size(A, 1),size(A, 2)];
      Real R2[2*size(A, 1),2*size(A, 2)];
      Real I[size(A, 1),size(A, 1)]=identity(size(A, 1));
      Real x[2*size(A,1)];
      Real c[2*size(A,1)];
      Real CC[size(A,1),2];
      Integer k;

    public
      output Real X[size(A, 1),size(A, 2)] "solution of the Lyapunov equation";

    algorithm
      if n > 1 then
        (R,U) := Matrices.rsf(transpose(A));
        C2 := transpose(U)*C*U;
        X := zeros(n, n);

    // Calculate the last 1 or 2 columns of X
        R22 := R;
        for i1 in 1:n loop
          R22[i1, i1] := R[i1, i1] + R[n, n];
        end for;
        if abs(R[n, n - 1]) < eps then
          X[:, n] := Matrices.solve(R22, C2[:, n]);
          k := n - 1;
        else
          R11 := R;
          R12 := zeros(n, n);
          R21 := zeros(n, n);
          for i1 in 1:n loop
            R11[i1, i1] := R[i1, i1] + R[n - 1, n - 1];
            R12[i1, i1] := R[n - 1, n];
            R21[i1, i1] := R[n, n - 1];
          end for;

    // solve 2nx2n equation for 2x2 Schur bump with Kronecker product and vec operator approach
          R2 := [R11,R12; R21,R22];
          c := cat(1, C2[:, n - 1], C2[:, n]);
          x := Matrices.solve(R2, c);
          X[:, n - 1] := x[1:n];
          X[:, n] := x[n + 1:2*n];
          k := n - 2;
        end if;

    // Calculate the rest of X

        while k > 1 loop
          R22 := R;
          for i1 in 1:n loop
            R22[i1, i1] := R[i1, i1] + R[k, k];
          end for;
          if abs(R[k, k - 1]) < eps then
            //real eigenvalue
    //        X[:, k] := Matrices.solve(R22, C2[:, k] - vector(X[:, k + 1:n]*transpose(matrix(R[k, k + 1:n]))));
            X[:, k] := Matrices.solve(R22, C2[:, k] - vector(X[:, k + 1:n]*matrix(R[k, k + 1:n])));
            k := k - 1;
          else
           // conjugated complex eigenvalues
            R11 := R;
            R12 := zeros(n, n);
            R21 := zeros(n, n);
            for i1 in 1:n loop
              R11[i1, i1] := R[i1, i1] + R[k - 1, k - 1];
              R12[i1, i1] := R[k - 1, k];
              R21[i1, i1] := R[k, k - 1];
            end for;
            R2 := [R11,R12; R21,R22];
            CC := C2[:,k-1:k] - X[:,k+1:n]*transpose(R[k-1:k,k+1:n]);
            c := cat(1, CC[:, 1], CC[:, 2]);
            x := Matrices.solve(R2, c);
            X[:, k - 1] := x[1:n];
            X[:, k] := x[n + 1:2*n];

            k := k - 2;
          end if;
        end while;// k=1 or k=0

    // if k=1 the first column (if there exist a real eigenvalue) has to be calculated separately
        if k == 1 then
          R22 := R;
          for i1 in 1:n loop
            R22[i1, i1] := R[i1, i1] + R[1, 1];
          end for;
          X[:, 1] := Matrices.solve(R22, C2[:, 1] - vector(X[:, 2:n]*matrix(R[1, 2:n])));
        end if;

    // transform X corresponding to the original form
        X := U*X*transpose(U);

      elseif n == 1 then
        X[1, 1] := C[1, 1]/(2*A[1, 1]);
      else
        X := fill(0, 0, 0);
      end if;

    end lyapunov;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

    function norm "Returns the norm of a matrix"
      extends Modelica.Icons.Function;
      input Real A[:,:] "Input matrix";
      input Real p(min=1) = 2
        "Type of p-norm (only allowed: 1, 2 or Modelica.Constants.inf)";
      output Real result=0.0 "p-norm of matrix A";

    algorithm
      if p == 1 then
        // column sum norm
        for i in 1:size(A, 2) loop
          result := max(result, sum(abs(A[:, i])));
        end for;
      elseif p == 2 then
        // largest singular value
        result := if min(size(A)) > 0 then max(
          Modelica.Math.Matrices.singularValues(A)) else 0;
      elseif p == 3 then
        // Frobenius norm
        result :=   Modelica_LinearSystems2.Math.Matrices.Internal.frobeniusNorm(A);
      elseif p == Modelica.Constants.inf then
        // row sum norm
        for i in 1:size(A, 1) loop
          result := max(result, sum(abs(A[i, :])));
        end for;
      else
        assert(false, "Optional argument \"p\" of function \"norm\" must be
    1, 2 or Modelica.Constants.inf");
      end if;
    end norm;

## 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 dcomposed 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 a function).

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, :] |  | input matrix

### Outputs

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

### Modelica definition

    function nullspace "Orthonormal nullspace of a matrix"
      extends Modelica.Icons.Function;

      input Real A[:,:] "input matrix";
      output Real Z[size(A, 2),:] "orthonormal nullspace of matrix A";
      output Integer nullity "nullity, i.e. the dimension of the nullspace";

    protected
      Real V[size(A, 2),size(A, 2)] "Right orthogonal matrix ";
      Real sigma[min(size(A, 1), size(A, 2))] "singular values";
      Integer rank "rank of matrix A";
      Real eps "tolerance for rank determination";
      Integer n=min(size(A, 1), size(A, 2));
      Integer i=n;

    algorithm
      (sigma,,V) := Modelica.Math.Matrices.singularValues(A);
      V := transpose(V);
      // rank computation
      eps := max(size(A, 1), size(A, 2))*max(sigma)*Modelica.Constants.eps;
      rank := 0;
      if n > 0 then
        while i > 0 loop
          if sigma[i] > eps then
            rank := i;
            i := 0;
          end if;
          i := i - 1;
        end while;
      end if;
      Z := V[:,rank + 1:size(A,2)];// nullspace computation
      nullity := size(A,2) - rank;// nullity

    end nullspace;

## 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

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

### Outputs

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

### Modelica definition

    function 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"
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;
      input Real A[:,size(A, 1)];
      input Real tau[size(A, 1) - 1] "scalar factors of the elementary reflectors";
      input Integer ilo=1
        "lowest index where the original matrix had been Hessenbergform - ilo must have the same value as in the previous call of DGEHRD";
      input Integer ihi
        "highest index where the original matrix had been Hessenbergform  - ihi must have the same value as in the previous call of DGEHRD";

      output Real Q[size(A, 1),size(A, 2)]
        "Orthogonal matrix as a result of elementary reflectors";
      output Integer info;

    algorithm
      (Q,info) := LAPACK.dorghr(
        A,
        ilo,
        ihi,
        tau);
    end orthogonalQ;

## Modelica_LinearSystems2.Math.Matrices.printMatrix

print matrix

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | M[:, :] |  |
Integer | significantDigits | 6 | Number of significant digits that are shown
String | name | "M" | Independent variable name used for printing

### Outputs

Type | Name | Description
--- | --- | ---
String | s |

### Modelica definition

    function printMatrix "print matrix"
      import Modelica_LinearSystems2.StateSpace;
      import Modelica.Utilities.Strings;

      input Real M[:,:];
      input Integer significantDigits=6
        "Number of significant digits that are shown";
      input String name="M" "Independent variable name used for printing";
      output String s="";
    protected
      String blanks=Strings.repeat(significantDigits);
      String space=Strings.repeat(8);
      String space2=Strings.repeat(3);
      Integer r=size(M, 1);
      Integer c=size(M, 2);

    algorithm
      if r == 0 or c == 0 then
        s := name + " = []";
      else
        s := "\n" + name + " = \n";
        for i in 1:r loop
          s := s + space;
          for j in 1:c loop
            if M[i, j] >= 0 then
              s := s + " ";
            end if;
            s := s + String(M[i, j], significantDigits=significantDigits) +
              Strings.repeat(significantDigits + 8 - Strings.length(String(abs(M[i,j]))));

          end for;
          s := s + "\n";
        end for;

      end if;
    end printMatrix;

## Modelica_LinearSystems2.Math.Matrices.QR

QR decomposition of a square matrix without column pivoting (A = Q*R)

### Inputs

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

### Outputs

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

### Modelica definition

    function QR
      "QR decomposition of a square matrix without column pivoting (A = Q*R)"

      input Real A[:,:] "Rectangular matrix with size(A,1) >= size(A,2)";
      output Real Q[size(A, 1),size(A, 2)]
        "Rectangular matrix with orthonormal columns such that Q*R=A[:,p]";
      output Real R[min(size(A,1),size(A, 2)),size(A, 2)]
        "Square upper triangular matrix";

      output Real tau[min(size(A, 1), size(A, 2))];
      output Real Q2[:,:];

    protected
      Integer nrow=size(A, 1);
      Integer ncol=size(A, 2);
      Integer minrowcol=min(nrow,ncol);
      Integer lwork=Internal.dgeqrf_workdim(A);

    algorithm
      if minrowcol > 0 then

        (Q,tau) := Modelica_LinearSystems2.Math.Matrices.LAPACK.dgeqrf(A, lwork);

      // determine R
        R := zeros(minrowcol, ncol);
        for i in 1:minrowcol loop
          for j in i:ncol loop
            R[i, j] := Q[i, j];
          end for;
        end for;

    //  Q2 := Modelica.Math.Matrices.LAPACK.dorgqr(Q, tau);
    //  Q2 := Modelica_LinearSystems2.Math.Matrices.LAPACK.dorgqr(Q, tau);
        Q2 := Modelica_LinearSystems2.Math.Matrices.LAPACK.dorgqr_x(Q, tau);
      else
    //    Q := fill(1, size(A, 1), 0);
        Q := fill(1, size(A, 1), size(A, 2));
    //    R := fill(0, 0, 0);
        R := fill(0, min(size(A,1),size(A, 2)),size(A, 2));
      end if;
    end QR;

## Modelica_LinearSystems2.Math.Matrices.rcond

Reciprocal condition number

### Information

#### Syntax

>     r = Matrices.rcond(A);

#### 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

    function rcond "Reciprocal condition number"
      extends Modelica.Icons.Function;
      input Real A[:,size(A,1)] "Square real matrix";
      input Boolean inf = false
        "Is true if infinity norm is used and false for 1-norm";
      output Real rcond "Reciprocal condition number of A";
      output Integer info "Information";
    protected
      Real LU[:,:] "LU factorization of matrix A, returned by dgetrf";
      Real anorm "Norm of matrix A";
      String normspec= if inf then "I" else "1" "specifies the norm 1 or inf";

    algorithm
      if min(size(A)) > 0 then
        (LU,,info) := Modelica.Math.Matrices.LAPACK.dgetrf(A);
        anorm := Modelica_LinearSystems2.Math.Matrices.LAPACK.dlange(A,normspec);
        (rcond,info) := Modelica_LinearSystems2.Math.Matrices.LAPACK.dgecon(LU,inf,anorm);
      else
        rcond := Modelica.Constants.inf;
        info := 0;
      end if;

    end rcond;

## 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

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

### Outputs

Type | Name | Description
--- | --- | ---
Real | T[size(A, 1), size(A, 2)] |
Real | Z[size(A, 1), size(A, 2)] |
Real | alphaReal[size(A, 1)] | Real part of eigenvalue=alphaReal+i*alphaImag
Real | alphaImag[size(A, 1)] | Imaginary part of eigenvalue=(alphaReal+i*alphaImag

### Modelica definition

    function rsf "Computes the real Schur form (RSF) of a square matrix"
      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Matrices.Internal;
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;
      input Real A[:,size(A, 1)];

    protected
      Integer n=size(A, 1);
      Integer i;
      Integer lwork;
      Real tau[max(0, size(A, 1) - 1)];

      Real Aout[size(A, 1),size(A, 2)];
      Real H[size(A, 1),size(A, 2)] "upper Hessenberg form of A";
      Real Q[size(A, 1),size(A, 2)]
        "represents the Hessenberg transformation as a product of the elementary reflectors";
      Integer info1;
      Integer info2;

    public
      output Real T[size(A, 1),size(A, 2)];
      output Real Z[size(A, 1),size(A, 2)];
      output Real alphaReal[size(A, 1)]
        "Real part of eigenvalue=alphaReal+i*alphaImag";
      output Real alphaImag[size(A, 1)]
        "Imaginary part of eigenvalue=(alphaReal+i*alphaImag";

    algorithm
      if size(A, 1) > 1 then

    // dgehrd reduces a real matrix A to upper Hessenberg form Aout by an orthogonal
    // similarity transformation:  Q' * A * Q = Aout. Q can be computed from
    // Aout and tau (dorghr)

        (Aout,tau,info1) := LAPACK.dgehrd(
          A,
          1,
          n);
        assert(info1 == 0, "The " + String(-info1) +
          "'th argument of LAPACK.dgehrd had an illegal value");
    // dorghr to compute Q

        (Q,info1) := LAPACK.dorghr(
          Aout,
          1,
          n,
          tau);
        assert(info1 == 0, "The " + String(-info1) +
          "'th argument of LAPACK.dorghr had an illegal value");

        H[1:2, :] := Aout[1:2, :];
        for i in 3:n loop
          H[i, i - 1:n] := Aout[i, i - 1:n];
        end for;
    // dhseqr computes the eigenvalues of a real upper Hessenberg matrix H,
    // the Schur form T of H that is also the Schur form of A as well as the matrix
    // Z containing the Schur vectors to get A = Q*H*Q' = (Z)*T*(Z)'

        lwork := max(Internal.dhseqr_workdim(H), 1);
        (alphaReal,alphaImag,info2,T,Z) := LAPACK.dhseqr(
          H,
          lwork,
          false,
          "V",
          Q);
        assert(info2 == 0, "The output info of LAPACK.dhseqr should be zero, else if\n
         info < 0:  if info = -i, the i-th argument of LAPACK.dhseqr had an illegal value\n
         info > 0:  if INFO = i, LAPACK.dhseqr failed to compute all of the
                     eigenvalues in a total of 30*n iterations;\n
                     elements 1:n-1 and i+1:n of WR and WI contain those
                     eigenvalues which have been successfully computed.\n");
      else
        T := A;
        if size(A, 1) > 0 then
          Z := [1];
          alphaReal := {1};
          alphaImag := {0};
        else
          Z := fill(
            1,
            0,
            0);
          alphaReal := fill(1, 0);
          alphaImag := fill(0, 0);
        end if;
      end if;

    end rsf;

## 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. usng 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

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

### Outputs

Type | Name | Description
--- | --- | ---
Real | S[size(A, 1), size(A, 2)] |
Real | QZ[size(A, 1), size(A, 2)] |
Real | alphaReal[size(A, 1)] | Real part of eigenvalue=alphaReal+i*alphaImag
Real | alphaImag[size(A, 1)] | Imaginary part of eigenvalue=(alphaReal+i*alphaImag

### Modelica definition

    function rsf2
      "Computes the real Schur form (RSF) of a square matrix but uses lapack.dgees"
      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Matrices.Internal;
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input Real A[:,size(A, 1)];

    public
      output Real S[size(A, 1),size(A, 2)];
      output Real QZ[size(A, 1),size(A, 2)];
      output Real alphaReal[size(A, 1)]
        "Real part of eigenvalue=alphaReal+i*alphaImag";
      output Real alphaImag[size(A, 1)]
        "Imaginary part of eigenvalue=(alphaReal+i*alphaImag";

    protected
      Integer info;

    algorithm
      if size(A, 1) > 1 then
        (S, QZ, alphaReal, alphaImag, info) := Modelica_LinearSystems2.Math.Matrices.LAPACK.dgees(A);
         assert(info == 0, "The output info of LAPACK.dgees should be zero, else if\n
         info < 0:  if info = -i, the i-th argument of dgees had an illegal value\n
         info > 0:  if INFO = i, and i is
                   <= N: the QR algorithm failed to compute all the
                         eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
                         contain those eigenvalues which have converged; if
                         JOBVS = 'V', VS contains the matrix which reduces A
                         to its partially converged Schur form.\n");
      else
        S := A;
        if size(A, 1) > 0 then
          QZ := [1];
          alphaReal := {1};
          alphaImag := {0};
        else
          QZ := fill(
            1,
            0,
            0);
          alphaReal := fill(1, 0);
          alphaImag := fill(0, 0);
        end if;
      end if;

    end rsf2;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

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

      extends Modelica.Icons.Function;
      input Real A[:,size(A, 1)] "Matrix A of A*x = b";
      input Real b[size(A, 1)] "Vector b of A*x = b";
      output Real x[size(b, 1)] "Vector x such that A*x = b";

    protected
      Integer info;
    algorithm
      if size(A, 1) > 0 then
        (x,info) := Modelica.Math.Matrices.LAPACK.dgesv_vec(A, b);
        assert(info == 0, "Solving a linear system of equations with function
    \"Matrices.solve\" is not possible, because the system has either
    no or infinitely many solutions (A is singular).");
      else
        x := fill(0, 0);
      end if;
    end solve;

## 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 a function).

### Inputs

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

### Outputs

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

### Modelica definition

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

      extends Modelica.Icons.Function;
      input Real A[:,size(A, 1)] "Matrix A of A*X = B";
      input Real B[size(A, 1),:] "Matrix B of A*X = B";
      output Real X[size(B, 1),size(B, 2)] "Matrix X such that A*X = B";

    protected
      Integer info;
    algorithm
      if size(A, 1) > 0 then
        (X,info) := Modelica.Math.Matrices.LAPACK.dgesv(A, B);
        assert(info == 0, "Solving a linear system of equations with function
    \"Matrices.solve2\" is not possible, because the system has either
    no or infinitely many solutions (A is singular).");
      else
        X := fill(
          0,
          size(B, 1),
          size(B, 2));
      end if;
    end solve2;

## 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

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

### Outputs

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

### Modelica definition

    function sylvester
      "Solution of continuous-time Sylvester equation A*X + X*B = C"
      import Modelica_LinearSystems2.Math.Matrices;

      input Real A[:,:] "Matrix A";
      input Real B[:,:] "Matrix B";
      input Real C[size(A, 1),size(B, 2)] "Matrix C";
      input Boolean aIsSchur=false "True if A has already real Schur form";
      input Boolean bIsSchur=false "True if B has already real Schur form";
      output Real X[size(A, 1),size(B, 2)] "Solution of Sylvester equation";

    protected
      Integer n=size(A, 1);
      Integer m=size(B, 1);
      Real S[size(A, 1),size(A, 2)];
      Real T[size(B, 1),size(B, 2)];
      Real U[size(A, 1),size(A, 1)];
      Real V[size(B, 1),size(B, 1)];
      Real Chat[size(C, 1),size(C, 2)];
      Real scale;
      Integer info;

    algorithm
      if n > 1 and m > 1 then
        if aIsSchur then
          S := A;
          U := identity(n);
        else
          (S,U) := Matrices.rsf2(A);
        end if;
        if bIsSchur then
          T := B;
          V := identity(m);
        else
          (T,V) := Matrices.rsf2(B);
        end if;

        Chat := if aIsSchur and bIsSchur then C else if aIsSchur then C*V else if
          bIsSchur then transpose(U)*C else transpose(U)*C*V;
        (X,scale,info) := Matrices.LAPACK.dtrsyl(S, T, Chat);
        assert(info == 0, "Solving of Sylvester equation with Matrices.sylvester was not sucessfull.\n
                        The value of info is " + String(info) + ", but should be zero. A value unequal to zero means:\n
                < 0: if INFO = -i, the i-th argument had an illegal value\n
                = 1: A and B have common or very close eigenvalues; perturbed
                     values were used to solve the equation (but the matrices
                     A and B are unchanged).");
        X := if aIsSchur and bIsSchur then scale*X else if aIsSchur then scale*X*
          transpose(V) else if bIsSchur then scale*U*X else scale*U*X*transpose(V);
      else
        X := fill(0, n, m);
      end if;

    end sylvester;

## 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

This function uses DGEHRD-LAPACK routine to calculate a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation: Q' * A * Q = H.

Inputs are

- square matrix A
- indices ilo and ihi. See DGHRD description for further details of ilo and ihi

Outputs are

- Hessenberg matrix H = Q' * A * Q
- matrix V which contains the (Householder) vectors vi to perform the elementary reflectors Hi = I - tau * v * v' which finally define the matrix Q = product(Hi)
- vector tau with scalar weighting factors of the elementary reflectors

       =====================================================================
        -- LAPACK routine (version 3.0) --
          Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
          Courant Institute, Argonne National Lab, and Rice University
          June 30, 1999

       Purpose
       =======

       DGEHRD reduces a real general matrix A to upper Hessenberg form H by
       an orthogonal similarity transformation:  Q' * A * Q = H .

       Arguments
       =========

       N       (input) INTEGER
               The order of the matrix A.  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. ILO and IHI are normally
               set by a previous call to DGEBAL; otherwise they should be
               set to 1 and N respectively. See Further Details.
               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 general matrix to be reduced.
               On exit, the upper triangle and the first subdiagonal of A
               are overwritten with the upper Hessenberg matrix H, and the
               elements below the first subdiagonal, with the array TAU,
               represent the orthogonal matrix Q as a product of elementary
               reflectors. See Further Details.

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

       TAU     (output) DOUBLE PRECISION array, dimension (N-1)
               The scalar factors of the elementary reflectors (see Further
               Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
               zero.

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

       LWORK   (input) INTEGER
               The length of the array WORK.  LWORK >= max(1,N).
               For optimum performance LWORK >= N*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.

       Further Details
       ===============

       The matrix Q is represented as a product of (ihi-ilo) elementary
       reflectors

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

       Each H(i) has the form

          H(i) = I - tau * v * v'

       where tau is a real scalar, and v is a real vector with
       v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
       exit in A(i+2:ihi,i), and tau in TAU(i).

       The contents of A are illustrated by the following example, with
       n = 7, ilo = 2 and ihi = 6:

       on entry,                        on exit,

       ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
       (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
       (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
       (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
       (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
       (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
       (                         a )    (                          a )

       where a denotes an element of the original matrix A, h denotes a
       modified element of the upper Hessenberg matrix H, and vi denotes an
       element of the vector defining H(i).

       =====================================================================

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, size(A, 1)] |  |
Integer | ilo | 1 | lowest index where the original matrix had been Hessenbergform
Integer | ihi |  |

### Outputs

Type | Name | Description
--- | --- | ---
Real | H[size(A, 1), size(A, 2)] | Upper Hessenberg form
Real | V[size(A, 1), size(A, 2)] | V=[v1,v2,..vn-1,0] with vi are vectors which define the elementary reflectors
Real | tau[size(A, 1) - 1] |
Integer | info |

### Modelica definition

    function toUpperHessenberg
      "transform a real general matrix A to upper Hessenberg form H by an orthogonal similarity transformation:  Q' * A * Q = H"
      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;
      input Real A[:,size(A, 1)];
      input Integer ilo=1
        "lowest index where the original matrix had been Hessenbergform";
      input Integer ihi;//=size(A, 1) "highest index where the original matrix had been Hessenbergform";
      output Real H[size(A, 1),size(A, 2)] "Upper Hessenberg form";
      output Real V[size(A, 1),size(A, 2)]
        "V=[v1,v2,..vn-1,0] with vi are vectors which define the elementary reflectors";

      output Real tau[size(A, 1) - 1];
      output Integer info;

    protected
      Integer n=size(A, 1);
      Real Aout[size(A, 1),size(A, 2)];
      Integer i;

    algorithm
      (Aout,tau,info) := LAPACK.dgehrd(
        A,
        ilo,
        ihi);
      H[1:2, 1:ihi] := Aout[1:2, 1:ihi];
      H[1:2, ihi + 1:n] := A[1:2, ihi + 1:n];

      for i in 3:n loop
        H[i, i - 1:ihi] := Aout[i, i - 1:ihi];
        H[i, ihi + 1:n] := A[i, ihi + 1:n];
      end for;

      for i in 1:min(n - 2, ihi) loop
        V[i + 1, i] := 1.0;
        V[i + 2:n, i] := Aout[i + 2:n, i];

      end for;
      V[n, n - 1] := 1;

    end toUpperHessenberg;

## Modelica_LinearSystems2.Math.Matrices.trace

tarce(A) is the sum of the diagonal elements of A

### Information

Extends from Modelica.Icons.Function (Icon for a function).

### Inputs

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

### Outputs

Type | Name | Description
--- | --- | ---
Real | result |

### Modelica definition

    function trace "tarce(A) is the sum of the diagonal elements of A"
      extends Modelica.Icons.Function;

      input Real A[:,size(A, 1)];
      output Real result;

    protected
      Integer n=size(A, 1);
      Real r;

    algorithm
      r := 0;
      if n > 0 then
        for i in 1:n loop
          r := r + A[i, i];
        end for;
      end if;
      result := r;
    end trace;