Skip to content

Commit 5b4ee3c

Browse files
authored
Merge pull request #120 from tomhaber/no_auto_shift_invert
remove automatic shift-invert when which==:SM
2 parents 6baea77 + 0b05c36 commit 5b4ee3c

File tree

3 files changed

+50
-21
lines changed

3 files changed

+50
-21
lines changed

src/Arpack.jl

Lines changed: 3 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -117,12 +117,6 @@ function _eigs(A, B;
117117
end
118118
isshift && which == :SM && @warn "Use of :SM in shift-and-invert mode is not recommended, use :LM to find eigenvalues closest to sigma"
119119

120-
if which==:SM && !isshift # transform into shift-and-invert method with sigma = 0
121-
isshift=true
122-
sigma=zero(T)
123-
which=:LM
124-
end
125-
126120
if (explicittransform==:auto)
127121
# Try to automatically detect if it is good to carry out an explicittransform
128122
if (isgeneral && (isshift || which==:LM))
@@ -173,6 +167,9 @@ function _eigs(A, B;
173167
end
174168

175169
whichstr = "LM"
170+
if which == :SM
171+
whichstr = "SM"
172+
end
176173
if which == :BE
177174
whichstr = "BE"
178175
end

src/libarpack.jl

Lines changed: 17 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -201,17 +201,20 @@ function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat,
201201
select = Vector{BlasInt}(undef, ncv)
202202
info = Ref{BlasInt}(0)
203203

204-
dmap = x -> abs.(x)
204+
dmap = if which == "LM" || which == "SM"
205+
abs
206+
elseif which == "LR" || which == "LA" || which == "BE" || which == "SR" || which == "SA"
207+
real
208+
elseif which == "LI" || which == "SI"
209+
abs imag # ARPACK returns largest,smallest abs(imaginary) (complex pairs come together)
210+
else
211+
error("unknown which string $which")
212+
end
213+
214+
rev = which[1] == 'L'
215+
205216
if iparam[7] == 3 # shift-and-invert
206-
dmap = x -> abs.(1 ./ (x .- sigma))
207-
elseif which == "LR" || which == "LA" || which == "BE"
208-
dmap = real
209-
elseif which == "SR" || which == "SA"
210-
dmap = x -> -real(x)
211-
elseif which == "LI"
212-
dmap = imag
213-
elseif which == "SI"
214-
dmap = x -> -imag(x)
217+
dmap = dmap (x -> 1 / (x - sigma))
215218
end
216219

217220
if cmplx
@@ -225,7 +228,7 @@ function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat,
225228
throw(XYEUPD_Exception(info[]))
226229
end
227230

228-
p = sortperm(dmap(d[1:nev]), rev=true)
231+
p = sortperm(d[1:nev], by=dmap, rev=rev)
229232
return ritzvec ? (d[p], v[1:n, p],iparam[5],iparam[3],iparam[9],resid) : (d[p],iparam[5],iparam[3],iparam[9],resid)
230233
elseif sym
231234
d = Vector{T}(undef, nev)
@@ -237,7 +240,7 @@ function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat,
237240
throw(XYEUPD_Exception(info[]))
238241
end
239242

240-
p = sortperm(dmap(d), rev=true)
243+
p = sortperm(d, by=dmap, rev=rev)
241244
return ritzvec ? (d[p], v[1:n, p],iparam[5],iparam[3],iparam[9],resid) : (d[p],iparam[5],iparam[3],iparam[9],resid)
242245
else
243246
dr = Vector{T}(undef, nev+1)
@@ -278,9 +281,9 @@ function eupd_wrapper(T, n::Integer, sym::Bool, cmplx::Bool, bmat,
278281
d = complex.(dr, di)
279282

280283
if j == nev+1
281-
p = sortperm(dmap(d[1:nev]), rev=true)
284+
p = sortperm(d[1:nev], by=dmap, rev=rev)
282285
else
283-
p = sortperm(dmap(d), rev=true)
286+
p = sortperm(d, by=dmap, rev=rev)
284287
p = p[1:nev]
285288
end
286289

test/runtests.jl

Lines changed: 30 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -116,7 +116,7 @@ end
116116
-1.0 0.0 0.0 0.0 0.0 0.0 7.0 1.0
117117
-1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 8.0
118118
]
119-
d, = eigs(A6965,which=:SM,nev=2,ncv=4,tol=eps())
119+
d, = eigs(A6965,which=:LM,nev=2,ncv=4,tol=eps(), sigma=0.0)
120120
@test d[1] 2.5346936860350002
121121
@test real(d[2]) 2.6159972444834976
122122
@test abs(imag(d[2])) 1.2917858749046127
@@ -413,6 +413,7 @@ end
413413
@test abs.(r.V[:, 1:n]'r0.V[:, 1:n]) I
414414
end
415415

416+
416417
# Regression test for #110.
417418
@testset "correct Krylov vector length check" begin
418419
m = 4
@@ -423,3 +424,31 @@ end
423424
@test_throws DimensionMismatch svds(a, nsv=1, v0 = ones(max(m, n)))
424425
@test_throws DimensionMismatch svds(a', nsv=1, v0 = ones(max(m, n)))
425426
end
427+
428+
@testset "ordering modes" begin
429+
N = 10
430+
nev = 4
431+
M = rand(N,N)
432+
433+
S = eigvals(M)
434+
435+
abs_imag = abs imag # ARPACK returns largest,smallest abs(imaginary) (complex pairs come together)
436+
437+
@testset "no shift-invert" begin
438+
for (which, sortby, rev) in [(:LM, abs, true), (:LR, real, true), (:LI, abs_imag, true),
439+
(:SM, abs, false), (:SR, real, false), (:SI, abs_imag, false)]
440+
d, _ = eigs(M, nev=nev, which=which)
441+
e = partialsort(S, 1:nev, by=sortby, rev=rev)
442+
@test sortby.(e) sortby.(d)
443+
end
444+
end
445+
446+
@testset "shift-invert" begin
447+
for (which, sortby, rev) in [(:LM, abs, true), (:LR, real, true), (:LI, abs_imag, true),
448+
(:SM, abs, false), (:SR, real, false), (:SI, abs_imag, false)]
449+
d, _ = eigs(M, nev=nev, which=which, sigma=0.0)
450+
e = S[partialsortperm(S, 1:nev, by=sortby inv, rev=rev)]
451+
@test sortby.(e) sortby.(d)
452+
end
453+
end
454+
end

0 commit comments

Comments
 (0)