Ensuring MKL compliance.

Project.toml

@@ -1,21 +1,22 @@
 name = "MatrixPencils"
 uuid = "48965c70-4690-11ea-1f13-43a2532b2fa8"
 authors = ["Andreas Varga <[email protected]>"]
-version = "1.7.11"
+version = "1.7.12"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [compat]
+LinearAlgebra = "1"
 Polynomials = "2.0, 3, 4"
+Random = "1"
 julia = "1.6"
 
 [extras]
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Random", "SparseArrays", "Test"]
+test = ["SparseArrays", "Test"]

README.md

@@ -34,7 +34,7 @@ The current version of the package includes the following functions:
 **Manipulation of general linear matrix pencils**
 
 * **pbalance!** Balancing linear matrix pencils.  
-* **pbalqual**  Balancing quality of a matrix pencils.  
+* **pbalqual**  Evaluation of the balancing quality of a matrix pencils.  
 * **preduceBF**  Reduction to the basic condensed form  `[B A-λE; D C]` with `E` upper triangular and nonsingular.
 * **klf**   Computation of the Kronecker-like form exhibiting the full Kronecker structure.
 * **klf_left**  Computation of the Kronecker-like form exhibiting the left and finite Kronecker structures.
@@ -84,6 +84,7 @@ The current version of the package includes the following functions:
 **Manipulation of linearizations of polynomial or rational matrices**
 
 * **lsbalance!** Scaling of a descriptor system based linearization. 
+* **lsbalqual** Evaluation of the scaling quality of descriptor system based linearizations.
 * **lsminreal** Computation of minimal order linearizations of the form `[A-λE B; C D]`.
 * **lsminreal2** Computation of minimal order linearizations of the form `[A-λE B; C D]` (potentially more efficient).
 * **lpsminreal**  Computation of strong minimal pencil based linearizations of the form `[A-λE B-λF; C-λG D-λH]`.
@@ -92,7 +93,6 @@ The current version of the package includes the following functions:
 * **lseval**   Evaluation of the value of the rational matrix corresponding to a descriptor system based linearization.
 * **lpseval**  Evaluation of the value of the rational matrix corresponding to a pencil based linearization.
 * **lps2ls**  Conversion of a pencil based linearization into a descriptor system based linearization.
-* **lsbalqual** Evaluation of the scaling quality of descriptor system based linearizations.
 
 **Manipulation of polynomial matrices** 
 

ReleaseNotes.md

@@ -1,5 +1,13 @@
 # Release Notes
 
