Don't throw in eigs on semidefinite Bs for generalized eigenproblems. (…

…#17873)

This makes it possible to solve problems with semidefinite B via explicit
shift.
(cherry picked from commit f95b8b1)

base/linalg/arnoldi.jl

@@ -160,7 +160,7 @@ function _eigs(A, B;
     T = eltype(A)
     iscmplx = T <: Complex
     isgeneral = B !== I
-    sym = issymmetric(A) && !iscmplx
+    sym = issymmetric(A) && issymmetric(B) && !iscmplx
     nevmax=sym ? n-1 : n-2
     if nevmax <= 0
         throw(ArgumentError("Input matrix A is too small. Use eigfact instead."))
@@ -178,9 +178,6 @@ function _eigs(A, B;
         ncv = ncvmin
     end
     ncv = BlasInt(min(ncv, n))
-    if isgeneral && !isposdef(B)
-        throw(PosDefException(0))
-    end
     bmat = isgeneral ? "G" : "I"
     isshift = sigma !== nothing
 
@@ -251,7 +248,7 @@ function _eigs(A, B;
             solveSI = x->x
         else                #    Shift-invert mode
             mode       = 3
-            F = factorize(sigma==zero(T) ? A : A - UniformScaling(sigma))
+            F = factorize(A - UniformScaling(sigma))
             solveSI = x -> F \ x
         end
     else                    # Generalized eigenproblem
@@ -262,7 +259,7 @@ function _eigs(A, B;
             solveSI = x -> F \ x
         else                #    Shift-invert mode
             mode       = 3
-            F = factorize(sigma==zero(T) ? A : A-sigma*B)
+            F = factorize(A - sigma*B)
             solveSI = x -> F \ x
         end
     end

test/linalg/arnoldi.jl

@@ -72,7 +72,6 @@ let
             @test_throws ArgumentError eigs(a+a.',which=:LI)
             @test_throws ArgumentError eigs(a,sigma=rand(Complex64))
         end
-        @test_throws Base.LinAlg.PosDefException eigs(a,b)
     end
 end
 
@@ -232,3 +231,12 @@ svds(rand(1:10, 10, 8))
 @test_throws MethodError eigs(big(rand(1:10, 10, 10)))
 @test_throws MethodError eigs(big(rand(1:10, 10, 10)), rand(1:10, 10, 10))
 @test_throws MethodError svds(big(rand(1:10, 10, 8)))
+
+# Symmetric generalized with singular B
+let
+    n = 10
+    k = 3
+    A = randn(n,n); A = A'A
+    B = randn(n,k);  B = B*B'
+    @test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(eigvals(A, B)[1:k])
+end