diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index f059577..038b91e 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -20,7 +20,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - "1.0"  # LTS
+          - "1.6"  # first supported version
           - "1"    # Latest Release
         os:
           - ubuntu-latest
diff --git a/Project.toml b/Project.toml
index c88423e..3bb255c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,10 +7,11 @@ version = "0.1.0"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
-julia = "1"
+julia = "1.6"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 
 [targets]
-test = ["Test"]
+test = ["Test","Random"]
diff --git a/README.md b/README.md
index 70338f0..292878b 100644
--- a/README.md
+++ b/README.md
@@ -2,35 +2,43 @@
 
 WARNING: Package still in development!
 
-This package proposes specialized functions for dense real skew-symmetric matrices i.e $A=-A^T$.
-It provides the structure [`SkewSymmetric`](@ref) and the classical linear operations on such 
+This package provides specialized functions for dense real skew-symmetric matrices i.e $A=-A^T$.
+It provides the type [`SkewHermitian`](@ref) and the classical linear operations on such
 matrices. The package fits in the framework given by the LinearAlgebra package.
 
-In particular, the package provides de following functions for $A::$ [`SkewSymmetric`](@ref) :
+In particular, the package provides de following functions for $A::$ [`SkewHermitian`](@ref) :
 
 -Tridiagonal reduction: [`hessenberg`](@ref)\
 -Eigensolvers: [`eigen`](@ref), [`eigvals`](@ref),[`eigmax`](@ref),[`eigmin`](@ref)\
 -SVD: [`svd`](@ref), [`svdvals`](@ref)\
 -Trigonometric functions:[`exp`](@ref), [`cis`](@ref),[`cos`](@ref),[`sin`](@ref),[`tan`](@ref),[`sinh`](@ref),[`cosh`](@ref),[`tanh`](@ref)
 
-The SkewSymmetric type uses the complete matrix representation as data. It doesn't verify automatically that the given matrix input is skew-symmetric. In particular, in-place methods could destroy the skew-symmetry. The provided function [`isskewsymmetric(A)`](@ref) allows to verify that A is indeed skew-symmetric.
-Here is a basic example to initialize a [`SkewSymmetric`](@ref)
-```
+The `SkewHermitian(A)` wraps an existing matrix `A`, which *must* already be skew-Hermitian,
+in the `SkewHermitian` type, which supports fast specialized operations noted above.  You
+can use the function `isskewhermitian(A)` to check whether `A` is skew-Hermitian (`A == -A'`).
+
+Alternatively, you can use the funcition `skewhermitian(A)` to take the skew-Hermitian *part*
+of `A`, given by `(A - A')/2`, and wrap it in a `SkewHermitian` view.  Alternatively, the
+function `skewhermitian!(A)` does the same operation in-place on `A`.
+
+Here is a basic example to initialize a [`SkewHermitian`](@ref)
+```jl
 julia> A = [0 2 -7 4; -2 0 -8 3; 7 8 0 1;-4 -3 -1 0]
 3×3 Matrix{Int64}:
   0  2 -7  4
  -2  0 -8  3
   7  8  0  1
   -4 -3 -1 0
-julia> A = SkewSymmetric(A)
-4×4 Matrix{Int64}:
+
+julia> isskewhermitian(A)
+true
+
+julia> A = SkewHermitian(A)
+4×4 SkewHermitian{Int64, Matrix{Int64}}:
   0   2  -7  4
  -2   0  -8  3
   7   8   0  1
  -4  -3  -1  0
- 
-julia> isskewsymmetric(A)
-true
 
 julia> tr(A)
 0
@@ -39,12 +47,12 @@ julia> det(A)
 81.0
 
 julia> inv(A)
-4×4 Matrix{Float64}:
- -0.0        0.111111  -0.333333     -0.888889
- -0.111111   0.0        0.444444      0.777778
-  0.333333  -0.444444   2.77556e-17   0.222222
-  0.888889  -0.777778  -0.222222      0.0
-  
+4×4 SkewHermitian{Float64, Matrix{Float64}}:
+  0.0        0.111111  -0.333333  -0.888889
+ -0.111111   0.0        0.444444   0.777778
+  0.333333  -0.444444   0.0        0.222222
+  0.888889  -0.777778  -0.222222   0.0
+
 julia> x=[1;2;3;4]
 4-element Vector{Int64}:
  1
@@ -59,7 +67,7 @@ julia> A\x
   0.3333333333333336
  -1.3333333333333333
   ```
-  
+
   The functions from the LinearAlgebra package can be used in the same fashion:
   ```
 julia> hessenberg(A)
@@ -68,22 +76,22 @@ julia> hessenberg(A)
  8.30662   0.0      -8.53382   ⋅
   ⋅        8.53382   0.0      1.08347
   ⋅         ⋅       -1.08347  0.0
-  
+
  julia> eigvals(A)
 4-element Vector{ComplexF64}:
   0.0 + 11.93445871397423im
   0.0 + 0.7541188264752853im
  -0.0 - 0.7541188264752877im
  -0.0 - 11.934458713974225im
- 
+
  julia> eigmax(A)
 -0.0 - 11.934458713974223im
-  
+
   ```
   \usepackage{amsymb}\
  ## Hessenberg/tridiagonal reduction
 The hessenberg reduction performs a reduction $A=QHQ^T$ where $Q=\prod_i I-\tau_i v_iv_i^T$ is an orthonormal matrix.
-The [`hessenberg`](@ref) function returns a structure of type [`SkewHessenberg`](@ref) containing the [`Tridiagonal`](@ref) reduction $H\in \mathbb{R}^{n\times n}$, the householder reflectors $v_i$ in a  [`UnitLowerTriangular`](@ref) $V\in \mathbb{R}^{n-1\times n-1}$  and the $n-2$ scalars $\tau_i$ associated to the reflectors. A function [`getQ`](@ref) is provided to retrieve the orthogonal transformation Q. 
+The [`hessenberg`](@ref) function returns a structure of type [`SkewHessenberg`](@ref) containing the [`Tridiagonal`](@ref) reduction $H\in \mathbb{R}^{n\times n}$, the householder reflectors $v_i$ in a  [`UnitLowerTriangular`](@ref) $V\in \mathbb{R}^{n-1\times n-1}$  and the $n-2$ scalars $\tau_i$ associated to the reflectors. A function [`getQ`](@ref) is provided to retrieve the orthogonal transformation Q.
 
   ```
   julia> H=hessenberg(A)
@@ -93,31 +101,31 @@ The [`hessenberg`](@ref) function returns a structure of type [`SkewHessenberg`]
   ⋅        8.53382   0.0      1.08347
   ⋅         ⋅       -1.08347  0.0
 3×3 UnitLowerTriangular{Float64, SubArray{Float64, 2, Matrix{Float64}, Tuple{UnitRange{Int64}, UnitRange{Int64}}, false}}:
-  1.0         ⋅         ⋅ 
+  1.0         ⋅         ⋅
  -0.679175   1.0        ⋅
   0.3881    -0.185238  1.0
 2-element Vector{Float64}:
  1.2407717061715382
  1.9336503566410876
- 
+
  julia> Q=getQ(H)
 4×4 Matrix{Float64}:
  1.0   0.0        0.0         0.0
  0.0  -0.240772  -0.95927    -0.14775
  0.0   0.842701  -0.282138    0.458534
  0.0  -0.481543  -0.0141069   0.876309
- 
+
  ```
 \usepackage{hyperref}
 
  ## Eigenvalues and eigenvectors
