Memory usage of complex hessenberg fixed (#34)

* tridiag.jl

* tridiag.jl

* tridiag.jl

* fixes tridiag.jl

* fixes tridiag.jl

* more fixes

* more fixes

* tridiag fixed

* fixes tridiag

* fixes tridiag

* fixes tridiag

* fixes tridiag

* fixes tridiag

* complex hessenberg

* complex hessenberg 2nd

* complex hessenberg 2nd

* complex hessenberg 2nd

* Memory complex hessenberg

* Memory complex hessenberg

src/complexhessenberg.jl

@@ -7,48 +7,51 @@
 end
 LA.hessenberg(A::SkewHermitian{<:Complex})=hessenberg!(copy(A))
 
-@views function complex_householder_reflector!(x,v,n)
+
+@views function complex_householder_reflector!(x,n)
+
     if n>1
-        xnorm=norm(x[2:end])
-        alpha=x[1]
-        alphar=real(alpha)
-        alphaim=imag(alpha)
-        if xnorm > 1e-15
-            if alphar>0
-                beta=-sqrt(alphar*alphar+alphaim*alphaim+xnorm*xnorm)
-            else
-                beta=sqrt(alphar*alphar+alphaim*alphaim+xnorm*xnorm)
-            end
-            tau=(beta-alphar)/beta-alphaim/beta*1im
-            alpha=1/(alpha-beta)
-            v[1]=1
-            @inbounds v[2:end]=x[2:end]
-            v[2:end] .*= alpha
-            alpha=beta
-        else
-            tau = convert(eltype(x),0)
-            v = zeros(eltype(x),n)
-            alpha=0
-        end
+        xnorm = norm(x[2:end])
     else
-        alpha=x[1]
-        alphar=real(alpha)
-        alphaim=imag(alpha)
+        xnorm = zero(real.(x[1]))
+    end
+    T=eltype(x)
+    alpha=x[1]
+    alphar=real.(alpha)
+    alphaim=imag.(alpha)
+    if xnorm > 1e-15 || n==1
+
         if alphar>0
-            beta=-sqrt(alphar*alphar+alphaim*alphaim)
+            beta=-sqrt(alphar*alphar+alphaim*alphaim+xnorm*xnorm)
         else
-            beta=sqrt(alphar*alphar+alphaim*alphaim)
+            beta=sqrt(alphar*alphar+alphaim*alphaim+xnorm*xnorm)
+        end
+        tau = complex.((beta-alphar)/beta,-alphaim/beta)
+        beta= convert(T,beta)
+        alpha = 1/(alpha-beta)
+        x[1] = convert(T,1)
+        alpha=convert(T,alpha)
+
+        if n>1
+            @inbounds x[2:n].*=alpha
         end
-        tau=(beta-alphar)/beta-alphaim/beta*1im
-        alpha=1/(alpha-beta)
-        v[1]=1
+
+
         alpha=beta
+
+    else
+        tau = convert(eltype(x),0)
+        x = zeros(eltype(x),n)
+        alpha=convert(eltype(x),0)
     end
-    return v,tau, alpha
+
+    return tau, alpha
+
 end
 @views function cger2!(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}
         return LA.BLAS.ger!(tau2, v, s, A)
     else
@@ -73,12 +76,11 @@ end
 @views function complexskewhess!(A::AbstractMatrix,tau::AbstractVector,E::AbstractVector)
     n = size(A,1)
     atmp = similar(A,n)
-    vtmp = similar(atmp)
     @inbounds (for i=1:n-1
-        v,stau,alpha = complex_householder_reflector!(A[i+1:end,i], vtmp[i+1:end],n-i)
-
+        stau,alpha = complex_householder_reflector!(A[i+1:end,i],n-i)
+        @views v=A[i+1:end,i]
         E[i] = alpha
-        A[i+1:end,i]=v
+
         complexleftHouseholder!(A[i+1:end,i+1:end],v,atmp[i+1:end],stau)
 
         s = mul!(atmp[i+1:end], A[i+1:end,i+1:end], v)
@@ -97,22 +99,30 @@ end
 
 
 
-@views function complexlatrd!(A::AbstractMatrix,E::AbstractVector,W::AbstractMatrix,V::AbstractVector,tau::AbstractVector,n::Number,nb::Number)
+
+@views function complexlatrd!(A::AbstractMatrix,E::AbstractVector,W::AbstractMatrix,tau::AbstractVector,tempconj::AbstractVector,n::Number,nb::Number)
 
     @inbounds(for i=1:nb
         if i>1
-
-            mul!(A[i:n,i],A[i:n,1:i-1],conj.(W[i,1:i-1]),1,1)
-            mul!(A[i:n,i],W[i:n,1:i-1],conj.(A[i,1:i-1]),-1,1)
+            @simd for j=1:i-1
+                tempconj[j] = conj.(W[i,j])
+            end
+            mul!(A[i:n,i],A[i:n,1:i-1],tempconj[1:i-1],1,1)
+            @simd for j=1:i-1
+                tempconj[j] = conj.(A[i,j])
+            end
+            mul!(A[i:n,i],W[i:n,1:i-1],tempconj[1:i-1],-1,1)
+
 
 
         end
 
         #Generate elementary reflector H(i) to annihilate A(i+2:n,i)
 
-        v,stau,alpha = complex_householder_reflector!(A[i+1:n,i],V[i:n-1],n-i)
+
+        stau,alpha = complex_householder_reflector!(A[i+1:n,i],n-i)
         E[i]   = real(alpha)
-        A[i+1:n,i] = v
+
 
         mul!(W[i+1:n,i],A[i+1:n,i+1:n], A[i+1:n,i],1,0)  
         if i>1
@@ -138,7 +148,9 @@ function set_nb2(n::Integer)
     elseif n<=100
         return 10
     else
-        return 20
+
+        return 60
+
     end
     return 1
 end
@@ -158,14 +170,18 @@ end
     tau = similar(A,n-1)
     W   = similar(A, n, nb)
     update = similar(A, n-nb, n-nb)
-    V   = similar(A, n-1)
+
+    tempconj=similar(A,nb)
+
 
     oldi = 0
 
     @inbounds(for i = 1:nb:n-nb-1
         size = n-i+1
 
-        complexlatrd!(A[i:n,i:n],E[i:i+nb-1],W,V,tau[i:i+nb-1],size,nb)
+
+        complexlatrd!(A[i:n,i:n],E[i:i+nb-1],W,tau[i:i+nb-1],tempconj,size,nb)
+
         mul!(update[1:n-nb-i+1,1:n-nb-i+1],A[i+nb:n,i:i+nb-1],adjoint(W[nb+1:size,:]))
 
         s = i+nb-1

test/runtests.jl

@@ -259,12 +259,18 @@ end
 end
 
 #=
+
 using BenchmarkTools
 n=1000
 A = SLA.skewhermitian(randn(n,n)+1im*randn(n,n))
 B = Hermitian(A.data*1im)
-#@btime hessenberg(B)
+
+C=Matrix(A)
+@btime hessenberg(B)
 @btime hessenberg(A)
+@btime hessenberg(C)
+a=1
+
 =#