## Modelica_LinearSystems2.Math.Matrices.Internal

### Package Content

Name | Description
--- | ---
carenls | Newton's method with exact line search for solving continuous algebraic riccati equation
dgeqp3_workdim | Calculate the optimal size of the WORK array in dgeqp3
dgeqrf_workdim | Calculate the optimal size of the WORK array in dgeqrf
dhseqr_workdim | Calculate the optimal size of the WORK array in dhseqr
eigenvalues2 | Compute eigenvalues and unnormalized eigenvectors
eigenvaluesHessenberg | Compute eigenvalues of an upper Hessenberg form matrix
findLocal_tk | Find a local minimizer tk to define the length of the step tk*Nk in carenls
haveZeroRow | Boolean output is true if at least one matrix row is zero vector
hessenberg2 | Compute an upper Hessenberg matrix by repeatedly applicated householder similarity transformation
hohoTrafoLowerHess | Compute the similarity transformation S*A*S of matrix A with householder matrix S = I - 2u*u' to compute a lower Hessenberg form
hohoTrafoUpperHess | Compute the similarity (Householder-) transformation S*A*S of matrix A with householder matrix S = I - 2u*u' to compute an upper Hessenberg form
householderSimilarityTransformation2 | Calculate the similarity transformation SAS of matrix A with householder matrix S = I - 2u*u'/(u'*u) to compute a lower Hessenberg form
multiplyWithOrthogonalQ_hr | overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', depending on inputs trans and side
multiplyWithOrthogonalQ_qr | overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', where Q is an orthogonal matrix of a QR factorization as returned by dgeqrf
QR | QR decomposition of a square matrix with column pivoting (A(:,p) = Q*R)
QR2 | QR decomposition of a square matrix with column pivoting (A(:,p) = Q*R). Uses dgeqpf instead of dgeqp3
readMatrixGain | Read a matrix from mat-file
reorderRSF | Reorders a real Schur factorization according to a given pattern of the eigenvalues
reorderRSF2 | Reorders a real Schur factorization for poleAssigmentMI design
k_care_u | calculate the upper bound of the CARE, i.e. Q + A'*X + X*A - X*G*X = 0 condition number using Lyapunov equations
frobeniusNorm | Return the Frobenius norm of a matrix

## Modelica_LinearSystems2.Math.Matrices.Internal.carenls

Newton's method with exact line search for solving continuous algebraic riccati equation

### Information

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

### 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, 2)] | identity(size(A, 1)) |
Real | X0[size(A, 1), size(A, 2)] |  |
Real | eps | Modelica_LinearSystems2.Math... |

### Outputs

Type | Name | Description
--- | --- | ---
Real | X[size(X0, 1), size(X0, 2)] |
Real | r |

### Modelica definition

    function carenls
      "Newton's method with exact line search for solving continuous algebraic riccati equation"
      extends Modelica.Icons.Function;

      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Polynomial;
      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, 2)]=identity(size(A, 1));
      input Real X0[size(A, 1),size(A, 2)];
      input Real eps=Modelica_LinearSystems2.Math.Matrices.Internal.frobeniusNorm(
                                            A)*1e-9;

      output Real X[size(X0, 1),size(X0, 2)];
      output Real r;

    protected
      Integer n=size(A, 1);
      Real G[size(A, 1),size(A, 2)]=B*Modelica.Math.Matrices.solve2(R, transpose(B));
      Real Xk[size(X, 1),size(X, 2)];
      Real Ak[size(A, 1),size(A, 2)];
      Real Rk[size(A, 1),size(A, 2)];
      Real Nk[size(A, 1),size(A, 2)];
      Real tk;
      Integer k;
      Complex xc[2];
      Boolean stop;

    algorithm
      if n > 1 then
        k := 0;
        stop := false;
        Xk := X0;
        while (not stop and k<10) loop
          k := k + 1;
          Ak := A - G*Xk;
          Rk := transpose(A)*Xk + Xk*A + Q - Xk*G*Xk;
          Nk := Matrices.lyapunov(Ak, -Rk);
          tk := Matrices.Internal.findLocal_tk(Rk, G, Nk);
          stop := eps > Modelica_LinearSystems2.Math.Matrices.Internal.frobeniusNorm(
                                               tk*Nk)/Modelica_LinearSystems2.Math.Matrices.Internal.frobeniusNorm(
                                                                             Xk);
          Xk := Xk + tk*Nk;
        end while;
        X := Xk;
        r := Modelica_LinearSystems2.Math.Matrices.Internal.frobeniusNorm(
                                    X*A + transpose(A)*X - X*G*X + Q);

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

    end carenls;

## Modelica_LinearSystems2.Math.Matrices.Internal.dgeqp3_workdim

Calculate the optimal size of the WORK array in dgeqp3

### Inputs

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

### Outputs

Type | Name | Description
--- | --- | ---
Integer | lwork |
Integer | info |

