Riccati Matrix Equation Solvers

Standard Riccati Matrix Equations

MatrixEquations.arecFunction
arec(A, G, Q = 0; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - XGX + Q = 0,

where G and Q are hermitian/symmetric matrices or uniform scaling operators. Scalar-valued G and Q are interpreted as appropriately sized uniform scaling operators G*I and Q*I. The Schur method of [1] is used.

To enhance the accuracy of computations, a block scaling of matrices G and Q is performed, if the default setting scaling = 'B' is used. This scaling is however performed only if norm(Q) > norm(G). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, experimental structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-GX.

Z = [U; V] is an orthogonal basis for the stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi, where Sx and Sxi are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx and scalinfo.Sxi, respectively.

Note: To solve the continuous-time algebraic Riccati equation

 A'X + XA - XBR^(-1)B'X + Q = 0,

with R a hermitian/symmetric matrix and B a compatible size matrix, G = BR^(-1)B' must be provided. This approach is not numerically suited when R is ill-conditioned and/or B has large norm.

Reference:

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

Example

julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> G = [1. 0. 0.; 0. 5. 0.; 0. 0. 10.]
3×3 Array{Float64,2}:
 1.0  0.0   0.0
 0.0  5.0   0.0
 0.0  0.0  10.0

julia> X, CLSEIG = arec(A,G,2I);

julia> X
3×3 Array{Float64,2}:
  0.459589   0.333603   -0.144406
  0.333603   0.65916    -0.0999216
 -0.144406  -0.0999216   0.340483

julia> A'*X+X*A-X*G*X+2I
3×3 Array{Float64,2}:
  2.22045e-16  4.44089e-16  -1.77636e-15
  4.44089e-16  6.66134e-16   1.11022e-16
 -1.77636e-15  1.11022e-16  -1.33227e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
 -4.411547592296008 + 2.4222082620381102im
 -4.411547592296008 - 2.4222082620381102im
 -4.337128244724371 + 0.0im

julia> eigvals(A-G*X)
3-element Array{Complex{Float64},1}:
 -4.4115475922960075 - 2.4222082620381076im
 -4.4115475922960075 + 2.4222082620381076im
  -4.337128244724374 + 0.0im
source
arec(A, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, orth = false, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - (XB+S)R^(-1)(B'X+S') + Q = 0,

where R and Q are hermitian/symmetric matrices or uniform scaling operators such that R is nonsingular. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)). The Schur method of [1] is used.

To enhance the accuracy of computations, a block scaling of matrices R, Q and S is performed, if the default setting scaling = 'B' is used. This scaling is however performed only if norm(Q) > norm(B)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D' or scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-BF.

F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S').

Z = [U; V; W] is a basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively. An orthogonal basis Z can be determined, with an increased computational cost, by setting orth = true.

Reference:

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

Example

julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
 1.0  2.0
 2.0  0.0
 0.0  1.0

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = arec(A,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.522588   0.303007  -0.327227
  0.303007   0.650895  -0.132608
 -0.327227  -0.132608   0.629825

julia> A'*X+X*A-X*B*inv(R)*B'*X+2I
3×3 Array{Float64,2}:
 -2.66454e-15  -1.55431e-15   8.88178e-16
 -1.55431e-15   2.22045e-15  -2.9976e-15
  9.99201e-16  -2.9976e-15    4.44089e-16

julia> CLSEIG
3-element Array{Complex{Float64},1}:
   -4.37703628399912 + 2.8107164873731247im
   -4.37703628399912 - 2.8107164873731247im
 -1.8663764577096091 + 0.0im

julia> eigvals(A-B*F)
3-element Array{Complex{Float64},1}:
  -4.377036283999118 - 2.8107164873731234im
  -4.377036283999118 + 2.8107164873731234im
 -1.8663764577096063 + 0.0im
source
arec(A, B, G, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, orth = false, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Computes X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the continuous-time algebraic Riccati equation

 A'X + XA - XGX - (XB+S)R^(-1)(B'X+S') + Q = 0,

where G, R and Q are hermitian/symmetric matrices or uniform scaling operators such that R is nonsingular. Scalar-valued G, R and Q are interpreted as appropriately sized uniform scaling operators G*I, R*I and Q*I. S, if not specified, is set to S = zeros(size(B)). For well conditioned R, the Schur method of [1] is used. For ill-conditioned R or if orth = true, the generalized Schur method of [2] is used.

To enhance the accuracy of computations, a block oriented scaling of matrices G, R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R)). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. If orth = true, two experimental scaling procedures can be activated using the options scaling = 'D' and scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'.

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) eigenvalues of A-BF-GX.

