Upgrade to V1.2.1

LICENSE.md

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2019 andreasvarga
+Copyright (c) 2019-2020 Andreas Varga
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

Project.toml

@@ -1,7 +1,7 @@
 name = "MatrixEquations"
 uuid = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
 authors = ["Andreas Varga <varga.andreas@gmail.com>"]
-version = "1.2.0"
+version = "1.2.1"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

docs/src/index.md

@@ -68,7 +68,12 @@ The collection of tools can be extended by adding new functionality, such as exp
 
 ## Release Notes
 
-### Versions 1.2.0
+### Version 1.2.1
+
+Patch release to address fallback issues to ensure compatibility to versions prior 1.3 of Julia,
+some enhancements of the 2x2 positive generalized Lyapunov equation solver, explicit handling of null dimension case in Riccati solvers.
+
+### Version 1.2.0
 
 Minor release targeting sensible (up to 50%) speed increase of various lower level solvers for Lyapunov and Sylvester equations. This goal has been achieved by the reduction of allocation burden using preallocation of small size work arrays, explicit forming of small order Kronecker product based coefficient matrices, performing updating operations with the 5-term `mul!` function introduced in `Julia 1.3` (compatibility with prior Julia versions ensured using calls to BLAS `gemm!`).  The functionality of lower level solvers has been strictly restricted to the basic real and complex data of types `BlasReal` and `BlasComplex`.
 

src/MatrixEquations.jl

@@ -34,4 +34,26 @@ include("sylvkr.jl")
 include("plyapunov.jl")
 include("meoperators.jl")
 include("condest.jl")
+# fallback for versions prior 1.3
+if VERSION < v"1.3.0" 
+    mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasReal} = 
+                            BLAS.gemm!('N', 'N', α, A, B, β, C)
+    mul!(C::StridedMatrix{T}, adjA::Transpose{T,<:StridedMatrix{T}}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasReal} = 
+                            BLAS.gemm!('T', 'N', α, parent(adjA), B, β, C)
+    mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, adjB::Transpose{T,<:StridedMatrix{T}}, α::T, β::T) where {T<:BlasReal} = 
+                            BLAS.gemm!('N', 'T', α, A, parent(adjB),  β, C)
+    mul!(C::StridedMatrix{T}, adjA::Transpose{T,<:StridedMatrix{T}}, adjB::Transpose{T,<:StridedMatrix{T}}, α::T, β::T) where {T<:BlasReal} = 
+                            BLAS.gemm!('T', 'T', α, parent(adjA), parent(adjB),  β, C)
+    # mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasComplex} = 
+    #                         BLAS.gemm!('N', 'N', α, A, B, β, C)
+    # mul!(C::StridedMatrix{T}, adjA::Adjoint{T,<:StridedMatrix{T}}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasComplex} = 
+    #                         BLAS.gemm!('C', 'N', α, parent(adjA), B, β, C)
+    # mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, adjB::Adjoint{T,<:StridedMatrix{T}}, α::T, β::T) where {T<:BlasComplex} = 
+    #                         BLAS.gemm!('N', 'C', α, A, parent(adjB),  β, C)
+    # mul!(C::StridedMatrix{T}, adjA::Adjoint{T,<:StridedMatrix{T}}, adjB::Adjoint{T,<:StridedMatrix{T}}, α::T, β::T) where {T<:BlasComplex} = 
+    #                         BLAS.gemm!('C', 'C', α, parent(adjA), parent(adjB),  β, C)
+    mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}) where {T<:BlasReal} = 
+         mul!(C,A,B,one(T),zero(T))
+ end
+
 end

src/lyapunov.jl