### Modelica definition

    function dgeqp3_workdim
      "Calculate the optimal size of the WORK array in dgeqp3"
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input Real A[:,:];

      output Integer lwork;
      output Integer info;

    protected
      Real work[:];

    algorithm
      (,,,info,work) := LAPACK.dgeqp3(A, -1);
      lwork := integer(work[1]);

    end dgeqp3_workdim;

## Modelica_LinearSystems2.Math.Matrices.Internal.dgeqrf_workdim

Calculate the optimal size of the WORK array in dgeqrf

### Inputs

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

### Outputs

Type | Name | Description
--- | --- | ---
Integer | lwork |
Integer | info |

### Modelica definition

    function dgeqrf_workdim
      "Calculate the optimal size of the WORK array in dgeqrf"
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input Real A[:,:];

      output Integer lwork;
      output Integer info;

    protected
      Real work[:];

    algorithm
      (,,info,work) := LAPACK.dgeqrf(A, -1);
      lwork := integer(work[1]);

    end dgeqrf_workdim;

## Modelica_LinearSystems2.Math.Matrices.Internal.dhseqr_workdim

Calculate the optimal size of the WORK array in dhseqr

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | H[:, :] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Integer | lwork |
Integer | info |

### Modelica definition

    function dhseqr_workdim
      "Calculate the optimal size of the WORK array in dhseqr"
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input Real H[:,:];

      output Integer lwork;
      output Integer info;

    protected
      Real work[:];

    algorithm
      if min(size(H, 1), size(H, 2)) > 0 then
        (,,info,,,work) := LAPACK.dhseqr(H, -1);
        lwork := integer(work[1]);
      else
        lwork := 1;
      end if;

    end dhseqr_workdim;

## Modelica_LinearSystems2.Math.Matrices.Internal.eigenvalues2

Compute eigenvalues and unnormalized eigenvectors

### Information

### Inputs

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

### Outputs

Type | Name | Description
--- | --- | ---
Real | alphaReal[size(A, 1)] | Real part of alpha (eigenvalue=(alphaReal+i*alphaImag))
Real | alphaImag[size(A, 1)] | Imaginary part of alpha (eigenvalue=(alphaReal+i*alphaImag))
Real | lEigenVectors[size(A, 1), size(A, 1)] | left eigenvectors of matrix A
Real | rEigenVectors[size(A, 1), size(A, 1)] | right eigenvectors of matrix A
Integer | info |

### Modelica definition

    function eigenvalues2
      "Compute eigenvalues and unnormalized eigenvectors"

      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Matrices.Internal;
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;

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

      output Real alphaReal[size(A, 1)]
        "Real part of alpha (eigenvalue=(alphaReal+i*alphaImag))";
      output Real alphaImag[size(A, 1)]
        "Imaginary part of alpha (eigenvalue=(alphaReal+i*alphaImag))";
      output Real lEigenVectors[size(A, 1),size(A, 1)]
        "left eigenvectors of matrix A";
      output Real rEigenVectors[size(A, 1),size(A, 1)]
        "right eigenvectors of matrix A";
      output Integer info=0;
    protected
      Integer n=size(A, 1);
      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];
      Real Q[size(A, 1),size(A, 2)];
      Real Zo[:,:];
      Real Ho[:,:];

      Integer lwork;

    algorithm
      (H,V,tau) := Modelica_LinearSystems2.Math.Matrices.toUpperHessenberg(
                                            A,1,n);
      Q := LAPACK.dorghr(V, 1, n, tau);

     lwork := Internal.dhseqr_workdim(H);

      if size(H, 1) > 0 then
       (alphaReal,alphaImag,info,Ho,Zo) := LAPACK.dhseqr(H, lwork, false, "V", Q);
      else
        alphaReal := fill(0, size(H, 1));
        alphaImag := fill(0, size(H, 1));
        Zo := fill(0, size(H, 1), size(H, 2));
        Ho := fill(0, size(H, 1), size(H, 2));
      end if;

      (lEigenVectors,rEigenVectors) := LAPACK.dtrevc(Ho, "B", "B", Zo);

    end eigenvalues2;

## Modelica_LinearSystems2.Math.Matrices.Internal.eigenvaluesHessenberg

Compute eigenvalues of an upper Hessenberg form matrix

### Information

This function uses DHSEQR Lapack-routine to calculate the eigenvalues of an upper Hessenberg form H. Therefore, H is reduced to Schur form T. The eigenvalues are obtained from the diagonal of T.

See Modelica_LinearSystems2.Math.Matrices.LAPACK.dhseqr for details

### Inputs

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

### Outputs

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