- 
- The package also provides eigensolvers for  [`SkewSymmetric`](@ref) matrices. The method to solve the eigenvalue problem is issued from: \url{C.Penke,A.Marek,C.Vorwerk,C.Draxl,P.Benner, \textit{ High Performance Solution of Skew-symmetric Eigenvalue Problems with Applications in Solving Bethe-Salpeter Eigenvalue Problem}, 2020.}
- 
+
+ The package also provides eigensolvers for  [`SkewHermitian`](@ref) matrices. The method to solve the eigenvalue problem is issued from: \url{C.Penke,A.Marek,C.Vorwerk,C.Draxl,P.Benner, \textit{ High Performance Solution of Skew-symmetric Eigenvalue Problems with Applications in Solving Bethe-Salpeter Eigenvalue Problem}, 2020.}
+
   The function [`eigen`](@ref) returns the eigenvalues plus the real part of the eigenvectors and the imaginary part separeted.
  ```
   julia> val,Qr,Qim=eigen(A)
- 
+
   julia> val
 4-element Vector{ComplexF64}:
   0.0 + 11.934458713974193im
@@ -138,7 +146,7 @@ julia> Qim
  -0.176712   0.0931315   0.0931315  -0.176712
   0.615785  -0.284619   -0.284619    0.615785
  -0.299303  -0.640561   -0.640561   -0.299303
- 
+
  ```
  The function [`eigvals`](@ref) provides de eigenvalues of $A$. The eigenvalues can be sorted and found partially with imaginary part in some given real range or by order.
  ```
@@ -153,7 +161,7 @@ julia> eigvals(A,0,15)
 2-element Vector{ComplexF64}:
  0.0 + 11.93445871397414im
  0.0 + 0.7541188264752858im
- 
+
 julia> eigvals(A,1:3)
 3-element Vector{ComplexF64}:
   0.0 + 11.93445871397423im
@@ -161,10 +169,10 @@ julia> eigvals(A,1:3)
  -0.0 - 0.7541188264752758im
  ```
  ## SVD
- 
+
  Specialized SVD using the eigenvalue decomposition is implemented for [`SkewSymmetric`](@ref) type.
  These functions can be used using the LinearAlgebra syntax.
- 
+
  ```
  julia> svd(A)
 SVD{ComplexF64, Float64, Matrix{ComplexF64}}
@@ -186,17 +194,17 @@ Vt factor:
  0.0-0.49111im    -0.176712-0.488014im   0.615785-0.143534im  -0.299303-0.00717668im
  0.0-0.508735im  -0.0931315+0.471107im   0.284619+0.138561im   0.640561+0.00692804im
  0.0-0.508735im   0.0931315+0.471107im  -0.284619+0.138561im  -0.640561+0.00692804im
- 
+
  julia> svdvals(A)
 4-element Vector{Float64}:
  11.93445871397423
  11.934458713974225
   0.7541188264752877
   0.7541188264752853
- 
+
  ```
  ## Trigonometric functions
- 
+
  The package implements special versions of the trigonometric functions using the eigenvalue decomposition. The provided functions are [`exp`](@ref), [`cis`](@ref),[`cos`](@ref),[`sin`](@ref),[`tan`](@ref),[`sinh`](@ref),[`cosh`](@ref),[`tanh`](@ref).
  ```
  julia> exp(A)
@@ -205,7 +213,7 @@ Vt factor:
  -0.697298   0.140338     0.677464   0.187414
   0.578289  -0.00844255   0.40033    0.710807
   0.279941  -0.559925     0.555524  -0.547275
-  
+
  julia> cis(A)
 4×4 Matrix{ComplexF64}:
    5.95183+0.0im       3.21734+1.80074im     -0.658082-3.53498im      -1.4454+5.61775im
@@ -226,9 +234,9 @@ julia> cosh(A)
  -0.756913   0.140338  0.334511  -0.186256
   0.154821   0.334511  0.40033    0.633165
   0.340045  -0.186256  0.633165  -0.547275
- 
+
  ```
- 
-  
-  
-  
+
+
+
+
diff --git a/src/SkewLinearAlgebra.jl b/src/SkewLinearAlgebra.jl
index 0f704a6..4ce0277 100644
--- a/src/SkewLinearAlgebra.jl
+++ b/src/SkewLinearAlgebra.jl
@@ -5,61 +5,51 @@ and types for skew-symmetricmatrices, i.e A=-A^T
 """
 module SkewLinearAlgebra
 
-
-import LinearAlgebra as LA
 using LinearAlgebra
-import Base: \, /, *, ^, +, -, ==, copy,copyto!, size,setindex!,getindex,display,conj,conj!,similar,
-        isreal,cos,sin,cosh,sinh,tanh,cis
-export 
+import LinearAlgebra as LA
+export
     #Types
-    SkewSymmetric,
+    SkewHermitian,
     SkewHessenberg,
     #functions
-    isskewsymmetric,
+    isskewhermitian,
+    skewhermitian,
+    skewhermitian!,
     getQ
 
-struct SkewSymmetric{T<:Real,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}
+struct SkewHermitian{T<:Number,S<:AbstractMatrix{<:T}} <: AbstractMatrix{T}
     data::S
 
-    function SkewSymmetric{T,S}(data) where {T,S<:AbstractMatrix{<:T}}
+    function SkewHermitian{T,S}(data) where {T,S<:AbstractMatrix{<:T}}
         LA.require_one_based_indexing(data)
-        new{T,S}(data)  
+        new{T,S}(data)
     end
 end
 
 """
-    SkewSymmetric(A)
-Transform matrix A in a Skewsymmetric structure. A is assumed to be correctly 
-build as a skew-symmetric matrix. 'isskewsymmetric(A)' allows to verify skew-symmetry
-"""
+    SkewHermitian(A) <: AbstractMatrix
 
-function SkewSymmetric(A::AbstractMatrix)
-    LA.checksquare(A)
-    n=size(A,1)
-    n>1 || throw("Skew-symmetric cannot be of size  1x1")
-    return skewsymmetric_type(typeof(A))(A)
-end
-
-skewsymmetric(A::AbstractMatrix) = SkewSymmetric(A)
-skewsymmetric(A::Number) = throw("Number cannot be skewsymmetric")
+Construct a `SkewHermitian` view of the skew-Hermitian matrix `A` (`A == -A'`),
+which allows one to exploit efficient operations for eigenvalues, exponentiation,
+and more.
 
+Takes "ownership" of the matrix `A`.  See also [`skewhermitian`](@ref), which takes the
+skew-hermitian part of `A`, and [`skewhermitian!`](@ref), which does this in-place,
+along with [`isskewhermitian`](@ref) which checks whether `A == -A'`.
 """
-    skewsymmetric_type(T::Type)
-The type of the object returned by `skewsymmetric(::T, ::Symbol)`. For matrices, this is an
-appropriately typed `SkewSymmetric`. If `skewsymmetric` is
-implemented for a custom type, so should be `skewsymmetric_type`, and vice versa.
+function SkewHermitian(A::AbstractMatrix)
+    isskewhermitian(A) || throw(ArgumentError("matrix `A` must be skew-Hermitian (equal `-A')"))
+    return SkewHermitian{eltype(A),typeof(A)}(A)
+end
+
 """
+    skewhermitian(A)
 
-function skewsymmetric_type(::Type{T}) where {S, T<:AbstractMatrix{S}}
-    return SkewSymmetric{Union{S, promote_op(transpose, S), skewsymmetric_type(S)}, T}
-end
-function skewsymmetric_type(::Type{T}) where {S<:Number, T<:AbstractMatrix{S}}
-    return SkewSymmetric{S, T}
-end
-function skewsymmetric_type(::Type{T}) where {S<:AbstractMatrix, T<:AbstractMatrix{S}}
-    return SkewSymmetric{AbstractMatrix, T}
-end
-skewsymmetric_type(::Type{T}) where {T<:Number} = throw("Number cannot be skewsymmetric")
+Returns the skew-Hermitian part of A, i.e. `(A-A')/2`.  See also
+[`skewhermitian!`](@ref), which does this in-place.
+"""
+skewhermitian(A::AbstractMatrix) = skewhermitian!(Base.copymutable(A))
+skewhermitian(a::Number) = imag(a)
 
 """
     getindex(A,i,j)
@@ -67,206 +57,189 @@ skewsymmetric_type(::Type{T}) where {T<:Number} = throw("Number cannot be skewsy
 Returns the value A(i,j)
 """
 