F is the stabilizing or anti-stabilizing gain matrix F = R^(-1)(B'X+S').

Z = [U; V; W] is a basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively. An orthogonal basis Z can be determined, with an increased computational cost, by setting orth = true.

Reference:

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

[2] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

source
MatrixEquations.aredFunction
ared(A, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the discrete-time algebraic Riccati equation

A'XA - X - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where R and Q are hermitian/symmetric matrices. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D', scaling = 'R' and scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable) eigenvalues of A-BF.

F is the stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = Sx*(V/U)*Sxi and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra

julia> A = [ 0. 1.; 0. 0. ]
2×2 Array{Float64,2}:
 0.0  1.0
 0.0  0.0

julia> B = [ 0.; sqrt(2.) ]
2-element Array{Float64,1}:
 0.0
 1.4142135623730951

julia> R = 1.
1.0

julia> Q = [ 1. -1.; -1. 1. ]
2×2 Array{Float64,2}:
  1.0  -1.0
 -1.0   1.0

julia> X, CLSEIG, F = ared(A,B,R,Q);

julia> X
2×2 Array{Float64,2}:
  1.0  -1.0
 -1.0   1.5

julia> A'*X*A-X-A'*X*B*inv(R+B'*X*B)*B'*X*A+Q
2×2 Array{Float64,2}:
  0.0          -3.33067e-16
 -3.33067e-16   8.88178e-16

julia> CLSEIG
2-element Array{Complex{Float64},1}:
 0.4999999999999998 - 0.0im
               -0.0 - 0.0im

julia> eigvals(A-B*F)
2-element Array{Float64,1}:
 -2.7755575615628914e-16
  0.5
source

Generalized Riccati Matrix Equations

MatrixEquations.garecFunction
garec(A, E, G, Q = 0; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - E'XGXE + Q = 0,

where G and Q are hermitian/symmetric matrices or uniform scaling operators and E is a nonsingular matrix. Scalar-valued G and Q are interpreted as appropriately sized uniform scaling operators G*I and Q*I. The generalized Schur method of [1] is used.

To enhance the accuracy of computations, a block scaling of matrices G and Q is performed, if the default setting scaling = 'B' is used. This scaling is however performed only if norm(Q) > norm(G). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-GXE,E).

Z = [U; V] is an orthogonal basis for the stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E, where Sx and Sxi are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx and scalinfo.Sxi, respectively.

Note: To solve the continuous-time algebraic Riccati equation

 A'XE + E'XA - E'XBR^(-1)B'XE + Q = 0,

with R a hermitian/symmetric matrix and B a compatible size matrix, G = BR^(-1)B' must be provided. This approach is not numerically suited when R is ill-conditioned and/or B has large norm.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

source
garec(A, E, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where R and Q are hermitian/symmetric matrices such that R is nonsingular, and E is a nonsingular matrix. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)). The generalized Schur method of [1] is used.

To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D' or scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E).

F is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> E = [10. 3. 0.; 0. 5. -1.; 0. 0. 10.]
3×3 Array{Float64,2}:
 10.0  3.0   0.0
  0.0  5.0  -1.0
  0.0  0.0  10.0

julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
 1.0  2.0
 2.0  0.0
 0.0  1.0

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = garec(A,E,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.0502214   0.0284089   -0.0303703
  0.0284089   0.111219    -0.00259162
 -0.0303703  -0.00259162   0.0618395

julia> A'*X*E+E'*X*A-E'*X*B*inv(R)*B'*X*E+2I
3×3 Array{Float64,2}:
  1.55431e-15  -1.9984e-15   -3.33067e-15
 -1.77636e-15   1.33227e-15  -3.33067e-15
 -2.88658e-15  -3.21965e-15   1.11022e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
  -0.6184265391601464 + 0.2913286844595737im
  -0.6184265391601464 - 0.2913286844595737im
 -0.21613059964451786 + 0.0im

julia> eigvals(A-B*F,E)
3-element Array{Complex{Float64},1}:
 -0.6184265391601462 - 0.29132868445957383im
 -0.6184265391601462 + 0.2913286844595739im
  -0.216130599644518 + 0.0im
source
garec(A, E, B, G, R, Q, S; scaling = 'B', pw2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized continuous-time algebraic Riccati equation

A'XE + E'XA - E'XGXE - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where G, Q and R are hermitian/symmetric matrices such that R is nonsingular, and E is a nonsingular matrix. Scalar-valued G, R and Q are interpreted as appropriately sized uniform scaling operators G*I, R*I and Q*I. The generalized Schur method of [1] is used.

To enhance the accuracy of computations, a block oriented scaling of matrices G, R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > max(norm(G), norm(B)^2/norm(R)). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D' or scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF-GXE,E).

F is the stabilizing/anti-stabilizing gain matrix F = R^(-1)(B'XE+S').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

source
MatrixEquations.garedFunction
gared(A, E, B, R, Q, S; scaling = 'B', pow2 = false, as = false, rtol::Real = nϵ, nrm = 1) -> (X, EVALS, F, Z, scalinfo)

Compute X, the hermitian/symmetric stabilizing solution (if as = false) or anti-stabilizing solution (if as = true) of the generalized discrete-time algebraic Riccati equation

A'XA - E'XE - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where R and Q are hermitian/symmetric matrices, and E ist non-singular. Scalar-valued R and Q are interpreted as appropriately sized uniform scaling operators R*I and Q*I. S, if not specified, is set to S = zeros(size(B)).

To enhance the accuracy of computations, a block oriented scaling of matrices R, Q and S is performed using the default setting scaling = 'B'. This scaling is performed only if norm(Q) > norm(B)^2/norm(R). A general, eigenvalue computation oriented scaling combined with a block scaling is used if scaling = 'G' is selected. An alternative, structure preserving scaling can be performed using the option scaling = 'S'. A symmetric matrix equilibration based scaling is employed if scaling = 'K', for which the underlying vector norm can be specified using the keyword argument nrm = p, where p = 1 is the default setting. Experimental structure preserving scalings can be performed using the options scaling = 'D', scaling = 'R' and scaling = 'T'. Scaling can be disabled with the choice scaling = 'N'. If pow2 = true, the scaling elements are enforced to the nearest power of 2 (default: pow2 = false).

By default, the lower bound for the 1-norm reciprocal condition number rtol is n*ϵ, where n is the order of A and ϵ is the machine epsilon of the element type of A.

EVALS is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair (A-BF,E).

F is the stabilizing/anti-stabilizing gain matrix F = (R+B'XB)^(-1)(B'XA+S').

Z = [U; V; W] is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that X = (Sx*(V/U)*Sxi)/E and F = -Sr*(W/U)*Sxi, where Sx, Sxi and Sr are diagonal scaling matrices contained in the named tuple scalinfo as scalinfo.Sx, scalinfo.Sxi and scalinfo.Sr, respectively.

Reference:

[1] W.F. Arnold, III and A.J. Laub, Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations, Proc. IEEE, 72:1746-1754, 1984.

Example

julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> E = [10. 3. 0.; 0. 5. -1.; 0. 0. 10.]
3×3 Array{Float64,2}:
 10.0  3.0   0.0
  0.0  5.0  -1.0
  0.0  0.0  10.0

julia> B = [1. 2.; 2. 0.; 0. 1.]
3×2 Array{Float64,2}:
 1.0  2.0
 2.0  0.0
 0.0  1.0

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = gared(A,E,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.065865   -0.0147205  -0.0100407
 -0.0147205   0.0885939   0.0101422
 -0.0100407   0.0101422   0.0234425

julia> A'*X*A-E'*X*E-A'*X*B*inv(R+B'*X*B)*B'*X*A+2I
3×3 Array{Float64,2}:
 -1.33227e-15  -2.48412e-15   1.38778e-16
 -2.498e-15    -4.44089e-16  -6.50521e-16
  1.80411e-16  -5.91541e-16  -1.33227e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
  -0.084235615751339 - 0.0im
  -0.190533552034239 - 0.0im
 -0.5238922629921539 - 0.0im

julia> eigvals(A-B*F,E)
3-element Array{Float64,1}:
 -0.5238922629921539
 -0.19053355203423886
 -0.08423561575133902
source