### Modelica definition

    function eigenvaluesHessenberg
      "Compute eigenvalues of an upper Hessenberg form matrix"
      import Modelica_LinearSystems2.Math.Matrices.Internal;
      import Modelica_LinearSystems2.Math.Matrices.LAPACK;
      input Real H[:,size(H, 1)];

      output Real alphaReal[size(H, 1)]
        "Real part of alpha (eigenvalue=(alphaReal+i*alphaImag))";
      output Real alphaImag[size(H, 1)]
        "Imaginary part of alpha (eigenvalue=(alphaReal+i*alphaImag))";
      output Integer info=0;
    protected
      Integer n=size(H, 1);
      Integer ilo=1;
      Integer ihi=n;
      Integer lwork=Internal.dhseqr_workdim(H);
      Real work[lwork];
      Real Z[n,n]=fill(
            0,
            n,
            n);

    algorithm
      if size(H, 1) > 0 then
        (alphaReal,alphaImag,info) := LAPACK.dhseqr(H, lwork);
      else
        alphaReal := fill(0, size(H, 1));
        alphaImag := fill(0, size(H, 1));
      end if;

    end eigenvaluesHessenberg;

## Modelica_LinearSystems2.Math.Matrices.Internal.findLocal_tk

Find a local minimizer tk to define the length of the step tk*Nk in carenls

### Information

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

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | Rk[:, size(Rk, 2)] |  |
Real | G[size(Rk, 1), size(Rk, 2)] |  |
Real | Nk[size(Rk, 1), size(Rk, 2)] |  |

### Outputs

Type | Name | Description
--- | --- | ---
Real | tk |

### Modelica definition

    function findLocal_tk
      "Find a local minimizer tk to define the length of the step tk*Nk in carenls"
      extends Modelica.Icons.Function;

      import Modelica_LinearSystems2.Math.Matrices;
      import Modelica_LinearSystems2.Math.Polynomial;
      import Modelica_LinearSystems2.Math.Complex;

      input Real Rk[:,size(Rk, 2)];
      input Real G[size(Rk, 1),size(Rk, 2)];
      input Real Nk[size(Rk, 1),size(Rk, 2)];

      output Real tk;

    protected
      Real Vk[size(Rk, 1),size(Rk, 2)];
      Real alpha_k;
      Real beta_k;
      Real gamma_k;
      Complex p[3];
      Boolean h;

    algorithm
      Vk := Nk*G*Nk;
      alpha_k := Matrices.trace(Rk*Rk);
      beta_k := Matrices.trace(Rk*Vk);
      gamma_k := Matrices.trace(Vk*Vk);

      if gamma_k > Modelica.Constants.eps then
        p := Polynomial.roots(Polynomial({4*gamma_k,6*beta_k,2*(alpha_k - 2*beta_k),
          -2*alpha_k}));
        h := false;
        for i1 in 1:3 loop
          if (abs(Complex.imag(p[i1])) < Modelica.Constants.eps) then
            if (abs(Complex.real(p[i1]) - 1) <= 1) then
              tk := Complex.real(p[i1]);
              h := true;
            end if;
          end if;
        end for;
        if not h then
          tk := 1;
        end if;

      else
        tk := 1;
      end if;

    end findLocal_tk;

## Modelica_LinearSystems2.Math.Matrices.Internal.haveZeroRow

Boolean output is true if at least one matrix row is zero vector

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, :] | fill(0, 0, 0) | input matrix

### Outputs

Type | Name | Description
--- | --- | ---
Boolean | result | is true if A has at least one zero row

### Modelica definition

    function haveZeroRow
      "Boolean output is true if at least one matrix row is zero vector"

      input Real A[:,:]=fill(
            0,
            0,
            0) "input matrix";
      output Boolean result "is true if A has at least one zero row";

    protected
      Boolean h;
      Integer i;
    algorithm

      i := 1;
      h := false;
      while i < size(A, 1) + 1 and not h loop
        h := Modelica.Math.Vectors.isEqual(
            zeros(size(A, 2)),
            A[i, :],
            100*Modelica.Constants.eps);
        i := i + 1;
      end while;
      result := h;

    end haveZeroRow;

## Modelica_LinearSystems2.Math.Matrices.Internal.hessenberg2

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.

In contrast to function Modelica_LinearSystems2.Math.Matrices.hess, function Modelica_LinearSystems2.Math.Matrices.hess3 does not use any LAPACK routine.

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | H[:, :] |  |
String | s | "u" |

### Outputs

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

### Modelica definition

    function hessenberg2
      "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[:,:];
      input String s="u";
      output Real Ht[size(H, 1),size(H, 2)];

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

    algorithm
      assert(s == "u" or s == "l",
        "parameter s should be equal to 'u' or 'l' to indicate upper or lower Hessenberg form");

      Ht := H;

      for ll in 1:q - 2 loop
        u := if s == "u" then cat(
          1,
          zeros(ll),
          cat(1, Vectors.householderVector(vector(Ht[ll + 1:q, ll]), cat(
            1,
            {1},
            zeros(q - ll - 1))))) else cat(
          1,
          cat(1, Vectors.householderVector(vector(Ht[1:q - ll, q - ll + 1]), cat(
            1,
            zeros(q - ll - 1),
            {1}))),
          zeros(ll));

        Ht := if s == "u" then
          Modelica_LinearSystems2.Math.Matrices.Internal.hohoTrafoUpperHess(
          Ht,
          u,
          ll) else
          Modelica_LinearSystems2.Math.Matrices.Internal.hohoTrafoLowerHess(
          Ht,
          u,
          ll);

      end for;

    end hessenberg2;

