lobpcg_test2.jl

lobpcg_test2.jl

using IterativeSolvers    
    
function run_test(;N=1000, nev=10)    
    A = [4.71217, 4.59973, 4.44981, 3.96257, 4.59973, 4.44981, 4.71217, 3.58777, 3.10053,    
         3.73769, 4.93705, 2.95061, 2.08857, 2.35093, 3.73769, 3.43785, 2.20101, 2.08857,    
         3.10053, 3.43785, 2.95061, 3.58777, 4.93705, 3.58777, 2.72573, 2.98809, 4.37485,    
         4.93705, 2.57581, 1.33897, 1.22653, 2.23849, 4.37485, 2.68825, 1.07661, 0.589372,
         1.22653, 2.98809, 3.92509, 1.93865, 1.07661, 1.33897, 2.72573, 3.92509, 2.68825,    
         2.57581, 3.58777, 4.93705, 4.71217, 3.47533, 3.36289, 4.37485, 3.32541, 1.71377,    
         1.22653, 1.86369, 3.62525, 3.06305, 1.07661, 0.214573, 0.476932, 1.86369, 4.37485,
         3.92509, 1.56385, 0.327013, 0.214573, 1.22653, 3.36289, 3.17549, 1.56385, 1.07661,
         1.71377, 3.47533, 3.92509, 3.06305, 3.32541, 4.71217, 4.86209, 3.21297, 2.35093,
         2.61329, 4.00005, 4.56225, 2.20101, 0.964172, 0.851732, 1.86369, 4.00005, 2.31345,
         0.701812, 0.214573, 0.851732, 2.61329, 3.55029, 1.56385, 0.701812, 0.964172,
         2.35093, 4.86209, 3.55029, 2.31345, 2.20101, 3.21297, 4.56225, 4.59973, 4.48729,
         4.44981, 2.83817, 2.35093, 2.98809, 4.74965, 4.18745, 2.20101, 1.33897, 1.60133,
         2.98809, 2.68825, 1.45141, 1.33897, 2.35093, 4.48729, 4.29989, 2.68825, 2.20101,
         2.83817, 4.59973, 4.18745, 4.44981, 4.97453, 4.82461, 3.58777, 3.47533, 4.48729,    
         4.93705, 3.32541, 2.83817, 3.47533, 4.18745, 3.32541, 3.58777, 4.97453, 4.93705,    
         4.82461]

    reference = sort(A)[1:nev]
    println("Size of matrix:  $(length(A))")
    println("Reference evals: $(reference)")
    println("Separation:      $(sort(A)[nev + 1] - sort(A)[nev])")

    faillog = nothing
    worked=0
    for i in 1:N
            res = lobpcg(Diagonal(A), false, nev, log=true)
            if res.λ ≈ reference
                worked += 1
            else
                faillog = res
            end
    end
    println("Success rate: $(worked / N * 100)%")

    if faillog !== nothing
        println("")
        println("A failing trace:")
        println(faillog.trace)
    end
end