-@inline function Base.getindex(A::SkewSymmetric, i::Integer, j::Integer)
-    @boundscheck checkbounds(A, i, j)
-    return @inbounds A.data[i,j]
-end 
+Base.@propagate_inbounds Base.getindex(A::SkewHermitian, i::Integer, j::Integer) = A.data[i,j]
 
-"""
-    setindex!(A,v,i,j)
-Set A(i,j)=v and A(j,i)=-v to conserve skew-symmetry
-"""
-
-function Base.setindex!(A::SkewSymmetric, v, i::Integer, j::Integer)
-    i!=j || throw("Cannot modify zero diagonal element")
-    setindex!(A.data, v, i, j)
-    setindex!(A.data, -v, j, i)
+Base.@propagate_inbounds function Base.setindex!(A::SkewHermitian, v, i::Integer, j::Integer)
+    if i == j
+        real(v) == 0 || throw(ArgumentError("diagonal elements must be zero"))
+    else
+        A.data[i,j] = v
+        A.data[j,i] = -v'
+    end
+    return v
 end
 
-Base.similar(A::SkewSymmetric, ::Type{T}) where {T} = SkewSymmetric(LA.similar(parent(A), T))
-Base.similar(A::SkewSymmetric) = SkewSymmetric(similar(parent(A)))
+Base.similar(A::SkewHermitian, ::Type{T}) where {T} = SkewHermitian(LA.similar(parent(A), T) .= 0)
+Base.similar(A::SkewHermitian) = SkewHermitian(similar(parent(A)) .= 0)
 
 # Conversion
-function Matrix(A::SkewSymmetric)
-    B = copy(A.data)
-    return B
-end
-Base.Array(A::SkewSymmetric) = convert(Matrix, A)
+Base.Matrix(A::SkewHermitian) = Matrix(A.data)
+Base.Array(A::SkewHermitian) = Matrix(A)
 
-Base.parent(A::SkewSymmetric) = A.data
-SkewSymmetric{T,S}(A::SkewSymmetric{T,S}) where {T,S<:AbstractMatrix{T}} = A
-SkewSymmetric{T,S}(A::SkewSymmetric) where {T,S<:AbstractMatrix{T}} = SkewSymmetric{T,S}(convert(S,A.data))
-Base.AbstractMatrix{T}(A::SkewSymmetric) where {T} = SkewSymmetric(convert(AbstractMatrix{T}, A.data))
+Base.parent(A::SkewHermitian) = A.data
+SkewHermitian{T,S}(A::SkewHermitian{T,S}) where {T,S} = A
+SkewHermitian{T,S}(A::SkewHermitian) where {T,S<:AbstractMatrix{T}} = SkewHermitian{T,S}(S(A.data))
+Base.AbstractMatrix{T}(A::SkewHermitian) where {T} = SkewHermitian(AbstractMatrix{T}(A.data))
 
-Base.copy(A::SkewSymmetric{T,S}) where {T,S} = (B = copy(A.data); SkewSymmetric{T,typeof(B)}(B))
-
-function copyto!(dest::SkewSymmetric, src::SkewSymmetric)
+Base.copy(A::SkewHermitian) = SkewHermitian(copy(A.data))
+function Base.copyto!(dest::SkewHermitian, src::SkewHermitian)
     copyto!(dest.data, src.data)
     return dest
 end
-Base.size(A::SkewSymmetric,n) = size(A.data,n)
-Base.size(A::SkewSymmetric) = size(A.data)
+function Base.copyto!(dest::SkewHermitian, src::AbstractMatrix)
+    isskewhermitian(src) || throw(ArgumentError("can only copy skew-Hermitian data to SkewHermitian"))
+    copyto!(dest.data, src)
+    return dest
+end
+Base.copyto!(dest::AbstractMatrix, src::SkewHermitian) = copyto!(dest, src.data)
 
+Base.size(A::SkewHermitian,n) = size(A.data,n)
+Base.size(A::SkewHermitian) = size(A.data)
 
 """
-    isskewsymmetric(A)
+    isskewhermitian(A)
 
-Verifies skew-symmetry of a matrix
+Returns whether `A` is skew-Hermitian, i.e. whether `A == -A'`.
 """
-
-function isskewsymmetric(A::SkewSymmetric) 
-    n=size(A,1)
-    for i=1:n
-        getindex(A,i,i) == 0 || return false
-        for j=1:i-1
-            getindex(A,i,j) == -getindex(A,j,i) ||return false
+function isskewhermitian(A::AbstractMatrix{<:Number})
+    axes(A,1) == axes(A,2) || throw(ArgumentError("axes $(axes(A,1)) and $(axex(A,2)) do not match"))
+    @inbounds for i in axes(A,1)
+        for j = firstindex(A, 1):i
+            A[i,j] == -A[j,i]' || return false
         end
     end
     return true
 end
