Add symmetric scaling - Githubissues

@michel2323 @frapac @sshin23 I discussed with Mike about symmetric scaling and I quickly implemented a Julia version of his MATLAB code:

function scaleSID(A::SparseMatrixCSC, iprint::Bool=false)
    n, m = size(A)
    if iprint
        println("\n scaleSID: Symmetric scaling of SID matrix.   n = $n   nnz = $(nnz(A))")
    end
    if n != m
        println("A must be square")
        return nothing, nothing
    end

    scale = maximum(abs.(A), dims=1)[:]
    Iz = findall(scale .== 0)
    lenz = length(Iz)
    if lenz > 0
        println("A contains $lenz empty columns")
        scale[Iz] .= 1
    end

    d = diag(A)            # dense column vector
    scale .= 1.0 ./ scale  # dense column vector
    d .= scale .* d        # dense diagonal of DAD
    scale .= sqrt.(scale)  # dense elements of D
    D = Diagonal(scale)    # sparse diagonal D
    AL = tril(A, -1)
    AL2 = copy(AL)
    mul!(AL2, AL, D)
    mul!(AL, D, AL2)
    DAD = AL + AL' + spdiagm(0 => d)

    if iprint
        jmin = argmin(scale)
        jmax = argmax(scale)
        smin = scale[jmin]
        smax = scale[jmax]
        println("\n\n  Min scale                     Max scale")
        println("  Col $jmin $smin    Col $jmax $smax")
    end

    return DAD, D
end

I tested with Krylov methods and it's quite efficient! It should not be hard to do a GPU version and add it in MadNLP.jl.

JuliaSmoothOptimizers / KrylovPreconditioners.jl

Add symmetric scaling #46