From da2fb1b971ff45be00dfaedf4e429bfdc4f84298 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Sun, 2 Feb 2025 10:26:57 -0500 Subject: [PATCH 01/15] resizings --- src/lapack.jl | 42 ++++++++++++++++++++++-------------------- 1 file changed, 22 insertions(+), 20 deletions(-) diff --git a/src/lapack.jl b/src/lapack.jl index a0c1c785..f698b654 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -2835,7 +2835,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in end end if n < size(A,2) - A[:,1:n] + reshape(resize!(vec(A), m * n), m, n) else A end @@ -2871,7 +2871,7 @@ for (orglq, orgqr, orgql, orgrq, ormlq, ormqr, ormql, ormrq, gemqrt, elty) in end end if n < size(A,2) - A[:,1:n] + reshape(resize!(vec(A), m * n), m, n) else A end @@ -3736,22 +3736,24 @@ for (trcon, trevc, trrfs, elty) in work, info, 1, 1) chklapackerror(info[]) + VLn = size(VL, 1) + VRn = size(VR, 1) #Decide what exactly to return if howmny == 'S' #compute selected eigenvectors if side == 'L' #left eigenvectors only - return select, VL[:,1:m[]] + return select, reshape(resize!(vec(VL), VLn * m[]), VLn, m[]) elseif side == 'R' #right eigenvectors only - return select, VR[:,1:m[]] + return select, reshape(resize!(vec(VR), VRn * m[]), VRn, m[]) else #side == 'B' #both eigenvectors - return select, VL[:,1:m[]], VR[:,1:m[]] + return select, reshape(resize!(vec(VL), VLn * m[]), VLn, m[]), reshape(resize!(vec(VR), VRn * m[]), VRn, m[]) end else #compute all eigenvectors if side == 'L' #left eigenvectors only - return VL[:,1:m[]] + return reshape(resize!(vec(VL), VLn * m[]), VLn, m[]) elseif side == 'R' #right eigenvectors only - return VR[:,1:m[]] + return reshape(resize!(vec(VR), VRn * m[]), VRn, m[]) else #side == 'B' #both eigenvectors - return VL[:,1:m[]], VR[:,1:m[]] + return reshape(resize!(vec(VL), VLn * m[]), VLn, m[]), reshape(resize!(vec(VR), VRn * m[]), VRn, m[]) end end end @@ -3873,22 +3875,21 @@ for (trcon, trevc, trrfs, elty, relty) in work, rwork, info, 1, 1) chklapackerror(info[]) - #Decide what exactly to return if howmny == 'S' #compute selected eigenvectors if side == 'L' #left eigenvectors only - return select, VL[:,1:m[]] + return select, reshape(resize!(vec(VL), VLn * m[]), VLn, m[]) elseif side == 'R' #right eigenvectors only - return select, VR[:,1:m[]] - else #side=='B' #both eigenvectors - return select, VL[:,1:m[]], VR[:,1:m[]] + return select, reshape(resize!(vec(VR), VRn * m[]), VRn, m[]) + else #side == 'B' #both eigenvectors + return select, reshape(resize!(vec(VL), VLn * m[]), VLn, m[]), reshape(resize!(vec(VR), VRn * m[]), VRn, m[]) end else #compute all eigenvectors if side == 'L' #left eigenvectors only - return VL[:,1:m[]] + return reshape(resize!(vec(VL), VLn * m[]), VLn, m[]) elseif side == 'R' #right eigenvectors only - return VR[:,1:m[]] - else #side=='B' #both eigenvectors - return VL[:,1:m[]], VR[:,1:m[]] + return reshape(resize!(vec(VR), VRn * m[]), VRn, m[]) + else #side == 'B' #both eigenvectors + return reshape(resize!(vec(VL), VLn * m[]), VLn, m[]), reshape(resize!(vec(VR), VRn * m[]), VRn, m[]) end end end @@ -4033,7 +4034,7 @@ for (stev, stebz, stegr, stein, elty) in w, iblock, isplit, work, iwork, info, 1, 1) chklapackerror(info[]) - w[1:m[]], iblock[1:m[]], isplit[1:nsplit[1]] + resize!(w, m[]), resize!(iblock, m[]), resize!(isplit, nsplit[1]) end function stegr!(jobz::AbstractChar, range::AbstractChar, dv::AbstractVector{$elty}, ev::AbstractVector{$elty}, vl::Real, vu::Real, il::Integer, iu::Integer) @@ -4056,7 +4057,8 @@ for (stev, stebz, stegr, stein, elty) in m = Ref{BlasInt}() w = similar(dv, $elty, n) ldz = jobz == 'N' ? 1 : n - Z = similar(dv, $elty, ldz, range == 'I' ? iu-il+1 : n) + Zn = range == 'I' ? iu-il+1 : n + Z = similar(dv, $elty, ldz * Zn) isuppz = similar(dv, BlasInt, 2*size(Z, 2)) work = Vector{$elty}(undef, 1) lwork = BlasInt(-1) @@ -4085,7 +4087,7 @@ for (stev, stebz, stegr, stein, elty) in resize!(iwork, liwork) end end - m[] == length(w) ? w : w[1:m[]], m[] == size(Z, 2) ? Z : Z[:,1:m[]] + m[] == length(w) ? w : resize!(w, m[]), m[] == size(Z, 2) ? Z : reshape(resize!(Z, ldz * m[]), ldz, m[]) end function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) From 631a3b8e0fe2b638dd5b9eed4c93c36a5f523eea Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Sun, 2 Feb 2025 14:09:02 -0500 Subject: [PATCH 02/15] fixes --- src/lapack.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/lapack.jl b/src/lapack.jl index f698b654..34530255 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -3875,6 +3875,9 @@ for (trcon, trevc, trrfs, elty, relty) in work, rwork, info, 1, 1) chklapackerror(info[]) + VLn = size(VL, 1) + VRn = size(VR, 1) + #Decide what exactly to return if howmny == 'S' #compute selected eigenvectors if side == 'L' #left eigenvectors only return select, reshape(resize!(vec(VL), VLn * m[]), VLn, m[]) From 78d1fc80c2a220acb44fa22fbc2144343a7d2fda Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Sun, 2 Feb 2025 14:34:35 -0500 Subject: [PATCH 03/15] simplify (~no-op or cheap if same size) --- src/lapack.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lapack.jl b/src/lapack.jl index 34530255..e0300219 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -4090,7 +4090,7 @@ for (stev, stebz, stegr, stein, elty) in resize!(iwork, liwork) end end - m[] == length(w) ? w : resize!(w, m[]), m[] == size(Z, 2) ? Z : reshape(resize!(Z, ldz * m[]), ldz, m[]) + resize!(w, m[]), reshape(resize!(Z, ldz * m[]), ldz, m[]) end function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) From ebaa8b2ad8afd7f1c228d5cc1a4df579f21669ad Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Tue, 11 Feb 2025 10:12:52 -0500 Subject: [PATCH 04/15] Test debugging info for CI --- test/eigen.jl | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/test/eigen.jl b/test/eigen.jl index a82c7454..51c0dda1 100644 --- a/test/eigen.jl +++ b/test/eigen.jl @@ -16,6 +16,7 @@ Random.seed!(12343219) areal = randn(n,n)/2 aimg = randn(n,n)/2 +# DEBUGGING INFO @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int) aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) asym = aa' + aa # symmetric indefinite @@ -24,6 +25,46 @@ aimg = randn(n,n)/2 (view(aa, 1:n, 1:n), view(asym, 1:n, 1:n), view(apd, 1:n, 1:n))) + @testset "non-symmetric eigen decomposition" begin + for T in (Hermitian(Tridiagonal(a), :U), Hermitian(Tridiagonal(a), :L)) + @info typeof(T) + f = eigen(T) + end + end + @testset "symmetric generalized eigenproblem" begin + if isa(a, Array) + asym_sg = asym[1:n1, 1:n1] + a_sg = a[:,n1+1:n2] + else + asym_sg = view(asym, 1:n1, 1:n1) + a_sg = view(a, 1:n, n1+1:n2) + end + ASG2 = a_sg'a_sg + + # matrices of different types (#14896) + D = Diagonal(ASG2) + for uplo in (:L, :U) + gd = eigen(Hermitian(Tridiagonal(ASG2), uplo), D) + @test Hermitian(Tridiagonal(ASG2), uplo) * gd.vectors ≈ D * gd.vectors * Diagonal(gd.values) + end + end + end +end + + + +@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int) + @info "eigen.jl" + @info eltya + aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) + @info repr(aa) + asym = aa' + aa # symmetric indefinite + apd = aa' * aa # symmetric positive-definite + for (a, asym, apd) in ((aa, asym, apd), + (view(aa, 1:n, 1:n), + view(asym, 1:n, 1:n), + view(apd, 1:n, 1:n))) + @info typeof(a) ε = εa = eps(abs(float(one(eltya)))) α = rand(eltya) @@ -33,6 +74,7 @@ aimg = randn(n,n)/2 @test eab.vectors == eigvecs(fill(α,1,1),fill(β,1,1)) @testset "non-symmetric eigen decomposition" begin + @info "non-symmetric eigen decomposition" d, v = eigen(a) for i in 1:size(a,2) @test a*v[:,i] ≈ d[i]*v[:,i] @@ -46,6 +88,7 @@ aimg = randn(n,n)/2 @test Array(f) ≈ a for T in (Tridiagonal(a), Hermitian(Tridiagonal(a), :U), Hermitian(Tridiagonal(a), :L)) + @info T f = eigen(T) d, v = f for i in 1:size(a,2) @@ -65,6 +108,7 @@ aimg = randn(n,n)/2 @test_throws DomainError eigmax(a - a') end @testset "symmetric generalized eigenproblem" begin + @info "symmetric generalized eigenproblem" if isa(a, Array) asym_sg = asym[1:n1, 1:n1] a_sg = a[:,n1+1:n2] @@ -100,6 +144,7 @@ aimg = randn(n,n)/2 # matrices of different types (#14896) D = Diagonal(ASG2) for uplo in (:L, :U) + @info uplo if eltya <: Real fs = eigen(Symmetric(asym_sg, uplo), ASG2) @test fs.values ≈ f.values From d12186a8a6a5ab8bcccf40ba89f88967c124dfd4 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Tue, 11 Feb 2025 12:14:57 -0500 Subject: [PATCH 05/15] remove debug lines --- test/eigen.jl | 45 --------------------------------------------- 1 file changed, 45 deletions(-) diff --git a/test/eigen.jl b/test/eigen.jl index 51c0dda1..a82c7454 100644 --- a/test/eigen.jl +++ b/test/eigen.jl @@ -16,7 +16,6 @@ Random.seed!(12343219) areal = randn(n,n)/2 aimg = randn(n,n)/2 -# DEBUGGING INFO @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int) aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) asym = aa' + aa # symmetric indefinite @@ -25,46 +24,6 @@ aimg = randn(n,n)/2 (view(aa, 1:n, 1:n), view(asym, 1:n, 1:n), view(apd, 1:n, 1:n))) - @testset "non-symmetric eigen decomposition" begin - for T in (Hermitian(Tridiagonal(a), :U), Hermitian(Tridiagonal(a), :L)) - @info typeof(T) - f = eigen(T) - end - end - @testset "symmetric generalized eigenproblem" begin - if isa(a, Array) - asym_sg = asym[1:n1, 1:n1] - a_sg = a[:,n1+1:n2] - else - asym_sg = view(asym, 1:n1, 1:n1) - a_sg = view(a, 1:n, n1+1:n2) - end - ASG2 = a_sg'a_sg - - # matrices of different types (#14896) - D = Diagonal(ASG2) - for uplo in (:L, :U) - gd = eigen(Hermitian(Tridiagonal(ASG2), uplo), D) - @test Hermitian(Tridiagonal(ASG2), uplo) * gd.vectors ≈ D * gd.vectors * Diagonal(gd.values) - end - end - end -end - - - -@testset for eltya in (Float32, Float64, ComplexF32, ComplexF64, Int) - @info "eigen.jl" - @info eltya - aa = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal) - @info repr(aa) - asym = aa' + aa # symmetric indefinite - apd = aa' * aa # symmetric positive-definite - for (a, asym, apd) in ((aa, asym, apd), - (view(aa, 1:n, 1:n), - view(asym, 1:n, 1:n), - view(apd, 1:n, 1:n))) - @info typeof(a) ε = εa = eps(abs(float(one(eltya)))) α = rand(eltya) @@ -74,7 +33,6 @@ end @test eab.vectors == eigvecs(fill(α,1,1),fill(β,1,1)) @testset "non-symmetric eigen decomposition" begin - @info "non-symmetric eigen decomposition" d, v = eigen(a) for i in 1:size(a,2) @test a*v[:,i] ≈ d[i]*v[:,i] @@ -88,7 +46,6 @@ end @test Array(f) ≈ a for T in (Tridiagonal(a), Hermitian(Tridiagonal(a), :U), Hermitian(Tridiagonal(a), :L)) - @info T f = eigen(T) d, v = f for i in 1:size(a,2) @@ -108,7 +65,6 @@ end @test_throws DomainError eigmax(a - a') end @testset "symmetric generalized eigenproblem" begin - @info "symmetric generalized eigenproblem" if isa(a, Array) asym_sg = asym[1:n1, 1:n1] a_sg = a[:,n1+1:n2] @@ -144,7 +100,6 @@ end # matrices of different types (#14896) D = Diagonal(ASG2) for uplo in (:L, :U) - @info uplo if eltya <: Real fs = eigen(Symmetric(asym_sg, uplo), ASG2) @test fs.values ≈ f.values From 06a05f414d7717ad33094a3016a9d739d976b6bc Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Tue, 11 Feb 2025 12:15:04 -0500 Subject: [PATCH 06/15] Try taking min --- src/lapack.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/lapack.jl b/src/lapack.jl index e0300219..8abe0181 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -4090,7 +4090,8 @@ for (stev, stebz, stegr, stein, elty) in resize!(iwork, liwork) end end - resize!(w, m[]), reshape(resize!(Z, ldz * m[]), ldz, m[]) + Zm = min(size(Z, 2), m[]) + resize!(w, m[]), reshape(resize!(Z, ldz * Zm), ldz, Zm) end function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) From e519b44eb0157efca352d6fe1ba63fe6358c73b0 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Tue, 11 Feb 2025 12:17:12 -0500 Subject: [PATCH 07/15] fix check --- src/lapack.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lapack.jl b/src/lapack.jl index 8abe0181..fb92f5df 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -4090,7 +4090,7 @@ for (stev, stebz, stegr, stein, elty) in resize!(iwork, liwork) end end - Zm = min(size(Z, 2), m[]) + Zm = min(Zn, m[]) resize!(w, m[]), reshape(resize!(Z, ldz * Zm), ldz, Zm) end From 4bf986bccf819836c292b2382530515be15c453f Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Tue, 11 Feb 2025 12:27:34 -0500 Subject: [PATCH 08/15] Also add for `w`. --- src/lapack.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/lapack.jl b/src/lapack.jl index fb92f5df..aa56a7fd 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -4090,8 +4090,9 @@ for (stev, stebz, stegr, stein, elty) in resize!(iwork, liwork) end end + wm = min(n, m[]) Zm = min(Zn, m[]) - resize!(w, m[]), reshape(resize!(Z, ldz * Zm), ldz, Zm) + resize!(w, wm), reshape(resize!(Z, ldz * Zm), ldz, Zm) end function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) From 4835a0fe8f57f04ee8dd73ae5b5426597a0cd307 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Tue, 11 Feb 2025 22:03:56 -0500 Subject: [PATCH 09/15] Add debug for CI --- src/lapack.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/src/lapack.jl b/src/lapack.jl index aa56a7fd..b90b2a99 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -4092,7 +4092,16 @@ for (stev, stebz, stegr, stein, elty) in end wm = min(n, m[]) Zm = min(Zn, m[]) - resize!(w, wm), reshape(resize!(Z, ldz * Zm), ldz, Zm) + try + return resize!(w, wm), reshape(resize!(Z, ldz * Zm), ldz, Zm) + catch ex + if isa(ex, ReadOnlyMemoryError) + @info "ReadOnlyMemoryError!" n Zn wm Zm m[] jobz range dv ev vl vu il iu + @info w + @info Z + end + rethrow() + end end function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) From ff355964c55a6987d7a0711f01f183d8d866e510 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 12 Feb 2025 04:40:08 -0500 Subject: [PATCH 10/15] Remove debug again --- src/lapack.jl | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/src/lapack.jl b/src/lapack.jl index b90b2a99..77c9e8a4 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -4092,16 +4092,7 @@ for (stev, stebz, stegr, stein, elty) in end wm = min(n, m[]) Zm = min(Zn, m[]) - try - return resize!(w, wm), reshape(resize!(Z, ldz * Zm), ldz, Zm) - catch ex - if isa(ex, ReadOnlyMemoryError) - @info "ReadOnlyMemoryError!" n Zn wm Zm m[] jobz range dv ev vl vu il iu - @info w - @info Z - end - rethrow() - end + return resize!(w, wm), reshape(resize!(Z, ldz * Zm), ldz, Zm) end function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) From 426bff6b4d57c1cf19f4ce3a48220cfbb7c61de8 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 12 Feb 2025 05:04:21 -0500 Subject: [PATCH 11/15] Add tests to cover cases where output is truncated --- test/lapack.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/test/lapack.jl b/test/lapack.jl index f05d7d99..a59b66f5 100644 --- a/test/lapack.jl +++ b/test/lapack.jl @@ -354,6 +354,10 @@ end @test_throws DimensionMismatch LAPACK.ormqr!('L','N',A,temp,B) @test_throws ArgumentError LAPACK.ormqr!('X','N',A,temp,B) @test_throws ArgumentError LAPACK.ormqr!('L','X',A,temp,B) + A = rand(elty,10,11) + A,tau = LAPACK.geqrf!(A) + B = copy(A) + @test LAPACK.orgqr!(B,tau) ≈ LAPACK.ormqr!('R','N',A,tau,Matrix{elty}(I, 10, 10)) A = rand(elty,10,10) A,tau = LAPACK.geqlf!(A) @@ -372,6 +376,11 @@ end @test_throws ArgumentError LAPACK.ormql!('X','N',A,temp,B) @test_throws ArgumentError LAPACK.ormql!('L','X',A,temp,B) + A = rand(elty,10,11) + A,tau = LAPACK.geqlf!(A) + B = copy(A) + @test LAPACK.orgql!(B,tau) ≈ LAPACK.ormql!('R','N',A,tau,Matrix{elty}(I, 10, 10)) + A = rand(elty,10,10) A,tau = LAPACK.gerqf!(A) @test_throws DimensionMismatch LAPACK.orgrq!(A,tau,11) From f3a923c5614c4eab0cab7dddafa28d06141750b4 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Wed, 12 Feb 2025 06:08:23 -0500 Subject: [PATCH 12/15] try to catch elsewhere? --- test/lapack.jl | 1 + test/tridiag.jl | 10 +++++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/test/lapack.jl b/test/lapack.jl index a59b66f5..da83e4a8 100644 --- a/test/lapack.jl +++ b/test/lapack.jl @@ -354,6 +354,7 @@ end @test_throws DimensionMismatch LAPACK.ormqr!('L','N',A,temp,B) @test_throws ArgumentError LAPACK.ormqr!('X','N',A,temp,B) @test_throws ArgumentError LAPACK.ormqr!('L','X',A,temp,B) + A = rand(elty,10,11) A,tau = LAPACK.geqrf!(A) B = copy(A) diff --git a/test/tridiag.jl b/test/tridiag.jl index 4b592a87..76118558 100644 --- a/test/tridiag.jl +++ b/test/tridiag.jl @@ -383,7 +383,15 @@ end w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a) evecs = LAPACK.stein!(b, a, w) - (e, v) = eigen(SymTridiagonal(b, a)) + try + (e, v) = eigen(SymTridiagonal(b, a)) + catch ex + if isa(ex, ReadOnlyMemoryError) + @info "rom error" evecs b a w + else + rethrow(ex) + end + end @test e ≈ w test_approx_eq_vecs(v, evecs) end From 605c2e6392fe112a08614995b4267615548a6456 Mon Sep 17 00:00:00 2001 From: Nicholas Bauer Date: Mon, 17 Feb 2025 21:06:52 -0500 Subject: [PATCH 13/15] duh --- src/lapack.jl | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/lapack.jl b/src/lapack.jl index 77c9e8a4..001cfc82 100644 --- a/src/lapack.jl +++ b/src/lapack.jl @@ -4062,7 +4062,7 @@ for (stev, stebz, stegr, stein, elty) in ldz = jobz == 'N' ? 1 : n Zn = range == 'I' ? iu-il+1 : n Z = similar(dv, $elty, ldz * Zn) - isuppz = similar(dv, BlasInt, 2*size(Z, 2)) + isuppz = similar(dv, BlasInt, 2 * Zn) work = Vector{$elty}(undef, 1) lwork = BlasInt(-1) iwork = Vector{BlasInt}(undef, 1) @@ -4090,9 +4090,7 @@ for (stev, stebz, stegr, stein, elty) in resize!(iwork, liwork) end end - wm = min(n, m[]) - Zm = min(Zn, m[]) - return resize!(w, wm), reshape(resize!(Z, ldz * Zm), ldz, Zm) + return resize!(w, m[]), reshape(resize!(Z, ldz * m[]), ldz, m[]) end function stein!(dv::AbstractVector{$elty}, ev_in::AbstractVector{$elty}, w_in::AbstractVector{$elty}, iblock_in::AbstractVector{BlasInt}, isplit_in::AbstractVector{BlasInt}) From 5dd14c06eb72e7de3c786e32e379059d33177849 Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 18 Feb 2025 11:05:30 +0100 Subject: [PATCH 14/15] remove try-catch --- test/tridiag.jl | 10 +--------- 1 file changed, 1 insertion(+), 9 deletions(-) diff --git a/test/tridiag.jl b/test/tridiag.jl index 76118558..4b592a87 100644 --- a/test/tridiag.jl +++ b/test/tridiag.jl @@ -383,15 +383,7 @@ end w, iblock, isplit = LAPACK.stebz!('V', 'B', -infinity, infinity, 0, 0, zero, b, a) evecs = LAPACK.stein!(b, a, w) - try - (e, v) = eigen(SymTridiagonal(b, a)) - catch ex - if isa(ex, ReadOnlyMemoryError) - @info "rom error" evecs b a w - else - rethrow(ex) - end - end + (e, v) = eigen(SymTridiagonal(b, a)) @test e ≈ w test_approx_eq_vecs(v, evecs) end From ee7ef2d461330a9819841a8005816414074f97fe Mon Sep 17 00:00:00 2001 From: Daniel Karrasch Date: Tue, 18 Feb 2025 13:24:17 +0100 Subject: [PATCH 15/15] include missing test case (all evecs, both sides) --- test/lapack.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/test/lapack.jl b/test/lapack.jl index da83e4a8..61846399 100644 --- a/test/lapack.jl +++ b/test/lapack.jl @@ -748,6 +748,11 @@ end select,Vln,Vrn = LAPACK.trevc!('B','S',select,copy(T)) @test Vrn ≈ v @test Vln ≈ Vl + Vl = LAPACK.trevc!('L','A',select,copy(T)) + Vr = LAPACK.trevc!('R','A',select,copy(T)) + Vla, Vra = LAPACK.trevc!('B','A',select,copy(T)) + @test Vr ≈ Vra + @test Vl ≈ Vla @test_throws ArgumentError LAPACK.trevc!('V','S',select,T) @test_throws ArgumentError LAPACK.trevc!('R','X',select,T) temp1010 = rand(elty,10,10)