+isskewhermitian(A::SkewHermitian) = true
+isskewhermitian(a::Number) = a == -a'
+
+"""
+    skewhermitian(A)
+
+Transforms `A` in-place to its skew-Hermitian part `(A-A')/2`,
+and returns a [`SkewHermitian`](@ref) view.
+"""
+function skewhermitian!(A::AbstractMatrix{<:Number})
+    LA.require_one_based_indexing(A)
+    n = LA.checksquare(A)
+    @inbounds for i in 1:n
+        A[i,i] = imag(A[i,i])
+        for j = 1:i-1
+            a = (A[i,j] - A[j,i]')/2
+            A[i,j] = a
+            A[j,i] = -a'
+        end
+    end
+    return SkewHermitian(A)
+end
 
 #Classic operators on a matrix
-Base.isreal(A::SkewSymmetric)=true
-Base.transpose(A::SkewSymmetric) = SkewSymmetric(-A.data)
-Base.adjoint(A::SkewSymmetric{<:Real}) = SkewSymmetric(-A.data)
-Base.adjoint(A::SkewSymmetric) = SkewSymmetric(-A.data)
-Base.real(A::SkewSymmetric{<:Real}) = A
-Base.real(A::SkewSymmetric) = A
-Base.imag(A::SkewSymmetric) = SkewSymmetric(LA.imag(A.data))
-
-Base.copy(A::SkewSymmetric) =SkewSymmetric(copy(parent(A)))
-Base.display(A::SkewSymmetric) = display(A.data)
-Base.conj(A::SkewSymmetric) = typeof(A)(A.data)
-Base.conj!(A::SkewSymmetric) = typeof(A)(A.data)
-LA.tr(A::SkewSymmetric) = 0
-
-
-LA.tril!(A::SkewSymmetric) = tril!(A.data)
-LA.tril(A::SkewSymmetric)  = tril(A.data)
-LA.triu!(A::SkewSymmetric) = triu!(A.data)
-LA.triu(A::SkewSymmetric)  = triu(A.data)
-LA.tril!(A::SkewSymmetric,k::Integer) = tril!(A.data,k)
-LA.tril(A::SkewSymmetric,k::Integer)  = tril(A.data,k)
-LA.triu!(A::SkewSymmetric,k::Integer) = triu!(A.data,k)
-LA.triu(A::SkewSymmetric,k::Integer)  = triu(A.data,k)
-
-
-function LA.dot(A::SkewSymmetric, B::SkewSymmetric)
+Base.isreal(A::SkewHermitian) = isreal(A.data)
+Base.transpose(A::SkewHermitian) = SkewHermitian(transpose(A.data))
+Base.adjoint(A::SkewHermitian) = SkewHermitian(A.data')
+Base.real(A::SkewHermitian{<:Real}) = A
+Base.real(A::SkewHermitian) = SkewHermitian(real(A.data))
+Base.imag(A::SkewHermitian) = SkewHermitian(imag(A.data))
+
+Base.copy(A::SkewHermitian) = SkewHermitian(copy(A))
+Base.conj(A::SkewHermitian) = SkewHermitian(conj(A.data))
+Base.conj!(A::SkewHermitian) = SkewHermitian(conj!(A.data))
+LA.tr(A::SkewHermitian{<:Real}) = zero(eltype(A))
+LA.tr(A::SkewHermitian) = tr(A.data)
+
+LA.tril!(A::SkewHermitian) = tril!(A.data)
+LA.tril(A::SkewHermitian)  = tril(A.data)
+LA.triu!(A::SkewHermitian) = triu!(A.data)
+LA.triu(A::SkewHermitian)  = triu(A.data)
+LA.tril!(A::SkewHermitian,k::Integer) = tril!(A.data,k)
+LA.tril(A::SkewHermitian,k::Integer)  = tril(A.data,k)
+LA.triu!(A::SkewHermitian,k::Integer) = triu!(A.data,k)
+LA.triu(A::SkewHermitian,k::Integer)  = triu(A.data,k)
+
+function LA.dot(A::SkewHermitian, B::SkewHermitian)
     n = size(A, 2)
     if n != size(B, 2)
-        throw(DimensionMismatch("A has dimensions $(size(A)) but B has dimensions $(size(B))"))
+        throw(DimensionMismatch("A has size $(size(A)) but B has size $(size(B))"))
     end
     dotprod = zero(dot(first(A), first(B)))
-    @inbounds for j = 1:n
-        for i = 1:(j - 1)
-            dotprod += 2 *(dot(A.data[i, j], B.data[i, j]))
-        end
+    @inbounds for j = 1:n, i = 1:j-1
+        dotprod += 2 * real(dot(A.data[i, j], B.data[i, j]))
     end
-    
     return dotprod
 end
 
-Base. -(A::SkewSymmetric) = SkewSymmetric(- A.data)
-
+Base.:-(A::SkewHermitian) = SkewHermitian(-A.data)
 
 for f in (:+, :-)
     @eval begin
-        $f(A::SkewSymmetric, B::SkewSymmetric) = SkewSymmetric($f(A.data, B.data))
+        Base.$f(A::SkewHermitian, B::SkewHermitian) = SkewHermitian($f(A.data, B.data))
    end
 end
 
 ## Matvec
-@inline function LA.mul!(y::StridedVector{T}, A::SkewSymmetric{T,<:StridedMatrix}, x::StridedVector{T},
-            α::Number, β::Number) where {T<:LA.BlasFloat}
-    alpha, beta = promote(α, β, zero(T))
-    if alpha isa Union{Bool,T} && beta isa Union{Bool,T}
-        return LA.BLAS.gemv!('N', alpha, A.data, x, beta, y)
-    else
-        return generic_matvecmul!(y, 'N', A, x, MulAddMul(α, β))
+LA.mul!(y::StridedVecOrMat, A::SkewHermitian, x::StridedVecOrMat, α::Number, β::Number) =
+    LA.mul!(y, A.data, x, α, β)
+LA.mul!(y::StridedVecOrMat, A::SkewHermitian, x::StridedVecOrMat) =
+    LA.mul!(y, A.data, x)
+LA.mul!(y::StridedVecOrMat, A::SkewHermitian, B::SkewHermitian, α::Number, β::Number) =
+    LA.mul!(y, A.data, B.data, α, β)
+LA.mul!(y::StridedVecOrMat, A::SkewHermitian, B::SkewHermitian) =
+    LA.mul!(y, A.data, B.data)
+LA.mul!(y::StridedVecOrMat, A::StridedMatrix, B::SkewHermitian, α::Number, β::Number) =
+    LA.mul!(y, A, B.data, α, β)
+LA.mul!(y::StridedVecOrMat, A::StridedMatrix, B::SkewHermitian) =
+    LA.mul!(y, A, B.data)
+
+function LA.dot(x::AbstractVector, A::SkewHermitian, y::AbstractVector)
+    LA.require_one_based_indexing(x, y)
+    (length(x) == length(y) == size(A, 1)) || throw(DimensionMismatch())
+    r = zero(eltype(x)) * zero(eltype(A)) * zero(eltype(y))
+    @inbounds for j = 1:length(y)
+        r += dot(x[j], A.data[j,j], y[j]) # zero if A is real
+        @simd for i = 1:j-1
+            Aij = A.data[i,j]
+            r += dot(x[i], Aij, y[j]) + dot(x[j], -Aij', y[i])
+        end
     end
+    return r
 end
-## Matmat
-@inline function LA.mul!(C::StridedMatrix{T}, A::SkewSymmetric{T,<:StridedMatrix}, B::StridedMatrix{T},
-            α::Number, β::Number) where {T<:LA.BlasFloat}
-    alpha, beta = promote(α, β, zero(T))
-    if alpha isa Union{Bool,T} && beta isa Union{Bool,T}
-        return LA.BLAS.gemm!('N', 'N', alpha, A.data, B, beta, C)
-    else
-        return generic_matmatmul!(C, 'N', 'N', A, B, MulAddMul(alpha, beta))
+
+# Scaling:
+for op in (:*, :/, :\)
+    if op in (:*, :/)
+        @eval Base.$op(A::SkewHermitian, x::Number) = $op(A.data, x)
+        @eval Base.$op(A::SkewHermitian, x::Real) = SkewHermitian($op(A.data, x))
     end
-end
-@inline function LA.mul!(C::StridedMatrix{T}, A::StridedMatrix{T}, B::SkewSymmetric{T,<:StridedMatrix},
-        α::Number, β::Number) where {T<:LA.BlasFloat}
-    alpha, beta = promote(α, β, zero(T))
-    if alpha isa Union{Bool,T} && beta isa Union{Bool,T}
-        return LA.BLAS.gemm!('N', 'N', alpha, B.data, A, beta, C)
-    else
-        return generic_matmatmul!(C, 'N', 'N', A, B, MulAddMul(alpha, beta))
+    if op in (:*, :\)
+        @eval Base.$op(x::Number, A::SkewHermitian) = $op(x, A.data)
+        @eval Base.$op(x::Real, A::SkewHermitian) = SkewHermitian($op(x, A.data))
     end
 end
-@inline function LA.mul!(C::StridedMatrix{T}, A::SkewSymmetric{T,<:StridedMatrix}, B::SkewSymmetric{T,<:StridedMatrix},
-        α::Number, β::Number) where {T<:LA.BlasFloat}
-    alpha, beta = promote(α, β, zero(T))
-    if alpha isa Union{Bool,T} && beta isa Union{Bool,T}
-        return LA.BLAS.gemm!('N', 'N', alpha, A.data, B.data, beta, C)
-    else
-        return generic_matmatmul!(C, 'N', 'N', A, B, MulAddMul(alpha, beta))
-    end
+function checkreal(x::Number)
+    isreal(x) || throw(ArgumentError("in-place scaling of SkewHermitian requires a real scalar"))
+    return real(a)
 end
+LA.rdiv!(A::SkewHermitian, b::Number) = LA.rdiv!(A.data, checkreal(b))
+LA.ldiv!(b::Number, A::SkewHermitian) = LA.ldiv!(checkreal(b), A.data)
+LA.rmul!(A::SkewHermitian, b::Number) = LA.rmul!(A.data, checkreal(b))
+LA.lmul!(b::Number, A::SkewHermitian) = LA.lmul!(checkreal(b), A.data)
 
+for f in (:det, :logdet, :lu, :lu!, :lq, :lq!, :qr, :qr!)
+    @eval LA.$f(A::SkewHermitian) = LA.$f(A.data)
+end
+@eval LA.inv(A::SkewHermitian) = skewhermitian!(LA.inv(A.data))
 
-Base. *(A::SkewSymmetric, B::SkewSymmetric) = LA.Symmetric(Base. *(A.data,  B.data))
-Base. *(A::SkewSymmetric, B::AbstractMatrix) = Base. *(A.data,  B)
-Base. *(A::AbstractMatrix, B::SkewSymmetric) = Base. *(A,  B.data)
+LA.kron(A::SkewHermitian,B::StridedMatrix) = kron(A.data,B)
+LA.kron(A::StridedMatrix,B::SkewHermitian) = kron(A,B.data)
 
-function LA.dot(x::AbstractVector, A::SkewSymmetric, y::AbstractVector)
-    LA.require_one_based_indexing(x, y)
-    (length(x) == length(y) == size(A, 1)) || throw(DimensionMismatch())
-    data = A.data
-    r = *( zero(eltype(x)),  zero(eltype(A)) , zero(eltype(y)))
-    @inbounds for j = 1:length(y)
-        @simd for i = 1:j-1
-            Aij = data[i,j]
-            r += dot(x[i], Aij, y[j]) + dot(x[j], -Aij, y[i])
-        end
-    end
-    
-    return r
+# to do: these functions should be specialized for SkewHermitian using eigen/eigvals
+for f in (:schur, :schur!)
+    @eval LA.$f(A::SkewHermitian) = LA.$f(A.data)
 end
+
 include("hessenberg.jl")
 include("eigen.jl")
 include("exp.jl")
-# Scaling with Number
-Base. *(A::SkewSymmetric, x::Number) = SkewSymmetric(A.data*x)
-Base. *(x::Number, A::SkewSymmetric) = SkewSymmetric(x*A.data)
-Base. /(A::SkewSymmetric, x::Number) = SkewSymmetric(A.data/x)
-Base. \(A::SkewSymmetric,b::AbstractVecOrMat) = \(A.data,b)
-
-LA.det(A::SkewSymmetric) = det(A.data)
-LA.logdet(A::SkewSymmetric) = logdet(A.data)
-LA.inv(A::SkewSymmetric)  = inv(A.data)
-LA.inv!(A::SkewSymmetric)  = inv!(A.data)
-
-LA.lu(A::SkewSymmetric)  = lu(A.data)
-LA.lu!(A::SkewSymmetric) = lu!(A.data)
-LA.lq(A::SkewSymmetric)  = lq(A.data)
-LA.lq!(A::SkewSymmetric) = lq!(A.data)
-LA.qr(A::SkewSymmetric)  = qr(A.data)
-LA.qr!(A::SkewSymmetric) = qr!(A.data)
-LA.schur(A::SkewSymmetric)=schur(A.data)
-LA.schur!(A::SkewSymmetric)=schur!(A.data)
-LA.rank(A::SkewSymmetric; atol::Real=0, rtol::Real=atol>0 ? 0 : n*ϵ)=rank(A.data;atol,rtol)
-LA.rank(A::SkewSymmetric, rtol::Real)=rank(A.data,rtol)
-
-LA.rdiv!(A::SkewSymmetric,b::Number) = rdiv!(A.data,b)
-LA.ldiv!(A::SkewSymmetric,b::Number) = ldiv!(A.data,b)
-LA.rmul!(A::SkewSymmetric,b::Number) = rmul!(A.data,b)
-LA.lmul!(A::SkewSymmetric,b::Number) = lmul!(A.data,b)
-
-
-
-LA.kron(A::SkewSymmetric,B::AbstractMatrix)=kron(A.data,B)
-LA.kron(A::AbstractMatrix,B::SkewSymmetric)=kron(A,B.data)
 
 end
 
diff --git a/src/eigen.jl b/src/eigen.jl
index 0abf83b..70afc89 100644
--- a/src/eigen.jl
+++ b/src/eigen.jl
@@ -1,105 +1,55 @@
-# This file is a part of Julia. License is MIT: https://julialang.org/license
+# Based on eigen.jl in Julia. License is MIT: https://julialang.org/license
 
 
-"""
-    eigvals!(A,sortby)
-Returns the eigenvalues of A where A is SkewSymmetric, 
-imaginary part of eigenvalues sorted by sortby.
-"""
-@views function LA.eigvals!(A::SkewSymmetric{<:LA.BlasReal,<:StridedMatrix}, sortby::Union{Function,Nothing}=nothing)
+@views function LA.eigvals!(A::SkewHermitian, sortby::Union{Function,Nothing}=nothing)
     vals = skeweigvals!(A)
     !isnothing(sortby) && sort!(vals, by=sortby)
-    return vals.*1im
+    return complex.(0, vals)
 end
-"""
-    eigvals!(A,irange)
-Returns the eigenvalues of A where A is SkewSymmetric, 
-irange specifies the indices of the eigenvalues to search for.
-"""
 
-@views function LA.eigvals!(A::SkewSymmetric{<:LA.BlasReal,<:StridedMatrix,}, irange::UnitRange)
-    vals= skeweigvals!(A,irange)
-    return vals.*1im
+@views function LA.eigvals!(A::SkewHermitian, irange::UnitRange)
+    vals = skeweigvals!(A,irange)
+    return complex.(0, vals)
 end
 
-"""
-    eigvals!(A,vl,vh)
-Returns the eigenvalues of A where A is SkewSymmetric, 
-[vl,vh] defines the range the imaginary part of the eigenvalues of A must be contained in.
-"""
-@views function LA.eigvals!(A::SkewSymmetric{<:LA.BlasReal,<:StridedMatrix}, vl::Real,vh::Real)
-    vals=skeweigvals!(A,-vh,-vl)
-    return vals.*1im
+@views function LA.eigvals!(A::SkewHermitian, vl::Real,vh::Real)
+    vals = skeweigvals!(A,-vh,-vl)
+    return complex.(0, vals)
 end
-"""
-    eigvals(A,sortby)
-Returns the eigenvalues of A where A is SkewSymmetric, 
-imaginary part of eigenvalues sorted by sortby.
-"""
 
-function LA.eigvals(A::SkewSymmetric; sortby::Union{Function,Nothing}=nothing)
-    T = eltype(A)
-    S = LA.eigtype(T)
-    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), sortby)
-end
-"""
-    eigvals(A,irange)
-Returns the eigenvalues of A where A is SkewSymmetric, 
-irange specifies the indices of the eigenvalues to search for.
-"""
-
-function LA.eigvals(A::SkewSymmetric, irange::UnitRange)
-    T = eltype(A)
-    S = LA.eigtype(T)
-    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), irange)
-end
-
-"""
-    eigvals(A,vl,vh)
-Returns the eigenvalues of A where A is SkewSymmetric, 
-[vl,vh] defines the range the imaginary part of the eigenvalues of A must be contained in.
-"""
+# no need to define LA.eigen(...) since the generic methods should work
 
