I am trying to compute the eigenvalues of a symmetric matrix in packed format, i.e. you stack the upper triangular entries in a vector (Packed Storage )
It seems like there exists a LAPACK function for that, e.g. spev
.
However, it seems like that function is not yet supported in Julia’s LAPACK package.
Another way would be to recreate the original matrix and then call eigen!
or LAPACK.syevr!
, like so:
using LinearAlgebra
# represent the matrix [1 2 4; 2 3 5; 4 5 6] in packed form
x_packed = [1;2;3;4;5;6]
function populate_upper_triangle!(A::AbstractMatrix, x::AbstractVector)
k = 0
for j in 1:size(A,1)
for i in 1:j
k += 1
A[i,j] = x[k]
end
end
nothing
end
function packed_eigen(x)
# compute original matrix size
n = Int(sqrt(0.25 + 2 * length(x)) - 0.5)
# allocate memory for orignal matrix
X = zeros(n, n)
populate_upper_triangle!(X, x)
# solve eigenvalue problem
(W, Z) = LAPACK.syevr!('V', 'A', 'U', X, 0.0, 0.0, 0, 0, -1.0);
return W
end
packed_eigen(x_packed)
My question is, what is the most efficient way to do this without allocating a lot of memory (I want to compute the eigenvalues hundreds of times)? I want to avoid
X = zeros(n, n)
populate_upper_triangle!(X, x)
and was envisioning something like reshape
to create a view-based upper triangular matrix that can be passed to LAPACK.syevr!
.