## Modelica_LinearSystems2.Math.Matrices.Internal.hohoTrafoLowerHess

Compute the similarity transformation S*A*S of matrix A with householder matrix S = I - 2u*u' to compute a lower Hessenberg form

### Information

This function calculates one step in the calculation of lower Hessenberg form. Therein it calculates the (n-r+1)'th column of the Hessenberg matrix which is of shape {fill(1,n-r-1),x,x}.

### Inputs

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

### Outputs

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

### Modelica definition

    function hohoTrafoLowerHess
      "Compute the similarity transformation S*A*S of matrix A with householder matrix S = I - 2u*u' to compute a lower Hessenberg form"

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

      input Real A[:,size(A, 1)];
      input Real u[size(A, 1)] "householder vector";
      input Integer r;
      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));   //S=u*u'/u'*u
      Integer i;

      Real P[na - r,na - r];
      Real A11[na - r,na - r]=A[1:na - r, 1:na - r];
      Real A21[r,na - r]=A[na - r + 1:na, 1:na - r];
      Real A22[r,r]=A[na - r + 1:na, na - r + 1:na];
      Real alpha;

    algorithm
      for i in 1:na loop
        S[i, i] := 1.0 + S[i, i];   //S=I-2u*u'
      end for;

      P := S[1:na - r, 1:na - r];
      alpha := P[na - r, :]*A[1:na - r, na - r + 1];

      SAS := [P*A11*P,[zeros(na - r - 1, r); matrix(alpha),zeros(1, r - 1)]; A21*P,
        A22];

    end hohoTrafoLowerHess;

## Modelica_LinearSystems2.Math.Matrices.Internal.hohoTrafoUpperHess

Compute the similarity (Householder-) transformation S*A*S of matrix A with householder matrix S = I - 2u*u' to compute an upper Hessenberg form

### Information

This function calculates one step in the calculation of upper Hessenberg form. Therein it calculates the r'th column of the Hessenberg matrix which is of shape {x,x,fill(1,n-r-1)}.

    </html>

### Inputs

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

### Outputs

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

### Modelica definition

    function hohoTrafoUpperHess
      "Compute the similarity (Householder-) transformation S*A*S of matrix A with householder matrix S = I - 2u*u' to compute an upper Hessenberg form"

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

      input Real A[:,size(A, 1)];
      input Real u[size(A, 1)] "householder vector";
      input Integer r;
      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));                                                             //S=u*u'/u'*u
      Integer i;

      Real P[na - r,na - r];
      Real A11[r,r]=A[1:r, 1:r];
      Real A12[r,na - r]=A[1:r, r + 1:na];
      Real A22[na - r,na - r]=A[r + 1:na, r + 1:na];

      Real alpha;

    algorithm
      for i in 1:na loop
        S[i, i] := 1.0 + S[i, i];
                                //S=I-2u*u'
      end for;

      P := S[r + 1:na, r + 1:na];
      alpha := P[1, :]*A[r + 1:na, r];

      SAS := [A11,A12*P; [zeros(1, r - 1),matrix(alpha); zeros(na - r - 1, r)],P*
        A22*P];

    end hohoTrafoUpperHess;

## Modelica_LinearSystems2.Math.Matrices.Internal.householderSimilarityTransformation2

Calculate the similarity transformation SAS of matrix A with householder matrix S = I - 2u*u'/(u'*u) to compute a lower Hessenberg form

### Information

It is assumed, that the input vector u is a Housholder vector of the shape

>     u = (u1, u2, ..., ur,0, ..., 0)

where r is an integer input. From

>     S = I - 2*u*u'/u'*u = [P, 0; 0, I]

with

>     dim(P) = r x r,   dim(I) = n-r x n-r

results

>     S*A*S = [P*A11*P, P*A12; A21*P, A22]

### Inputs

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

### Outputs

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

### Modelica definition

    function householderSimilarityTransformation2
      "Calculate the similarity transformation SAS of matrix A with householder matrix S = I - 2u*u'/(u'*u) to compute a lower Hessenberg form"

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

      input Real A[:,size(A, 1)];
      input Real u[size(A, 1)] "householder vector";
      input Integer r;
      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));                                                             //S=u*u'/u'*u
      Integer i;

    algorithm
      assert(r <= na,
        "Input r in function \"Matrices.Internal.householderSimilarityTransformation2\" must fulfill r<=size(A,1)");
      for i in 1:na loop
        S[i, i] := 1.0 + S[i, i];   //S=I-2u*u'
      end for;

      SAS := if r < size(A, 1) then [S[1:r, 1:r]*A[1:r, 1:r]*S[1:r, 1:r],[
        zeros(r - 1, 1),A[1:r - 1, r + 2:na]; matrix(S[r, 1:r]*A[1:r, r + 1]),
        transpose(matrix(A[r, r + 2:na]))]; A[r + 1:na, 1:r]*S[1:r, 1:r],A[
        r + 1:na, r + 1:na]] else S*A*S;

    end householderSimilarityTransformation2;