-function LA.eigvals(A::SkewSymmetric, vl::Real, vh::Real)
-    T = eltype(A)
-    S = LA.eigtype(T)
-    eigvals!(S != T ? convert(AbstractMatrix{S}, A) : copy(A), vl, vh)
-end
-
-eigmax(A::SkewSymmetric{<:Real,<:StridedMatrix}) = eigvals(A, size(A, 1):size(A, 1))[1]
-eigmin(A::SkewSymmetric{<:Real,<:StridedMatrix}) = eigvals(A, 1:1)[1]
-
-@views function skeweigvals!(S::SkewSymmetric)
+@views function skeweigvals!(S::SkewHermitian)
     n = size(S.data,1)
     E = sktrd!(S)[2]
     H = SymTridiagonal(zeros(n),E)
     vals = eigvals!(H)
-    return -vals
-    
+    return vals .= .-vals
+
 end
 
-@views function skeweigvals!(S::SkewSymmetric,irange::UnitRange)
+@views function skeweigvals!(S::SkewHermitian,irange::UnitRange)
     n = size(S.data,1)
     E = sktrd!(S)[2]
     H = SymTridiagonal(zeros(n),E)
     vals = eigvals!(H,irange)
-    return -vals
-    
+    return vals .= .-vals
+
 end
 
