RFC: the generalized eigenvalue problem (#6758)

base/linalg/arnoldi.jl

@@ -2,51 +2,67 @@ 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")
     else
         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

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)

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``.

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