diff --git a/Project.toml b/Project.toml index 0562463..5fc119d 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "1.4.1" [deps] Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97" +KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" @@ -14,9 +15,10 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] julia = "1" +Arpack = "0.3.2" +KrylovKit = "0.5.2" MathOptInterface = "0.9.14" TimerOutputs = "0.5.0" -Arpack = "0.3.2" [extras] diff --git a/src/ProxSDP.jl b/src/ProxSDP.jl index fb838f7..bb3bf92 100644 --- a/src/ProxSDP.jl +++ b/src/ProxSDP.jl @@ -1,13 +1,15 @@ module ProxSDP + using Arpack + using KrylovKit using MathOptInterface using TimerOutputs - using Arpack + using Printf using SparseArrays using LinearAlgebra - import Random + import LinearAlgebra: BlasInt include("structs.jl") diff --git a/src/eigsolver.jl b/src/eigsolver.jl index 38a06a1..7d1e773 100644 --- a/src/eigsolver.jl +++ b/src/eigsolver.jl @@ -729,22 +729,47 @@ end function eig!(arc::ARPACKAlloc, A::Symmetric{T,Matrix{T}}, nev::Integer, opt::Options)::Nothing where T up_ncv = false - for i in 1:1 - # Initialize parameters and do memory allocation - @timeit "update_arc" _update_arc!(arc, A, nev, opt, up_ncv)::Nothing + if opt.eigsolver == 1 + for i in 1:1 + # Initialize parameters and do memory allocation + @timeit "update_arc" _update_arc!(arc, A, nev, opt, up_ncv)::Nothing + + # Top level reverse communication interface to solve real double precision symmetric problems. + @timeit "saupd" _saupd!(arc)::Nothing + + if arc.nev > arc.iparam[5] + up_ncv = true + else + break + end + end - # Top level reverse communication interface to solve real double precision symmetric problems. - @timeit "saupd" _saupd!(arc)::Nothing + # Post processing routine (eigenvalues and eigenvectors purification) + @timeit "seupd" _seupd!(arc)::Nothing + else + @timeit "update_arc" _update_arc!(arc, A, nev, opt, up_ncv)::Nothing + @timeit "krylovkit" begin + vals, vecs, info = eigsolve( + A, arc.resid, arc.nev, :LR,Lanczos(krylovdim = arc.ncv, tol = 1e-12, maxiter=100) + ) - if arc.nev > arc.iparam[5] - up_ncv = true - else - break + if info.converged == 0 + @show arc.nev, info + end + fill!(arc.d, 0.0) + fill!(arc.v, 0.0) + for i in 1:min(arc.nev, info.converged) + arc.d[i] = vals[i] + for j in 1:arc.n + arc.v[j,i] = vecs[i][j] + end + end + arc.converged = true + arc.arpackerror = false + # if info.converged == 1 + # arc.converged = false + # end end end - - # Post processing routine (eigenvalues and eigenvectors purification) - @timeit "seupd" _seupd!(arc)::Nothing - return nothing end \ No newline at end of file diff --git a/src/prox_operators.jl b/src/prox_operators.jl index c05a263..b35cb50 100644 --- a/src/prox_operators.jl +++ b/src/prox_operators.jl @@ -27,15 +27,19 @@ function psd_projection!(v::Vector{Float64}, a::AuxiliaryData, cones::ConicSets, a.m[idx][1] = max(0., a.m[idx][1]) p.min_eig[idx] = a.m[idx][1] - elseif !opt.full_eig_decomp && p.target_rank[idx] <= opt.max_target_rank_krylov_eigs && sdp.sq_side > opt.min_size_krylov_eigs + elseif !opt.full_eig_decomp && + p.target_rank[idx] <= opt.max_target_rank_krylov_eigs && + sdp.sq_side > opt.min_size_krylov_eigs && + mod(p.iter, opt.full_eig_freq) > opt.full_eig_len # for full from time to time @timeit "eigs" begin eig!(arc_list[idx], a.m[idx], p.target_rank[idx], opt) + if hasconverged(arc_list[idx]) fill!(a.m[idx].data, 0.) for i in 1:p.target_rank[idx] if unsafe_getvalues(arc_list[idx])[i] > 0. p.current_rank[idx] += 1 - vec = unsafe_getvectors(arc_list[idx])[:, i] + vec = view(unsafe_getvectors(arc_list[idx]), :, i) LinearAlgebra.BLAS.gemm!('N', 'T', unsafe_getvalues(arc_list[idx])[i], vec, vec, 1., a.m[idx].data) end end @@ -74,7 +78,8 @@ function full_eig!(a::AuxiliaryData, idx::Int, opt::Options, p::Params)::Nothing fill!(a.m[idx].data, 0.) for i in 1:length(fact.values) if fact.values[i] > 0. - LinearAlgebra.BLAS.gemm!('N', 'T', fact.values[i], fact.vectors[:, i], fact.vectors[:, i], 1., a.m[idx].data) + v = view(fact.vectors, :, i) + LinearAlgebra.BLAS.gemm!('N', 'T', fact.values[i], v, v, 1., a.m[idx].data) if fact.values[i] > opt.tol_psd p.current_rank[idx] += 1 end diff --git a/src/structs.jl b/src/structs.jl index 087412f..ecb7234 100644 --- a/src/structs.jl +++ b/src/structs.jl @@ -51,9 +51,9 @@ Base.@kwdef mutable struct Options # Linesearch parameters max_linsearch_steps::Int = 5000 - delta::Float64 = .999 + delta::Float64 = .9999 initial_theta::Float64 = 1. - linsearch_decay::Float64 = .9 + linsearch_decay::Float64 = .75 # Spectral decomposition parameters full_eig_decomp::Bool = false @@ -63,6 +63,13 @@ Base.@kwdef mutable struct Options rank_increment::Int = 1 # 0=multiply, 1 = add rank_increment_factor::Int = 1 # 0 multiply, 1 = add + # eigsolver selection + #= + 1: Arpack [dsaupd] (tipically non-deterministic) + 2: KrylovKit [eigsolve/Lanczos] (DEFAULT) + =# + eigsolver::Int = 2 + # Arpack arpack_tol::Float64 = 1e-10 #= @@ -73,7 +80,7 @@ Base.@kwdef mutable struct Options =# arpack_resid_init::Int = 3 arpack_resid_seed::Int = 1234 - arpack_reset_resid::Bool = true # true for determinism + arpack_reset_resid::Bool = false # true for determinism arpack_max_iter::Int = 10_000 arpack_min_lanczos::Int = 25 # larger is more stable to converge and more deterministic @@ -83,6 +90,9 @@ Base.@kwdef mutable struct Options reduce_rank::Bool = false rank_slack::Int = 3 + full_eig_freq::Int = 10000000 + full_eig_len::Int = 0 + # equilibration parameters equilibration::Bool = true equilibration_iters::Int = 1000