-@views function skeweigvals!(S::SkewSymmetric,vl::Real,vh::Real)
+@views function skeweigvals!(S::SkewHermitian,vl::Real,vh::Real)
     n = size(S.data,1)
     E = sktrd!(S)[2]
     H = SymTridiagonal(zeros(n),E)
     vals = eigvals!(H,vl,vh)
-    return -vals
+    return vals .= .-vals
 end
 @views function WYform(A,tau,Q)
     n = size(A,1)
     W = zeros(n-1,n-2)
     Yt = zeros(n-2,n-1)
     temp = zeros(n-1)
-    
+
     for i = 1:n-2
         t = tau[i]
         v = A[i+1:n,i]
@@ -124,7 +74,6 @@ end
     for i=1:k
         prevlastv=max(i,prevlastv)
         if tau[i]==0
-            println("coucou\n")
             for j=1:i
                 T[j,i]=0
             end
@@ -161,6 +110,7 @@ function set_nb2(n::Integer)
     end
     return 1
 end
+
 @views function dormqr(n,k,Q,V,tau)
     nb=set_nb2(n)
     T=similar(Q,nb,nb)
@@ -176,18 +126,14 @@ end
     end
 end
 
-
-
-
-
-@views function skeweigen!(S::SkewSymmetric)
+@views function skeweigen!(S::SkewHermitian)
     n = size(S.data,1)
     tau,E = sktrd!(S)
     A = S.data
     s = similar(A,n)
     H = SymTridiagonal(zeros(n),E)
     trisol = eigen!(H)
-    
+
     vals  = trisol.values*1im
     vals .*= -1
     Qdiag = trisol.vectors
@@ -195,10 +141,10 @@ end
     Qr   = similar(A,n,(n+1)÷2)
     Qim  = similar(A,n,n÷2)
     temp = similar(A,n,n)
-    
+
     Q  = diagm(ones(n))
     LA.LAPACK.ormqr!('L','N',A[2:n,1:n-2],tau,Q[2:end,2:end])
-    
+
     Q1 = similar(A,(n+1)÷2,n)
     Q2 = similar(A,n÷2,n)
     @inbounds for j=1:n
@@ -206,7 +152,7 @@ end
             k=(i+1)÷2
             Q1[k,j] = Qdiag[i,j]
             Q2[k,j] = Qdiag[i+1,j]
-        end   
+        end
     end
 
     c = 1
@@ -218,7 +164,7 @@ end
         end
         c *= (-1)
     end
-    
+
     if n%2==1
         k=(n+1)÷2
         @simd for j=1:n
@@ -228,58 +174,33 @@ end
     end
     mul!(temp,Qr,Q1) #temp is Qr
     mul!(Qdiag,Qim,Q2) #Qdiag is Qim
-    
-    return vals,temp,Qdiag
 
-    
+    return vals,temp,Qdiag
 end
-"""
-    eigen!(A)
 
-Returns [val,Re(Q),Im(Q)], containing the eigenvalues in vals, 
-the real part of the eigenvectors in Re(Q) and the Imaginary part of the eigenvectors in Im(Q)
-"""
-LA.eigen!(A::SkewSymmetric{<:LA.BlasReal,<:StridedMatrix}) = skeweigen!(A)
+# FIXME: return an Eigen structure, not a tuple
+LA.eigen!(A::SkewHermitian) = skeweigen!(A)
 
-"""
-    eigen(A)
+LA.eigen(A::SkewHermitian) =
+    LA.eigen!(copyto!(similar(A, LA.eigtype(eltype(A))), A))
 
-Returns [val,Re(Q),Im(Q)], containing the eigenvalues in vals, 
-the real part of the eigenvectors in Re(Q) and the Imaginary part of the eigenvectors in Im(Q)
-"""
-function LA.eigen(A::SkewSymmetric)
-    T = eltype(A)
-    S = LA.eigtype(T)
-    return eigen!(S != T ? convert(AbstractMatrix{S}, A) : copy(A))
-end
-
-@views function LA.svdvals!(A::SkewSymmetric)
+@views function LA.svdvals!(A::SkewHermitian)
     n=size(A,1)
     vals = skeweigvals!(A)
-    @simd for i=1:n
-        vals[i]=abs(vals[i])
-    end
-    sort!(vals;rev=true)
-    return vals
+    vals .= abs.(vals)
+    return sort!(vals; rev=true)
 end
 
-function LA.svdvals(A::SkewSymmetric)
-    return svdvals!(copy(A))
-end
-@views function LA.svd!(A::SkewSymmetric)
+@views function LA.svd!(A::SkewHermitian)
     n=size(A,1)
     eig,Qr,Qim=eigen!(A)
-    vals=similar(A,n)
     U=Qim*1im
     U+=Qr
-    @simd for i=1:n
-        vals[i]=imag(eig[i])
-    end
+    vals = imag.(eig)
     I=sortperm(vals;by=abs,rev=true)
     permute!(vals,I)
     Base.permutecols!!(U,I)
-    V=copy(U)*1im
-    V.*= -1
+    V = U .* -1im
     @inbounds for i=1:n
         if vals[i] < 0
             vals[i]=-vals[i]
@@ -290,6 +211,5 @@ end
     end
     return LA.SVD(U,vals,adjoint(V))
 end
-function LA.svd(A::SkewSymmetric)
-    return svd!(copy(A))
-end
+
+LA.svd(A::SkewHermitian) = svd!(copy(A))
\ No newline at end of file
diff --git a/src/exp.jl b/src/exp.jl
index f6bf5d5..7c2ce60 100644
--- a/src/exp.jl
+++ b/src/exp.jl
@@ -1,10 +1,8 @@
 # This file is a part of Julia. License is MIT: https://julialang.org/license
 
-@views function skewexpm(A::SkewSymmetric)
+function skewexp!(A::SkewHermitian)
     n = size(A,1)
-    if n == 1
-        return exp(A.data)
-    end
+    n == 1 && return fill(exp(A.data[1,1]), 1,1)
     vals,Qr,Qim = skeweigen!(A)
 
     temp2 = similar(A,n,n)
@@ -14,12 +12,11 @@
     Sin=similar(A,n)
 
     @simd for i=1:n
-        @inbounds Cos[i]=cos(imag(vals[i]))
-        @inbounds Sin[i]=sin(imag(vals[i]))
+        @inbounds Sin[i],Cos[i]=sincos(imag(vals[i]))
     end
     C=Diagonal(Cos)
     S=Diagonal(Sin)
-    
+
     mul!(Q1,Qr,C)
     mul!(Q2,Qim,S)
     Q1 .-= Q2
@@ -32,36 +29,27 @@
     return temp2
 end
 
-
-"""
-    exp(A)
-Returns the matrix exponential of A skew-symmetric using the eigenvalue decomposition.
-"""
-@views function Base.exp(A::SkewSymmetric)
-    return skewexpm(copy(A))
+function Base.exp(A::SkewHermitian)
+    return skewexp!(copy(A))
 end
 
-@views function skewcis(A::SkewSymmetric)
+function skewcis!(A::SkewHermitian)
     n = size(A,1)
-    if n == 1
-        return exp(A.data*1im)
-    end
+    n == 1 && return fill(cis(A.data[1,1]), 1,1)
+
     vals,Qr,Qim = skeweigen!(A)
     Q = Qim.*1im
     Q.+=Qr
     temp = similar(Q,n,n)
     temp2 = similar(Q,n,n)
-    eig=similar(A,n)
-    @simd for i=1:n
-        @inbounds eig[i]=exp(-imag(vals[i]))
-    end
+    eig = @. exp(-imag(vals))
     E=Diagonal(eig)
     mul!(temp,Q,E)
     mul!(temp2,temp,adjoint(Q))
     return temp2
 end
 
