Maximize test coverage and fix bugs

Author: mohamed82008

src/lobpcg.jl

@@ -445,11 +445,11 @@ LOBPCGIterator(A, largest::Bool, X, P=nothing, C=nothing) = LOBPCGIterator(A, no
     to a subspace; must overload `A_mul_B!`;
 """
 function LOBPCGIterator(A, B, largest::Bool, X, P=nothing, C=nothing)
-    constr! = Constraint(C, B, X, BWrapper())
     precond! = RPreconditioner(P, X)
-    return LOBPCGIterator(A, B, largest, X, constr!, precond!)
+    constr! = Constraint(C, B, X)
+    return LOBPCGIterator(A, B, largest, X, precond!, constr!)
 end
-function LOBPCGIterator(A, B, largest::Bool, X, constr!::Constraint, precond!::RPreconditioner)
+function LOBPCGIterator(A, B, largest::Bool, X, precond!::RPreconditioner, constr!::Constraint)
     T = eltype(X)
     nev = size(X, 2)
     if B isa Void
@@ -515,7 +515,7 @@ function LOBPCGIterator(A, B, largest::Bool, X, nev::Int, P=nothing, C=nothing)
     end
     constr! = Constraint(Y, BY, X, NotBWrapper())
     precond! = RPreconditioner(P, X)
-    return LOBPCGIterator(A, B, largest, X, constr!, precond!)
+    return LOBPCGIterator(A, B, largest, X, precond!, constr!)
 end
 
 function ortho_AB_mul_X!(blocks::Blocks, ortho!, A, B, bs=-1)
@@ -929,7 +929,7 @@ function lobpcg(A, B, largest::Bool, X0, nev::Int;
 
     sizeX = min(nev, sizeX)
     X = X0[:, 1:sizeX]
-    iterator = LOBPCGIterator(A, B, largest, X, nev, C, P)
+    iterator = LOBPCGIterator(A, B, largest, X, nev, P, C)
 
     r = EmptyLOBPCGResults(X, nev, tol, maxiter)
     rnext = lobpcg!(iterator, log=log, tol=tol, maxiter=maxiter, not_zeros=not_zeros)

test/lobpcg.jl

@@ -286,15 +286,15 @@ end
             end
         end
     end
-    @testset "nev = 3, block size = 2" begin
+    @testset "nev = 3, block size = $block_size" for block_size in (1, 2)
         n = 10
         @testset "Simple eigenvalue problem" begin
             @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128)
                 @testset "largest = $largest" for largest in (true, false)
                     A = rand(T, n, n)
                     A = A' * A + I
                     tol = eps(real(T))^0.4
-                    X0 = rand(T, n, 2)
+                    X0 = rand(T, n, block_size)
                     r = lobpcg(A, largest, X0, 3, tol=tol, maxiter=Inf, log=true)
                     λ, X = r.λ, r.X
                     @test max_err(A*X - X*diagm(λ)) ≤ tol
@@ -311,13 +311,53 @@ end
                     B = B' * B + I
                     tol = eps(real(T))^0.4
 
-                    X0 = rand(T, n, 2)
+                    X0 = rand(T, n, block_size)
                     r = lobpcg(A, B, largest, X0, 3, tol=tol, maxiter=Inf, log=true)
                     λ, X = r.λ, r.X
                     @test max_err(A*X - B*X*diagm(λ)) ≤ tol
                     @test all(isapprox.(Ac_mul_B(X, B*X), eye(3), atol=2*n*tol))
                 end
             end
         end
+        @testset "Constraint" begin
+            @testset "Simple eigenvalue problem" begin
+                @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128)
+                    @testset "largest = $largest" for largest in (true, false)
+                        A = rand(T, n, n)
+                        A = A' * A + I
+                        tol = √eps(real(T))
+                        r = lobpcg(A, largest, 1; tol=tol, maxiter=Inf, log=false)
+                        λ1, X1 = r.λ, r.X
+
+                        X0 = rand(T, n, block_size)
+                        r = lobpcg(A, largest, X0, 3, C=copy(r.X), tol=tol, maxiter=Inf, log=true)
+                        λ2, X2 = r.λ, r.X
+                        @test max_err(A*X2 - X2*diagm(λ2)) ≤ tol
+                        @test all(isapprox.(Ac_mul_B(X2, X2), eye(3), atol=2*n*tol))
+                        @test all(isapprox.(real(Ac_mul_B(X1, X2)), 0, atol=2*n*tol))
+                    end
+                end
+            end
+            @testset "Generalized eigenvalue problem" begin
+                @testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128)
+                    @testset "largest = $largest" for largest in (true, false)
+                        A = rand(T, n, n)
+                        A = A' * A + 2I
+                        B = rand(T, n, n)
+                        B = B' * B + 2I
+                        tol = eps(real(T))^0.4
+                        r = lobpcg(A, B, largest, 1; tol=tol, maxiter=Inf, log=false)
+                        λ1, X1 = r.λ, r.X
+
+                        X0 = rand(T, n, block_size)
+                        r = lobpcg(A, B, largest, X0, 2, C=copy(r.X), tol=tol, maxiter=Inf, log=true)
+                        λ2, X2 = r.λ, r.X
+                        @test max_err(A*X2 - B*X2*diagm(λ2)) ≤ tol
+                        @test all(isapprox.(Ac_mul_B(X2, B*X2), eye(2), atol=2*n*tol))
+                        @test all(isapprox.(real(Ac_mul_B(X1, B*X2)), 0, atol=2*n*tol))
+                    end
+                end
+            end
+        end        
     end
 end