From e9939e6b905b620d9fbd40780d9627a8c7d98f6e Mon Sep 17 00:00:00 2001 From: "Viral B. Shah" Date: Sat, 10 May 2014 16:24:59 +0530 Subject: [PATCH] Solve generalized eigen-problem with eigs (Fix #6758) Remove complexOP keyword argument, as I believe it was not used correctly. Update docs and add tests. --- base/linalg/arnoldi.jl | 72 ++++++++++++++++++++++++++---------------- base/linalg/arpack.jl | 50 +++++++++++++++++++++++++++-- doc/stdlib/linalg.rst | 5 ++- test/arpack.jl | 35 ++++++++++++++++---- 4 files changed, 121 insertions(+), 41 deletions(-) diff --git a/base/linalg/arnoldi.jl b/base/linalg/arnoldi.jl index a3ded63eb69ae..17a843aa1e758 100644 --- a/base/linalg/arnoldi.jl +++ b/base/linalg/arnoldi.jl @@ -2,18 +2,25 @@ using .ARPACK ## eigs -function eigs(A;nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM", - tol=0.0, maxiter::Integer=1000, sigma=nothing, v0::Vector=zeros((0,)), - ritzvec::Bool=true, complexOP::Bool=false) +eigs(A; args...) = eigs(A, :Identity; args...) + +function eigs(A, B; + nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM", + tol=0.0, maxiter::Integer=1000, sigma=nothing, v0::Vector=zeros((0,)), + ritzvec::Bool=true) + n = chksquare(A) (n <= 6) && (nev = min(n-1, nev)) ncv = blas_int(min(max(2*nev+2, ncv), n)) sym = issym(A) T = eltype(A) - cmplx = T<:Complex - bmat = "I" - + iscmplx = T<:Complex + isgeneral = B !== :Identity + bmat = isgeneral ? "G" : "I" + isshift = sigma !== nothing + sigma = isshift ? sigma : zero(T) + if !isempty(v0) length(v0)==n || throw(DimensionMismatch("")) eltype(v0)==T || error("Starting vector must have eltype $T") @@ -21,32 +28,41 @@ function eigs(A;nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM", v0=zeros(T,(0,)) end - if sigma===nothing # normal iteration - mode = 1 - sigma = 0 - linop(x) = A * x - else # inverse iteration with level shift - C = lufact(sigma==0 ? A : A - sigma*eye(A)) - if cmplx - mode = 3 - linop(x) = C\x - else - if !complexOP - mode = 3 - linop(x) = real(C\x) - else - mode = 4 - linop(x) = imag(C\x) - end - end + # Refer to ex-*.doc files in ARPACK/DOCUMENTS for calling sequence + matvecA(x) = A * x + if !isgeneral # Standard problem + matvecB(x) = x + if !isshift # Regular mode + mode = 1 + solveSI(x) = x + else # Shift-invert mode + mode = 3 + solveSI(x) = lufact(sigma==0 ? A : A-sigma*eye(A)) \ x + end + else # Generalized eigen problem + matvecB(x) = B * x + if !isshift # Regular inverse mode + mode = 2 + solveSI(x) = lufact(B) \ x + else # Shift-invert mode + mode = 3 + solveSI(x) = lufact(sigma==0 ? A : A-sigma*B) \ x + end end - + # Compute the Ritz values and Ritz vectors (resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, TOL) = - ARPACK.aupd_wrapper(T, linop, n, sym, cmplx, bmat, nev, ncv, which, tol, maxiter, mode, v0) + ARPACK.aupd_wrapper(T, matvecA, matvecB, solveSI, n, sym, iscmplx, bmat, nev, ncv, which, tol, maxiter, mode, v0) # Postprocessing to get eigenvalues and eigenvectors - ARPACK.eupd_wrapper(T, n, sym, cmplx, bmat, nev, which, ritzvec, TOL, - resid, ncv, v, ldv, sigma, iparam, ipntr, workd, workl, lworkl, rwork) + if ritzvec + (eval, evec, nconv, niter, nmult, resid) = + ARPACK.eupd_wrapper(T, n, sym, iscmplx, bmat, nev, which, ritzvec, TOL, + resid, ncv, v, ldv, sigma, iparam, ipntr, workd, workl, lworkl, rwork) + else + (eval, nconv, niter, nmult, resid) = + ARPACK.eupd_wrapper(T, n, sym, iscmplx, bmat, nev, which, ritzvec, TOL, + resid, ncv, v, ldv, sigma, iparam, ipntr, workd, workl, lworkl, rwork) + end end diff --git a/base/linalg/arpack.jl b/base/linalg/arpack.jl index b7dcf839d186a..ba54f1242b55e 100644 --- a/base/linalg/arpack.jl +++ b/base/linalg/arpack.jl @@ -4,7 +4,7 @@ import ..LinAlg: BlasInt, blas_int, ARPACKException ## aupd and eupd wrappers -function aupd_wrapper(T, linop::Function, n::Integer, +function aupd_wrapper(T, matvecA::Function, matvecB::Function, solveSI::Function, n::Integer, sym::Bool, cmplx::Bool, bmat::ASCIIString, nev::Integer, ncv::Integer, which::ASCIIString, tol::Real, maxiter::Integer, mode::Integer, v0::Vector) @@ -47,11 +47,55 @@ function aupd_wrapper(T, linop::Function, n::Integer, naupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n, iparam, ipntr, workd, workl, lworkl, info) end + + # Check for warnings and errors + # Refer to ex-*.doc files in ARPACK/DOCUMENTS for calling sequence if info[1] == 3; warn("try eigs/svds with a larger value for ncv"); end if info[1] == 1; warn("maximum number of iterations reached; check nconv for number of converged eigenvalues"); end if info[1] < 0; throw(ARPACKException(info[1])); end - if (ido[1] != -1 && ido[1] != 1); break; end - workd[ipntr[2]+zernm1] = linop(getindex(workd, ipntr[1]+zernm1)) + + load_idx = ipntr[1]+zernm1 + store_idx = ipntr[2]+zernm1 + x = workd[load_idx] + if ido[1] == -1 + if mode == 1 + workd[store_idx] = matvecA(x) + elseif mode == 2 + if sym + temp = matvecA(x) + workd[load_idx] = temp # overwrite as per Remark 5 in dsaupd.f + workd[store_idx] = solveSI(temp) + else + workd[store_idx] = solveSI(matvecA(x)) + end + elseif mode == 3 + if bmat == "I" + workd[store_idx] = solveSI(x) + elseif bmat == "G" + workd[store_idx] = solveSI(matvecB(x)) + end + end + elseif ido[1] == 1 + if mode == 1 + workd[store_idx] = matvecA(x) + elseif mode == 2 + if sym + temp = matvecA(x) + workd[load_idx] = temp + workd[store_idx] = solveSI(temp) + else + workd[store_idx] = solveSI(matvecA(x)) + end + elseif mode == 3 + workd[store_idx] = solveSI(workd[ipntr[3]+zernm1]) + end + elseif ido[1] == 2 + workd[store_idx] = matvecB(x) + elseif ido[1] == 99 + break + else + error("Internal ARPACK error") + end end return (resid, v, n, iparam, ipntr, workd, workl, lworkl, rwork, TOL) diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 1b482f5d45c66..b6a4af7e6aa52 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -494,9 +494,9 @@ Linear algebra functions in Julia are largely implemented by calling functions f The conjugate transposition operator (``'``). -.. function:: eigs(A; nev=6, which="LM", tol=0.0, maxiter=1000, sigma=nothing, ritzvec=true, op_part=:real,v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid) +.. function:: eigs(A, [B,]; nev=6, which="LM", tol=0.0, maxiter=1000, sigma=nothing, ritzvec=true, v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid) - ``eigs`` computes eigenvalues ``d`` of ``A`` using Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively. The following keyword arguments are supported: + ``eigs`` computes eigenvalues ``d`` of ``A`` using Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively. If ``B`` is provided, the generalized eigen-problem is solved. The following keyword arguments are supported: * ``nev``: Number of eigenvalues * ``which``: type of eigenvalues to compute. See the note below. @@ -518,7 +518,6 @@ Linear algebra functions in Julia are largely implemented by calling functions f * ``maxiter``: Maximum number of iterations * ``sigma``: Specifies the level shift used in inverse iteration. If ``nothing`` (default), defaults to ordinary (forward) iterations. Otherwise, find eigenvalues close to ``sigma`` using shift and invert iterations. * ``ritzvec``: Returns the Ritz vectors ``v`` (eigenvectors) if ``true`` - * ``op_part``: which part of linear operator to use for real ``A`` (``:real``, ``:imag``) * ``v0``: starting vector from which to start the iterations ``eigs`` returns the ``nev`` requested eigenvalues in ``d``, the corresponding Ritz vectors ``v`` (only if ``ritzvec=true``), the number of converged eigenvalues ``nconv``, the number of iterations ``niter`` and the number of matrix vector multiplications ``nmult``, as well as the final residual vector ``resid``. diff --git a/test/arpack.jl b/test/arpack.jl index 7ed33135843e4..7df010d532bca 100644 --- a/test/arpack.jl +++ b/test/arpack.jl @@ -1,33 +1,54 @@ - begin - local n,a,asym,d,v +begin + srand(1234) + local n,a,asym,b,bsym,d,v n = 10 areal = randn(n,n) + breal = randn(n,n) acmplx = complex(randn(n,n), randn(n,n)) + bcmplx = complex(randn(n,n), randn(n,n)) for elty in (Float32, Float64, Complex64, Complex128) if elty == Complex64 || elty == Complex128 a = acmplx + b = bcmplx else a = areal + b = breal end a = convert(Matrix{elty}, a) asym = a' + a # symmetric indefinite apd = a'*a # symmetric positive-definite + b = convert(Matrix{elty}, b) + bsym = b' + b + bpd = b'*b + + if elty == Complex64 || elty == Float32 + testtol = 1e-3 + else + testtol = 1e-6 + end + (d,v) = eigs(a, nev=3) @test_approx_eq a*v[:,2] d[2]*v[:,2] + (d,v) = eigs(a, b, nev=3, tol=1e-8) + @test_approx_eq_eps a*v[:,2] d[2]*b*v[:,2] testtol (d,v) = eigs(asym, nev=3) @test_approx_eq asym*v[:,1] d[1]*v[:,1] -# @test_approx_eq eigs(asym; nev=1, sigma=d[3])[1][1] d[3] + @test_approx_eq eigs(asym; nev=1, sigma=d[3])[1][1] d[3] (d,v) = eigs(apd, nev=3) @test_approx_eq apd*v[:,3] d[3]*v[:,3] -# @test_approx_eq eigs(apd; nev=1, sigma=d[3])[1][1] d[3] + @test_approx_eq eigs(apd; nev=1, sigma=d[3])[1][1] d[3] + (d,v) = eigs(apd, bpd, nev=3, tol=1e-8) + @test_approx_eq_eps apd*v[:,2] d[2]*bpd*v[:,2] testtol - # test (shift-and-)invert mode - (d,v) = eigs(apd, nev=3, sigma=0) - @test_approx_eq apd*v[:,3] d[3]*v[:,3] + # test (shift-and-)invert mode + (d,v) = eigs(apd, nev=3, sigma=0) + @test_approx_eq apd*v[:,3] d[3]*v[:,3] + (d,v) = eigs(apd, bpd, nev=3, sigma=0, tol=1e-8) + @test_approx_eq_eps apd*v[:,1] d[1]*bpd*v[:,1] testtol end end