-@views function skewcos(A::SkewSymmetric)
+@views function skewcos!(A::SkewHermitian)
     n = size(A,1)
     if n == 1
         return exp(A.data*1im)
@@ -71,13 +59,10 @@ end
     temp2 = similar(A,n,n)
     Q1=similar(A,n,n)
     Q2=similar(A,n,n)
-    eig=similar(A,n)
 
-    @simd for i=1:n
-        @inbounds eig[i]=exp(-imag(vals[i]))
-    end
+    eig = @. exp(-imag(vals))
     E=Diagonal(eig)
-    
+
     mul!(Q1,Qr,E)
     mul!(Q2,Qim,E)
     mul!(temp2,Q1,transpose(Qr))
@@ -85,7 +70,7 @@ end
     Q1 .+= temp2
     return Q1
 end
-@views function skewsin(A::SkewSymmetric)
+@views function skewsin!(A::SkewHermitian)
     n = size(A,1)
     if n == 1
         return exp(A.data*1im)
@@ -95,13 +80,10 @@ end
     temp2 = similar(A,n,n)
     Q1=similar(A,n,n)
     Q2=similar(A,n,n)
-    eig=similar(A,n)
 
-    @simd for i=1:n
-        @inbounds eig[i]=exp(-imag(vals[i]))
-    end
+    eig = @. exp(-imag(vals))
     E=Diagonal(eig)
-    
+
     mul!(Q1,Qr,E)
     mul!(Q2,Qim,E)
     mul!(temp2,Q1,transpose(Qim))
@@ -109,37 +91,33 @@ end
     Q1 .-= temp2
     return Q1
 end
-@views function Base.cis(A::SkewSymmetric)
-    return skewcis(copy(A))
-end
-@views function Base.cos(A::SkewSymmetric)
-    return skewcos(copy(A))
-end
-@views function Base.sin(A::SkewSymmetric)
-    return skewsin(copy(A))
-end
+Base.cis(A::SkewHermitian) = skewcis!(copy(A))
+Base.cos(A::SkewHermitian) = skewcos!(copy(A))
+Base.sin(A::SkewHermitian) = skewsin!(copy(A))
 
-function Base.tan(A::SkewSymmetric)
+function Base.tan(A::SkewHermitian)
     E=cis(A)
     S=imag(E)
     C=real(E)
     return C\S
 end
-function Base.sinh(A::SkewSymmetric)
-    S = exp(A)
-    S .-= transpose(S)
-    S .*= 0.5
-    return S
+Base.sinh(A::SkewHermitian) = skewhermitian!(exp(A))
+Base.cosh(A::SkewHermitian) = hermitian!(exp(A))
 
-end
-function Base.cosh(A::SkewSymmetric)
-    C = exp(A)
-    C .+= transpose(C)
-    C .*= 0.5
-    return C
-end
-function Base.tanh(A::SkewSymmetric)
-    E=exp(*(A,2))
+function Base.tanh(A::SkewHermitian)
+    E = skewexp!(2A)
     return (E+LA.I)\(E-LA.I)
-    
+end
+
+# someday this should be in LinearAlgebra: https://github.com/JuliaLang/julia/pull/31836
+function hermitian!(A::AbstractMatrix{<:Number})
+    LA.require_one_based_indexing(A)
+    n = LA.checksquare(A)
+    @inbounds for i in 1:n
+        A[i,i] = real(A[i,i])
+        for j = 1:i-1
+            A[i,j] = A[j,i] = (A[i,j] + A[j,i]')/2
+        end
+    end
+    return LA.Hermitian(A)
 end
diff --git a/src/hessenberg.jl b/src/hessenberg.jl
index 0a428e0..e2f9025 100644
--- a/src/hessenberg.jl
+++ b/src/hessenberg.jl
@@ -1,54 +1,40 @@
-# This file is a part of Julia. License is MIT: https://julialang.org/license
+# Based on hessenberg.jl in Julia. License is MIT: https://julialang.org/license
 
 """
-    SkewHessenberg 
-Type returned by hessenberg(A) when A is skew-symmetric. 
+    SkewHessenberg
+
+Type returned by hessenberg(A) when A is skew-symmetric.
 Contains the tridiagonal reduction in H, the reflectors in V that is UnitLowerTriangular.
 τ contains the scalars relative to the householder reflector i.e Q_i= I - τ_i v_i v_i^T.
 """
-
-struct SkewHessenberg{SH<:Union{Tridiagonal,UpperHessenberg},S<:UnitLowerTriangular,tau<:AbstractVector}
+struct SkewHessenberg{SH<:Tridiagonal,S<:UnitLowerTriangular,tau<:AbstractVector}
     H::SH # Tridiagonal
     V::S # reflector data in unit lower triangular
     τ::tau # more Q (reflector) data
 end
 
 SkewHessenberg(F::SkewHessenberg) = F
+Base.copy(F::SkewHessenberg) = SkewHessenberg(copy(F.H), copy(F.V), copy(F.τ))
+Base.size(F::SkewHessenberg, d) = size(F.H, d)
+Base.size(F::SkewHessenberg) = size(F.H)
 
-copy(F::SkewHessenberg{<:Any,<:UpperHessenberg}) = SkewHessenberg(copy(F.H),copy(F.V), copy(F.τ))
-copy(F::SkewHessenberg{<:Any,<:Tridiagonal}) = SkewHessenberg(copy(F.H),copy(F.V), copy(F.τ))
-size(F::SkewHessenberg, d) = size(F.H, d)
-size(F::SkewHessenberg) = size(F.H)
-
-"""
-    hessenberg!(A)
-
-Returns the tridiagonal reduction of A skew-symmetric.
-The result is returned as a SkewHessenberg structure.
-"""
-
-@views function LA.hessenberg!(A::SkewSymmetric{<:LA.BlasFloat})
+@views function LA.hessenberg!(A::SkewHermitian)
     tau,E = sktrd!(A)
     n = size(A,1)
     return SkewHessenberg(Tridiagonal(E,zeros(n),-E),UnitLowerTriangular(A.data[2:end,1:end-1]),tau)
 end
 
-"""
-    hessenberg(A)
-
-Returns the tridiagonal reduction of A skew-symmetric.
-The result is returned as a SkewHessenberg structure.
-"""
-
-function LA.hessenberg(A::SkewSymmetric{<:LA.BlasFloat})
-    return hessenberg!(copy(A))
+function Base.show(io::IO, mime::MIME"text/plain", F::SkewHessenberg)
+    summary(io, F)
+    println(io, "\nV factor:")
+    show(io, mime, F.V)
+    println(io, "\nH factor:")
+    show(io, mime, F.H)
 end
 
-function Base.display(F::SkewHessenberg)
-    display(F.H)
-    display(F.V)
-    display(F.τ)
-end
+# fixme: to get the Q factor, F.Q should return a type that
+# implicitly multiplies by reflectors, and you should convert it
+# to a matrix with Matrix(F.Q), eliminating getQ.
 
 """
     getQ(H)
@@ -80,7 +66,7 @@ end
     tau = 2/((norm(v)^2))
     return v,tau
 end
-@inline @views function ger2!(tau::Number , v::StridedVector{T} , s::StridedVector{T}, 
+@inline @views function ger2!(tau::Number , v::StridedVector{T} , s::StridedVector{T},
     A::StridedMatrix{T}) where {T<:LA.BlasFloat}
     tau2 = promote(tau, zero(T))[1]
     if tau2 isa Union{Bool,T}
@@ -151,12 +137,12 @@ end
         end
     end
     """
-    
+
     for j=1:n
         @inbounds axpy!(x[j],A[j+1:n,j],y[j+1:n])
         @inbounds y[j] -= dot(A[j+1:n,j],x[j+1:n])
     end
-    
+
     """
     nb=10
     @inbounds for j=1:n
@@ -177,7 +163,7 @@ end
                 @inbounds y[i] += temp1*A[i,j]
                 @inbounds temp2 += A[i,j]*x[i]
             end
-        end 
+        end
         y[j] -= temp2
     end
     """
@@ -202,14 +188,14 @@ end
             mul!(A[i:n,i],A[i:n,1:i-1],W[i,1:i-1],1,1)
             mul!(A[i:n,i],W[i:n,1:i-1],A[i,1:i-1],-1,1)
         end
-        
+
         #Generate elementary reflector H(i) to annihilate A(i+2:n,i)
-        
+
         v,stau = householder_reflector!(A[i+1:n,i],V[i:n-1],n-i)
         A[i+1,i] -= stau*dot(v,A[i+1:n,i])
         E[i]   = A[i+1,i]
         A[i+1:end,i] = v
-        
+
         mul!(W[i+1:n,i],A[i+1:n,i+1:n], A[i+1:n,i])  #Key point 60% of running time of sktrd!
         #skmv!(A[i+1:n,i+1:n], A[i+1:n,i],W[i+1:n,i],n-i)
         if i>1
@@ -219,14 +205,14 @@ end
             mul!(W[i+1:n,i],W[i+1:n,1:i-1],W[1:i-1,i],-1,1)
         end
         W[i+1:n,i] .*= stau
-        
+
         alpha = -stau*dot(W[i+1:n,i],A[i+1:n,i])/2
         W[i+1:n,i].+=alpha.*A[i+1:n,i]
         tau[i] = stau
-        
-        
+
+
     end)