@@ -1745,24 +1745,3 @@ function lyapds!(A::Matrix{T1},E::Union{Matrix{T1},UniformScaling{Bool}}, C::Mat
       end
    end
 end
-# fallback for versions prior 1.3
-if VERSION < v"1.3.0" 
-   mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasReal} = 
-                           BLAS.gemm!('N', 'N', α, A, B, β, C)
-   mul!(C::StridedMatrix{T}, adjA::Transpose{T,<:StridedMatrix{T}}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasReal} = 
-                           BLAS.gemm!('T', 'N', α, parent(adjA), B, β, C)
-   mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, adjB::Transpose{T,<:StridedMatrix{T}}, α::T, β::T) where {T<:BlasReal} = 
-                           BLAS.gemm!('N', 'T', α, A, parent(adjB),  β, C)
-   mul!(C::StridedMatrix{T}, adjA::Transpose{T,<:StridedMatrix{T}}, adjB::Transpose{T,<:StridedMatrix{T}}, α::T, β::T) where {T<:BlasReal} = 
-                           BLAS.gemm!('T', 'T', α, parent(adjA), parent(adjB),  β, C)
-   # mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasComplex} = 
-   #                         BLAS.gemm!('N', 'N', α, A, B, β, C)
-   # mul!(C::StridedMatrix{T}, adjA::Adjoint{T,<:StridedMatrix{T}}, B::StridedMatrix{T}, α::T, β::T) where {T<:BlasComplex} = 
-   #                         BLAS.gemm!('C', 'N', α, parent(adjA), B, β, C)
-   # mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, adjB::Adjoint{T,<:StridedMatrix{T}}, α::T, β::T) where {T<:BlasComplex} = 
-   #                         BLAS.gemm!('N', 'C', α, A, parent(adjB),  β, C)
-   # mul!(C::StridedMatrix{T}, adjA::Adjoint{T,<:StridedMatrix{T}}, adjB::Adjoint{T,<:StridedMatrix{T}}, α::T, β::T) where {T<:BlasComplex} = 
-   #                         BLAS.gemm!('C', 'C', α, parent(adjA), parent(adjB),  β, C)
-   mul!(C::AbstractMatrix{T}, A::AbstractMatrix{T}, B::AbstractMatrix{T}) where {T<:BlasReal} = 
-        mul!(C,A,B,one(T),zero(T))
-end

src/plyapunov.jl

@@ -2161,6 +2161,7 @@ function pglyap2!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, R::AbstractMatri
    # This function is based on the SLICOT routine SG03BX, which implements the
    # generalization of the method due to Hammarling ([1], section 6) for Lyapunov
    # equations of order 2. A more detailed description is given in [2].
+   # The 2x2 matrix is allowed to be a full matrix.
 
    # [1] Hammarling S. J.
    #     Numerical solution of the stable, non-negative definite Lyapunov equation.
@@ -2215,11 +2216,20 @@ function pglyap2!(A::AbstractMatrix{T1}, E::AbstractMatrix{T1}, R::AbstractMatri
       RR[1,1] = RR[2,2]
       RR[2,2] = V
    end
-   scale1, scale2, LAMR, W, LAMI = LapackUtil.lag2(AA,EE,small)
+   if iszero(EE[2,1])
+      # handle the case when EE is upper triangular (required by lag2)
+      scale1, scale2, LAMR, W, LAMI = LapackUtil.lag2(AA,EE,small)
+   else
+      # handle the case when EE is full
+      E1 = copy(EE)
+      G, E1[1,1] = givens(E1[1,1],E1[2,1],1,2)
+      lmul!(G,view(E1,:,2)); E1[2,1] = ZERO
+      scale1, scale2, LAMR, W, LAMI = LapackUtil.lag2(lmul!(G,copy(AA)),E1,small)
+   end
    LAMI == ZERO && error("The pair (A,E) has real generalized eigenvalues")
-   # Compute right orthogonal transformation matrix Q.
+   # Compute right orthogonal transformation matrix Q (modified to cope with nonzero E[2,1])
    @inbounds CR, CI, SR, SI, L =  cgivensc2( scale1*AA[1,1] - EE[1,1]*LAMR,
-                                   -EE[1,1]*LAMI, scale1*AA[2,1], ZERO, small )
+                                   -EE[1,1]*LAMI, scale1*AA[2,1] - EE[2,1]*LAMR, -EE[2,1]*LAMI, small )
    QR = [ CR SR; -SR CR ]
    QI = [ -CI  -SI; -SI CI ]
    #  A := Q * A

src/riccati.jl