## Modelica_LinearSystems2.Math.Matrices.Internal.multiplyWithOrthogonalQ_hr

overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', depending on inputs trans and side

### Information

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

          .. Scalar Arguments ..
          ..
          .. Array Arguments ..
          ..

       Purpose
       =======

         DORMHR overwrites the general real M-by-N matrix C with

                       SIDE = 'L'     SIDE = 'R'
       TRANS = 'N':      Q * C          C * Q
       TRANS = 'T':      Q**T * C       C * Q**T

       where Q is a real orthogonal matrix of order nq, with nq = m if
       SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
       IHI-ILO elementary reflectors, as returned by DGEHRD:

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

       Arguments
       =========

       SIDE    (input) CHARACTER*1
               = 'L': apply Q or Q**T from the Left;
               = 'R': apply Q or Q**T from the Right.

       TRANS   (input) CHARACTER*1
               = 'N':  No transpose, apply Q;
               = 'T':  Transpose, apply Q**T.

       M       (input) INTEGER
               The number of rows of the matrix C. M >= 0.

       N       (input) INTEGER
               The number of columns of the matrix C. 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).
               If SIDE = 'L', then 1 <= ILO <= IHI <= M, if M > 0, and
               ILO = 1 and IHI = 0, if M = 0;
               if SIDE = 'R', then 1 <= ILO <= IHI <= N, if N > 0, and
               ILO = 1 and IHI = 0, if N = 0.

       A       (input) DOUBLE PRECISION array, dimension
                                    (LDA,M) if SIDE = 'L'
                                    (LDA,N) if SIDE = 'R'
               The vectors which define the elementary reflectors, as
               returned by DGEHRD.

       LDA     (input) INTEGER
               The leading dimension of the array A.
               LDA >= max(1,M) if SIDE = 'L'; LDA >= max(1,N) if SIDE = 'R'.

       TAU     (input) DOUBLE PRECISION array, dimension
                                    (M-1) if SIDE = 'L'
                                    (N-1) if SIDE = 'R'
               TAU(i) must contain the scalar factor of the elementary
               reflector H(i), as returned by DGEHRD.

       C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
               On entry, the M-by-N matrix C.
               On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

       LDC     (input) INTEGER
               The leading dimension of the array C. LDC >= max(1,M).

       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.
               If SIDE = 'L', LWORK >= max(1,N);
               if SIDE = 'R', LWORK >= max(1,M).
               For optimum performance LWORK >= N*NB if SIDE = 'L', and
               LWORK >= M*NB if SIDE = 'R', 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 | C[:, :] |  |
Real | A[:, :] |  |
Real | tau[size(A, 2) - 1] |  |
String | side | "L" |
String | trans | "N" |
Integer | ilo | 1 | lowest index where the original matrix had been Hessenbergform
Integer | ihi | size(A, 2) | highest index where the original matrix had been Hessenbergform

### Outputs

Type | Name | Description
--- | --- | ---
Real | Cout[size(C, 1), size(C, 2)] | contains the Hessenberg form in the upper triangle and the first subdiagonal and below the first subdiagonal it contains the elementary reflectors which represents (with array tau) as a product the orthogonal matrix Q
Integer | info |

### Modelica definition

    function multiplyWithOrthogonalQ_hr
      "overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', depending on inputs trans and side"

      input Real C[:,:];
      input Real A[:,:];
      input Real tau[size(A, 2) - 1];
      input String side="L";
      input String trans="N";
      input Integer ilo=1
        "lowest index where the original matrix had been Hessenbergform";
      input Integer ihi=size(A, 2)
        "highest index where the original matrix had been Hessenbergform";
      output Real Cout[size(C, 1),size(C, 2)]=C
        "contains the Hessenberg form in the upper triangle and the first subdiagonal and below the first subdiagonal it contains the elementary reflectors which represents (with array tau) as a product the orthogonal matrix Q";

      output Integer info;
    protected
      Integer m=size(C, 1);
      Integer n=size(C, 2);
      Integer lda=max(1, size(A, 2));
      Integer ldc=max(1, m);
      Integer lwork=2*size(A, 2);
      Real work[lwork];

    external "Fortran 77" dormhr(
          side,
          trans,
          m,
          n,
          ilo,
          ihi,
          A,
          lda,
          tau,
          Cout,
          ldc,
          work,
          lwork,
          info);

    end multiplyWithOrthogonalQ_hr;

## Modelica_LinearSystems2.Math.Matrices.Internal.multiplyWithOrthogonalQ_qr

overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', where Q is an orthogonal matrix of a QR factorization as returned by dgeqrf