-    return 
+    return
 end
 function set_nb(n::Integer)
     if n<=12
@@ -239,7 +225,7 @@ function set_nb(n::Integer)
     return 1
 end
 
-@views function sktrd!(S::SkewSymmetric)
+@views function sktrd!(S::SkewHermitian{<:Real})
     #println("begin\n")
     n = size(S.data,1)
 
@@ -257,7 +243,7 @@ end
     V   = similar(A, n-1)
 
     oldi = 0
-    
+
     @inbounds(for i = 1:nb:n-nb-2
         size = n-i+1
 
@@ -274,16 +260,16 @@ end
             @inbounds A[s+j,s+j] = 0
         end
         """
-        
+
         for k = 1:n-s
             A[s+k,s+k] = 0
             @simd for j = k+1:n-s
                 A[s+j,s+k] += update[j,k]-update[k,j]
                 A[s+k,s+j] = - A[s+j,s+k]
             end
-            
+
         end
-        
+
         """
         N=n-nb-i+1
         @inbounds (for j=1:N
@@ -297,9 +283,9 @@ end
                     A[s+t,k1] += A[s+t,k2]*temp1-W[nb+t,l]*temp2
                 end
             end
-            
+
             @simd for t=j+1:N
-                A[k1,s+t]=-A[s+t,k1] 
+                A[k1,s+t]=-A[s+t,k1]
             end
         end)
         """
@@ -316,6 +302,6 @@ end
     E[end] = A[end,end-1]
 
     return tau, E
-    
+
 end
 
diff --git a/test/runtests.jl b/test/runtests.jl
index d01933c..291e9d0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,22 +1,19 @@
-using LinearAlgebra
+using LinearAlgebra, Random
 import SkewLinearAlgebra as SLA
 using Test
 
+Random.seed!(314159) # use same pseudorandom stream for every test
+
 @testset "SkewLinearAlgebra.jl" begin
     for n in [2,20,153,200]
-        A=randn(n,n)
-        for i=1:n
-            A[i,i]=0
-            for j=1:i-1
-                A[j,i]=-A[i,j]
-            end
-        end
-
-        A=SLA.SkewSymmetric(A)
-        @test SLA.isskewsymmetric(A)==true
+        @show n
+        A=SLA.skewhermitian(randn(n,n))
+        @test SLA.isskewhermitian(A)
+        @test SLA.isskewhermitian(A.data)
         B=2*Matrix(A)
+        @test SLA.isskewhermitian(B)
 
-        @test A==copy(A)
+        @test A==copy(A)::SLA.SkewHermitian
         @test size(A)==size(A.data)
         @test size(A,1)==size(A.data,1)
         @test size(A,2)==size(A.data,2)
@@ -26,12 +23,15 @@ using Test
         @test A*A == Symmetric(A.data*A.data)
         @test A*B == A.data*B
         @test B*A == B*A.data
+        if iseven(n) # for odd n, a skew-Hermitian matrix is singular
+            @test inv(A)::SLA.SkewHermitian ≈ inv(A.data)
+        end
         @test (A*2).data ==A.data*2
         @test (2*A).data ==2*A.data
         @test (A/2).data == A.data/2
         C=A+A
         @test C.data==A.data+A.data
-        B=SLA.SkewSymmetric(B)
+        B=SLA.SkewHermitian(B)
         C=A-B
         @test C.data==-A.data
         B=triu(A)
@@ -56,7 +56,7 @@ using Test
         @test getindex(A,n-1,n)==-3
         @test parent(A)== A.data
 
-        
+
         x=randn(n)
         y=zeros(n)
         mul!(y,A,x,2,0)
@@ -76,14 +76,14 @@ using Test
         @test C==2*A.data*B
         mul!(C,B,A,2,0)
         @test C==2*B*A.data
-        B=SLA.SkewSymmetric(B)
+        B=SLA.SkewHermitian(B)
         mul!(C,B,A,2,0)
         @test C==2*B.data*A.data
         A.data[n,n]=4
-        @test SLA.isskewsymmetric(A)==false
+        @test SLA.isskewhermitian(A.data)==false
         A.data[n,n]=0
         A.data[n,1]=4
-        @test SLA.isskewsymmetric(A)==false
+        @test SLA.isskewhermitian(A.data)==false
         #LU=lu(A)
         #@test LU.L*LU.U≈A.data
         LQ=lq(A)
@@ -96,14 +96,7 @@ using Test
 end
 @testset "hessenberg.jl" begin
     for n in [2,20,153,200]
-        A=randn(n,n)
-        for i=1:n
-            A[i,i]=0
-            for j=1:i-1
-                A[j,i]=-A[i,j]
-            end
-        end
-        A=SLA.SkewSymmetric(A)
+        A=SLA.skewhermitian(randn(n,n))
         B=Matrix(A)
         HA=hessenberg(A)
         HB=hessenberg(B)
@@ -117,7 +110,7 @@ end
     A=zeros(4,4)
     A[2:4,1]=ones(3)
     A[1,2:4]=-ones(3)
-    A=SLA.SkewSymmetric(A)
+    A=SLA.SkewHermitian(A)
     B=Matrix(A)
     HA=hessenberg(A)
     HB=hessenberg(B)
@@ -126,16 +119,9 @@ end
 end
 @testset "eigen.jl" begin
     for n in [2,20,153,200]
-        A=randn(n,n)
-        for i=1:n
-            A[i,i]=0
-            for j=1:i-1
-                A[j,i]=-A[i,j]
-            end
-        end
-        A=SLA.SkewSymmetric(A)
+        A=SLA.skewhermitian(randn(n,n))
         B=Matrix(A)
-        
+
         valA = imag(eigvals(A))
         valB = imag(eigvals(B))
         sort!(valA)
@@ -158,14 +144,7 @@ end
 @testset "exp.jl" begin
 
     for n in [2,20,153,200]
-        A=randn(n,n)
-        for i=1:n
-            A[i,i]=0
-            for j=1:i-1
-                A[j,i]=-A[i,j]
-            end
-        end
-        A=SLA.SkewSymmetric(A)
+        A=SLA.skewhermitian(randn(n,n))
         B=Matrix(A)
         @test exp(B)≈exp(A)
         @test cis(A)≈exp(Hermitian(A.data*1im))