@@ -88,6 +88,8 @@ function arec(A::AbstractMatrix, G::Union{AbstractMatrix,UniformScaling,Real,Com
     eltype(G) == T || (typeof(G) <: AbstractMatrix ? G = convert(Matrix{T},G) : G = convert(T,G.λ)*I)
     eltype(Q) == T || (typeof(Q) <: AbstractMatrix ? Q = convert(Matrix{T},Q) : Q = convert(T,Q.λ)*I)
 
+    n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0) )
+
     S = schur([A  -G; -Q  -copy(A')])
 
     as ? select = real(S.values) .> 0 : select = real(S.values) .< 0
@@ -282,6 +284,9 @@ function arec(A::AbstractMatrix, B::AbstractVecOrMat, G::Union{AbstractMatrix,Un
          S = convert(Matrix{T},S)
       end
    end
+
+   n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0), zeros(T,m,0) )
+
    S0flag = iszero(S)
    typeof(R) <: UniformScaling && (R = Matrix{T}(R,m,m))
    SR = schur(R)
@@ -375,6 +380,9 @@ function garec(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, G::Un
     eident || eltype(E) == T || (E = convert(Matrix{T},E))
     eltype(G) == T || (typeof(G) <: AbstractMatrix ? G = convert(Matrix{T},G) : G = convert(T,G.λ)*I)
     eltype(Q) == T || (typeof(Q) <: AbstractMatrix ? Q = convert(Matrix{T},Q) : Q = convert(T,Q.λ)*I)
+
+    n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0) )
+
     if !eident
        Et = LinearAlgebra.LAPACK.getrf!(copy(E))
        LinearAlgebra.LAPACK.gecon!('1',Et[1],opnorm(E,1))  < epsm && error("E must be non-singular")
@@ -562,11 +570,11 @@ function garec(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, B::Ab
     eltype(Q) == T || (typeof(Q) <: AbstractMatrix ? Q = convert(Matrix{T},Q) : Q = convert(T,Q.λ)*I)
     eltype(R) == T || (typeof(R) <: AbstractMatrix ? R = convert(Matrix{T},R) : R = convert(T,R.λ)*I)
     eltype(S) == T || (typeof(S) <: AbstractVector ? S = convert(Vector{T},S) : S = convert(Matrix{T},S))
+
+    n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0), zeros(T,m,0) )
+
     cond(R)*epsm < 1 || error("R must be non-singular")
-    if !eident
-       Et = LinearAlgebra.LAPACK.getrf!(copy(E))
-       LinearAlgebra.LAPACK.gecon!('1',Et[1],opnorm(E,1))  < epsm && error("E must be non-singular")
-    end
+
     #  Method:  A stable/ant-stable deflating subspace Z1 = [Z11; Z21; Z31] of the pencil
     #               [  A  -G    B ]      [ E  0  0 ]
     #      L -s P = [ -Q  -A'  -S ]  - s [ 0  E' 0 ]
@@ -611,7 +619,7 @@ end
 
 
 """
-    ared(A, B, R, Q, S; as = false, rtol::Real = n) -> (X, EVALS, F, Z)
+    ared(A, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)
 
 Compute `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or 
 anti-stabilizing solution (if `as = true`) of the discrete-time algebraic Riccati equation
@@ -685,7 +693,7 @@ function ared(A::AbstractMatrix, B::AbstractVecOrMat, R::Union{AbstractMatrix,Un
     gared(A, I, B, R, Q, S; as = as, rtol = rtol)
 end
 """
-    gared(A, E, B, R, Q, S; as = false, rtol::Real = n) -> (X, EVALS, F, Z)
+    gared(A, E, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)
 
 Compute `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or 
 anti-stabilizing solution (if `as = true`) of the generalized discrete-time
@@ -807,6 +815,9 @@ function gared(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, B::Ab
     eltype(Q) == T || (typeof(Q) <: AbstractMatrix ? Q = convert(Matrix{T},Q) : Q = convert(T,Q.λ)*I)
     eltype(R) == T || (typeof(R) <: AbstractMatrix ? R = convert(Matrix{T},R) : R = convert(T,R.λ)*I)
     eltype(S) == T || (typeof(S) <: AbstractVector ? S = convert(Vector{T},S) : S = convert(Matrix{T},S))
+
+    n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0), zeros(T,m,0) )
+
     if !eident
       Et = LinearAlgebra.LAPACK.getrf!(copy(E))
       LinearAlgebra.LAPACK.gecon!('1',Et[1],opnorm(E,1))  < epsm && error("E must be non-singular")