### Information

       Purpose
       =======

       DORMQR overwrites the general real M-by-N matrix C with

                       SIDE = 'L'     SIDE = 'R'
       TRANS = 'N':      Q * C          C * Q
       TRANS = 'T':      Q**T * C       C * Q**T

       where Q is a real orthogonal matrix defined as the product of k
       elementary reflectors

             Q = H(1) H(2) . . . H(k)

       as returned by DGEQRF. Q is of order M if SIDE = 'L' and of order N
       if SIDE = 'R'.

       Arguments
       =========

       SIDE    (input) CHARACTER*1
               = 'L': apply Q or Q**T from the Left;
               = 'R': apply Q or Q**T from the Right.

       TRANS   (input) CHARACTER*1
               = 'N':  No transpose, apply Q;
               = 'T':  Transpose, apply Q**T.

       M       (input) INTEGER
               The number of rows of the matrix C. M >= 0.

       N       (input) INTEGER
               The number of columns of the matrix C. N >= 0.

       K       (input) INTEGER
               The number of elementary reflectors whose product defines
               the matrix Q.
               If SIDE = 'L', M >= K >= 0;
               if SIDE = 'R', N >= K >= 0.

       A       (input) DOUBLE PRECISION array, dimension (LDA,K)
               The i-th column must contain the vector which defines the
               elementary reflector H(i), for i = 1,2,...,k, as returned by
               DGEQRF in the first k columns of its array argument A.
               A is modified by the routine but restored on exit.

       LDA     (input) INTEGER
               The leading dimension of the array A.
               If SIDE = 'L', LDA >= max(1,M);
               if SIDE = 'R', LDA >= max(1,N).

       TAU     (input) DOUBLE PRECISION array, dimension (K)
               TAU(i) must contain the scalar factor of the elementary
               reflector H(i), as returned by DGEQRF.

       C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
               On entry, the M-by-N matrix C.
               On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.

       LDC     (input) INTEGER
               The leading dimension of the array C. LDC >= max(1,M).

       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.
               If SIDE = 'L', LWORK >= max(1,N);
               if SIDE = 'R', LWORK >= max(1,M).
               For optimum performance LWORK >= N*NB if SIDE = 'L', and
               LWORK >= M*NB if SIDE = 'R', 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 | C[:, :] |  |
Real | Q[:, :] |  |
Real | tau[:] |  |
String | side | "L" |
String | trans | "N" |

### Outputs

Type | Name | Description
--- | --- | ---
Real | Cout[size(C, 1), size(C, 2)] | contains Q*C or Q**T*C or C*Q**T or C*Q
Integer | info |

### Modelica definition

    function multiplyWithOrthogonalQ_qr
      "overwrites the general real M-by-N matrix C with Q * C or C * Q or Q' * C or C * Q', where Q is an orthogonal matrix of a QR factorization as returned by dgeqrf"

      input Real C[:,:];
      input Real Q[:,:];
      input Real tau[:];
      input String side="L";
      input String trans="N";

      output Real Cout[size(C, 1),size(C, 2)]=C
        "contains Q*C or Q**T*C or C*Q**T or C*Q";
      output Integer info;
    protected
      Integer m=size(C, 1);
      Integer n=size(C, 2);
      Integer k=if side == "L" then m else n;
      Integer lda=if side == "L" then max(1, m) else max(1, n);
      Integer ldc=max(1, m);
      Integer lwork=if side == "L" then max(1, n) else max(1, m);
      Real work[lwork];

    external "Fortran 77" dormqr(
          side,
          trans,
          m,
          n,
          k,
          Q,
          lda,
          tau,
          Cout,
          ldc,
          work,
          lwork,
          info);

    end multiplyWithOrthogonalQ_qr;

## Modelica_LinearSystems2.Math.Matrices.Internal.QR

QR decomposition of a square matrix with column pivoting (A(:,p) = Q*R)

### Information

#### Syntax

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

#### Description

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

> Q*R = A[:,p]

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

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

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

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

where R1 is a regular, upper triangular matrix.

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

#### Example

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

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

### 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[size(A, 2), size(A, 2)] | Square upper triangular matrix
Integer | p[size(A, 2)] | Column permutation vector