+## Version 1.7.12
+Version bump to allow using both `OpenBLAS` and `MKL`.  
+The Julia  `_dlanv2` replaces (temporarily) the wrapper for `dlanv2` 
+(this will be changed back after[MKL#138](https://github.com/JuliaLinearAlgebra/MKL.jl/issues/138) is fixed)
+
+## Version 1.7.11
+Version bump to include `Polynomials.jl` v4.0. 
+
 ## Version 1.7.10
 This patch release includes the new function `pbalqual` and fixes some small bugs.
 

docs/src/index.md

@@ -36,7 +36,7 @@ The current version of the package includes the following functions:
 | Function | Description |
 | :--- | :--- |
 | **[`pbalance!`](@ref)** |Balancing arbitrary matrix pencils.  |
-| **[`pbalqual`](@ref)** |Balancing quality of a matrix pencils.  |
+| **[`pbalqual`](@ref)** |Evaluation of the balancing quality of a matrix pencils.  |
 | **[`preduceBF`](@ref)** | Reduction to the basic condensed form  `[B A-λE; D C]` with `E` upper triangular and nonsingular. |
 | **[`klf`](@ref)** |   Computation of the Kronecker-like form exhibiting the full Kronecker structure |
 | **[`klf_right`](@ref)** |   Computation of the Kronecker-like form exhibiting the right and finite Kronecker structures |
@@ -100,6 +100,7 @@ The current version of the package includes the following functions:
 | Function | Description |
 | :--- | :--- |
 | **[`lsbalance!`](@ref)** | Scaling of a descriptor system based linearization |
+| **[`lsbalqual`](@ref)** | Evaluation of the scaling quality of descriptor system based linearizations |
 | **[`lsminreal`](@ref)** | Computation of irreducible descriptor system based linearizations  |
 | **[`lsminreal2`](@ref)** | Computation of irreducible descriptor system based linearizations (potentially more efficient)|
 | **[`lpsminreal`](@ref)** | Computation of strongly minimal pencil based linearizations |
@@ -108,8 +109,6 @@ The current version of the package includes the following functions:
 | **[`lseval`](@ref)** | Evaluation of the value of the rational matrix corresponding to a descriptor system based linearization  |
 | **[`lpseval`](@ref)** | Evaluation of the value of the rational matrix corresponding to a pencil based linearization |
 | **[`lps2ls`](@ref)** | Conversion of a pencil based linearization into a descriptor system based linearization |
-| **[`lsbalqual`](@ref)** | Evaluation of the scaling quality of descriptor system based linearizations |
-
 
 **Manipulation of polynomial matrices** 
 

src/gsfstab.jl

@@ -1045,14 +1045,96 @@ function ordeigvals(A::AbstractMatrix{T}) where T
                αi[i+1] = -ei
             else
                # an arbitrary 2x2 block
-               αr[i], αi[i], αr[i+1], αi[i+1] = lanv2(A[i,i], A[i,i+1], A[i+1,i], A[i+1,i+1])
+               # αr[i], αi[i], αr[i+1], αi[i+1] = lanv2(A[i,i], A[i,i+1], A[i+1,i], A[i+1,i+1])
+               αr[i], αi[i], αr[i+1], αi[i+1] = _lanv2(A[i,i], A[i,i+1], A[i+1,i], A[i+1,i+1])
             end
             i +=1
          end
       end
       return iszero(αi) ? αr : Complex.(αr,αi)
    end
 end
+function _lanv2(a::T,b::T,c::T,d::T) where {T <: Real}
+   # compute the eigenvalues of a real 2x2 [a b; c d] in e1r, e1i, e2r, e2i 
+   # extracted from the translation from LAPACK::dlanv2 by Ralph Smith
+   # Copyright:
+   # Univ. of Tennessee
+   # Univ. of California Berkeley
+   # Univ. of Colorado Denver
+   # NAG Ltd.
+   ZERO = zero(T)
+   ONE = one(T)
+   sgn(x) = (x < 0) ? -ONE : ONE # fortran sign differs from Julia
+   half = ONE / 2
+   small = 4eps(T) # how big discriminant must be for easy reality check
+   if c==0
+   elseif b==0
+       # swap rows/cols
+       a,b,c,d = d,-c,ZERO,a
+   elseif ((a-d) == 0) && (b*c < 0)
+       # nothing to do
+   else
+       asubd = a-d
+       p = half*asubd
+       bcmax = max(abs(b),abs(c))
+       bcmis = min(abs(b),abs(c)) * sgn(b) * sgn(c)
+       scale = max(abs(p), bcmax)
+       z = (p / scale) * p + (bcmax / scale) * bcmis
+       # if z is of order machine accuracy: postpone decision
+       if z >= small
+           # real eigenvalues
+           z = p + sqrt(scale) * sqrt(z) * sgn(p)
+           a = d + z
+           d -= (bcmax / z) * bcmis
+           b -= c
+           c = ZERO
+       else
+           # complex or almost equal real eigenvalues
+           σ = b + c
+           τ = hypot(σ, asubd)
+           cs = sqrt(half * (ONE + abs(σ) / τ))
+           sn = -(p / (τ * cs)) * sgn(σ)
+           # apply rotations
+           aa = a*cs + b*sn
+           bb = -a*sn + b*cs
+           cc = c*cs + d*sn
+           dd = -c*sn + d*cs
+           a = aa*cs + cc*sn
+           b = bb*cs + dd*sn
+           c = -aa*sn + cc*cs
+           d = -bb*sn + dd*cs
+           midad = half * (a+d)
+           a = midad
+           d = a
+           if (c != 0)
+               if (b != 0)
+                   if b*c >= 0
+                       # real eigenvalues
+                       sab = sqrt(abs(b))
+                       sac = sqrt(abs(c))
+                       p = sab*sac*sgn(c)
+                       a = midad + p
+                       d = midad - p
+                       b -= c
+                       c = 0
+                   end
+               else
+                   b,c = -c,ZERO
+               end
+           end
+       end
+   end
+
+   w1r,w2r = a, d
+   if c==0
+       w1i,w2i = ZERO,ZERO
+   else
+       rti = sqrt(abs(b))*sqrt(abs(c))
+       w1i,w2i = rti,-rti
+   end
+   return w1r,w1i,w2r,w2i
+end
+
 """
     ordeigvals(A, B) -> (ev, α, β)
 

src/lapackutil2.jl

@@ -1,17 +1,30 @@
 module LapackUtil2
 
 #const liblapack = VERSION < v"1.7" ? Base.liblapack_name : "libblastrampoline" * (Sys.iswindows() ? "-5" : "")
-const liblapack = Base.liblapack_name
 
 import LinearAlgebra.BLAS.@blasfunc
 
-import LinearAlgebra: BlasFloat, BlasReal, BlasComplex, BlasInt, LAPACKException, 
+using LinearAlgebra
+import LinearAlgebra: libblastrampoline, BlasFloat, BlasReal, BlasComplex, BlasInt, LAPACKException, 
        DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare
 
 using Base: iszero, has_offset_axes
+using LinearAlgebra.LAPACK
 
 export larfg!, larfgl!, larf!
 
+#const liblapack = Base.liblapack_name
+
+@static if VERSION < v"1.7"
+    using LinearAlgebra.LAPACK: liblapack
+elseif VERSION < v"1.9"
+    const liblapack = "libblastrampoline"
+else
+    const liblapack = LinearAlgebra.libblastrampoline
+end
+
+@show liblapack
+
 function chkside(side::AbstractChar)
     # Check that left/right hand side multiply is correctly specified
     side == 'L' || side == 'R' ||
@@ -25,6 +38,7 @@ end
 
 for (fn, elty) in ((:dlanv2_, :Float64),
                    (:slanv2_, :Float32))
+
     @eval begin
         function lanv2(A::$elty, B::$elty, C::$elty, D::$elty)
            """