### Modelica definition

    function QR
      "QR decomposition of a square matrix with column pivoting (A(:,p) = Q*R)"

      extends Modelica.Icons.Function;
      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[size(A, 2),size(A, 2)] "Square upper triangular matrix";
      output Integer p[size(A, 2)] "Column permutation vector";

    protected
      Integer nrow=size(A, 1);
      Integer ncol=size(A, 2);
      Real tau[ncol];
    algorithm
      assert(nrow >= ncol, "\nInput matrix A[" + String(nrow) + "," + String(ncol)
         + "] has more columns as rows.
    This is not allowed when calling Modelica.Matrices.QR(A).");
      if ncol > 0 then
        (Q,tau,p) := Modelica.Math.Matrices.LAPACK.dgeqpf(A);

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

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

## Modelica_LinearSystems2.Math.Matrices.Internal.QR2

QR decomposition of a square matrix with column pivoting (A(:,p) = Q*R). Uses dgeqpf instead of dgeqp3

### Information

#### Syntax

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

#### Description

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

> Q*R = A[:,p]

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

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

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

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

where R1 is a regular, upper triangular matrix.

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

#### Example

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

### Inputs

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[size(A, 2), size(A, 2)] | Square upper triangular matrix
Integer | p[size(A, 2)] | Column permutation vector
Real | tau[min(size(A, 1), size(A, 2))] |

### Modelica definition

    function QR2
      "QR decomposition of a square matrix with column pivoting (A(:,p) = Q*R). Uses dgeqpf instead of dgeqp3"

      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[size(A, 2),size(A, 2)] "Square upper triangular matrix";
      output Integer p[size(A, 2)] "Column permutation vector";
      output Real tau[min(size(A, 1), size(A, 2))];

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

    algorithm
      assert(nrow >= ncol, "\nInput matrix A[" + String(nrow) + "," + String(ncol)
         + "] has more columns as rows.
    This is not allowed when calling Modelica.Matrices.QR(A).");
      if ncol > 0 then

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

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

        Q := Modelica.Math.Matrices.LAPACK.dorgqr(Q, tau);
      else
        Q := fill(
          1,
          size(A, 1),
          0);
        R := fill(
          0,
          0,
          0);
      end if;
    end QR2;

## Modelica_LinearSystems2.Math.Matrices.Internal.readMatrixGain

Read a matrix from mat-file

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
String | fileName | "matrixGain.mat" | Name of the matrix data file
String | matrixName | "K" | Name of the matrix
Integer | m |  |
Integer | n |  |

### Outputs

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

### Modelica definition

    encapsulated function readMatrixGain "Read a matrix from mat-file"

      import Modelica;
      import Modelica_LinearSystems2;

     input String fileName="matrixGain.mat" "Name of the matrix data file";
     input String matrixName="K" "Name of the matrix";
     input Integer m;
     input Integer n;

    public
       output Real K[m,n];

    algorithm
     K := Modelica_LinearSystems2.Internal.Streams.readMatrixInternal(fileName, matrixName, m, n);

    end readMatrixGain;

## Modelica_LinearSystems2.Math.Matrices.Internal.reorderRSF

Reorders a real Schur factorization according to a given pattern of the eigenvalues

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Boolean | iscontinuous |  |
Real | T[:, :] |  |
Real | Q[:, size(T, 2)] |  |
Real | alphaReal[size(T, 1)] |  | Real part of eigenvalue=alphaReal+i*alphaImag
Real | alphaImag[size(T, 1)] |  | Imaginary part of eigenvalue=(alphaReal+i*alphaImag

### Outputs

Type | Name | Description
--- | --- | ---
Real | To[size(T, 1), size(T, 2)] |
Real | Qo[size(T, 1), size(T, 2)] |
Real | wr[size(T, 2)] |
Real | wi[size(T, 2)] |

### Modelica definition

    function reorderRSF
      "Reorders a real Schur factorization according to a given pattern of the eigenvalues"

      import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input Boolean iscontinuous;
      input Real T[:,:];
      input Real Q[:,size(T, 2)];
      input Real alphaReal[size(T, 1)]
        "Real part of eigenvalue=alphaReal+i*alphaImag";
      input Real alphaImag[size(T, 1)]
        "Imaginary part of eigenvalue=(alphaReal+i*alphaImag";

      output Real To[size(T, 1),size(T, 2)];
      output Real Qo[size(T, 1),size(T, 2)];
      output Real wr[size(T, 2)];
      output Real wi[size(T, 2)];

    protected
      Integer n=size(T, 2);
      Boolean select[:]=fill(false, size(T, 2));
      Integer i;
    algorithm
      if iscontinuous then
        for i in 1:n loop
          if alphaReal[i] < 0 then
            select[i] := true;
          end if;
        end for;
      else
        for i in 1:n loop
          if alphaReal[i]^2 + alphaImag[i]^2 < 1 then
            select[i] := true;
          end if;
        end for;
      end if;

      (To,Qo,wr,wi) := LAPACK.dtrsen(
          "E",
          "V",
          select,
          T,
          Q);

    end reorderRSF;

## Modelica_LinearSystems2.Math.Matrices.Internal.reorderRSF2

Reorders a real Schur factorization for poleAssigmentMI design

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | T[:, :] |  | upper quasi-triangular matrix in Schur canonical form
Real | Q[:, size(T, 2)] |  | matrix of Schur vectors
Real | alphaReal[size(T, 1)] |  | Real part of eigenvalue=alphaReal+i*alphaImag
Real | alphaImag[size(T, 1)] |  | Imaginary part of eigenvalue=(alphaReal+i*alphaImag
Real | alpha |  | maximum admissible value for real parts(continuous) or for moduli (discrete) of the eigenvalues of A which will not be modified by the eigenvalue assignment algorithm

### Outputs

Type | Name | Description
--- | --- | ---
Real | To[size(T, 1), size(T, 2)] |
Real | Qo[size(T, 1), size(T, 2)] |
Real | wr[size(T, 2)] |
Real | wi[size(T, 2)] |

### Modelica definition

    function reorderRSF2
      "Reorders a real Schur factorization for poleAssigmentMI design"

      import Modelica_LinearSystems2.Math.Matrices.LAPACK;

      input Real T[:,:] "upper quasi-triangular matrix in Schur canonical form";
      input Real Q[:,size(T, 2)] "matrix of Schur vectors";
      input Real alphaReal[size(T, 1)]
        "Real part of eigenvalue=alphaReal+i*alphaImag";
      input Real alphaImag[size(T, 1)]
        "Imaginary part of eigenvalue=(alphaReal+i*alphaImag";
      input Real alpha
        "maximum admissible value for real parts(continuous) or for moduli (discrete) of the eigenvalues of A which will not be modified by the eigenvalue assignment algorithm";

      output Real To[size(T, 1),size(T, 2)];
      output Real Qo[size(T, 1),size(T, 2)];
      output Real wr[size(T, 2)];
      output Real wi[size(T, 2)];

    protected
      Integer n=size(T, 2);
      Boolean select[n]=fill(false, n);
      Integer i;
    algorithm
      for i in 1:n loop
        if alphaReal[i] < alpha then
          select[i] := true;
        end if;
      end for;

      (To,Qo,wr,wi) := LAPACK.dtrsen(
          "E",
          "V",
          select,
          T,
          Q);

    end reorderRSF2;

## Modelica_LinearSystems2.Math.Matrices.Internal.k_care_u

calculate the upper bound of the CARE, i.e. Q + A'*X + X*A - X*G*X = 0 condition number using Lyapunov equations

### Information

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

### Inputs

Type | Name | Default | Description
--- | --- | --- | ---
Real | A[:, size(A, 1)] |  | care-matrix A
Real | Q[:, size(Q, 1)] |  | care-matrix Q
Real | G[size(A, 1), size(A, 2)] |  | care-matrix G
Real | X[size(A, 1), size(A, 2)] |  | solution of care

### Outputs

Type | Name | Description
--- | --- | ---
Real | ku | upper bound of the care condition number

### Modelica definition

    function k_care_u
      "calculate the upper bound of the CARE, i.e. Q + A'*X + X*A - X*G*X = 0  condition number using Lyapunov equations"
      extends Modelica.Icons.Function;

      import Modelica_LinearSystems2.Math.Matrices;

      input Real A[:,size(A, 1)] "care-matrix A";
      input Real Q[:,size(Q, 1)] "care-matrix Q";
      input Real G[size(A, 1),size(A, 2)] "care-matrix G";
      input Real X[size(A, 1),size(A, 2)] "solution of care";

      output Real ku "upper bound of the care condition number";

    protected
      Real Z0[size(A, 1),size(A, 2)]=Matrices.lyapunov(A - G*X, -identity(size(A, 1)))
        "solution of lyapunov equation H'*Z0+Z0*H=-I, H=A-G*X";
      Real Z2[size(A, 1),size(A, 2)]=Matrices.lyapunov(A - G*X, -X*X)
        "solution of lyapunov equation H'*Z2+Z2*H=-X*X, H=A-G*X";
      Real normA=Matrices.norm(A, 2) "spectral norm of matrix A (largest sv)";
      Real normQ=Matrices.norm(Q, 2) "spectral norm of matrix Q (largest sv)";
      Real normG=Matrices.norm(G, 2) "spectral norm of matrix G (largest sv)";
      Real normX=Matrices.norm(X, 2) "spectral norm of matrix X (largest sv)";
      Real normZ0=Matrices.norm(Z0, 2) "spectral norm of matrix Z0 (largest sv)";
      Real normZ2=Matrices.norm(Z2, 2) "spectral norm of matrix Z2 (largest sv)";
    algorithm
      ku := (normZ0*normQ + 2*sqrt(normZ0*normZ2)*normA + normZ2*normG)/normX;
    end k_care_u;

## Modelica_LinearSystems2.Math.Matrices.Internal.frobeniusNorm

Return the Frobenius norm of a matrix

### Information

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

### Inputs

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

### Outputs

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

### Modelica definition

    function frobeniusNorm "Return the Frobenius norm of a matrix"
      extends Modelica.Icons.Function;
      input Real A[:,:] "Input matrix";
      output Real result=0.0 "frobenius norm of matrix A";

    algorithm
      for i1 in 1:size(A, 1) loop
        for i2 in 1:size(A, 2) loop
          result := result + A[i1, i2]*A[i1, i2];
        end for;
      end for;
      result := sqrt(result);
    end frobeniusNorm;