View the Project on GitHub rjkyng/SDDM2023

  • Tutorial for SDDM solver experiments
  • Guide to testing other solvers
  • Examples and discussion of our APIs
  • Examples and discussion of our APIs

    On this page, we discuss how to use our APIs provided in Laplacians.jl and SDDM2023.git to reproduce our experiments.

    Generating problem instances

    In this section we present how to generate the matrices for each problem tested in our experiments, using the Laplacians.jl package. For the Poisson grid problem, it is not always possible to generate the matrix with the desired number of non-zero entries exactly. We solve for the number of variables exactly and round it to the closest integer.

    Running our solvers

    First, set up the approximate Cholesky algorithm and solver. Suppose you are solving a linear system in M. If M is strictly SDDM, then:

    solver = approxchol_sddm(M, verbose=false,params=ApproxCholParams(:deg), tol=1e-8)

    Otherwise, if M is a Laplacian, then you need to pass in the adjecency matrix corresponding to M:

    a, d = adj(M)

    a is the adjecency matrix and d is $M\mathbf{1}$ which is an all-zero vector in this case. Then set up the approximate Cholesky algorithm and the solver as follows:

    solver = approxchol_lap(a, verbose=false,params=ApproxCholParams(:deg), tol=1e-8)

    params allows you to control some parameters of the approximate Cholesky algorithm. In particular, :deg makes the approximate Cholesky algorithm eliminates vertices in an order that is greedy on the unweighted degree of the vertices. Note that this ordering is dynamic. You can pass in split and merge parameters. split controls how many copies are the original edges of the graph splitted into. By default the edges are not splitted (and multi-edges are not allowed). If you pass in the split parameter, it should be at least 1. merge controls how many copies of each edge do you want to keep during the algorithm. Only pass in merge parameter if you also pass in the split parameter and merge should be at least the value of split. By default, merge is set to infinity, in other words, if a split is specified (to be at least 1) while merge is not, then no merging is performed.

    To set up the approximate cholesky parameters to greedy ordering and no multi-edge:

    params = ApproxCholParams(:deg)

    To set up the approximate cholesky parameters to greedy ordering and split the edges into 2 copies and no merging:

    params = ApproxCholParams(:deg, 2)

    To set up the approximate cholesky parameters to greedy ordering and split the edges into 2 copies and merge all multi-edges up to at most 2 copies:

    params = ApproxCholParams(:deg, 2, 2)

    Once the solver is setup, you can pass in the right hand side b and solve the linear system:

    x = solver(b)


    We recommend you to benchmark the performances of our solvers through another interface that we provide at SDDM2023/julia-files/compareSolvers.jl. Again, you need to setup the solvers first.

    ac_deg_sddm = function(M; verbose=false, args...)
        approxchol_sddm(M; params=ApproxCholParams(:deg), verbose=verbose, args...)
    test_ac_deg_sddm = SolverTest(ac_deg_sddm, "ac_deg_sddm")
    ac_wdeg_sddm = function(M; verbose=false, args...)
        approxchol_sddm(M; params=ApproxCholParams(:wdeg), verbose=verbose, args...)
    test_ac_wdeg_sddm = SolverTest(ac_wdeg_sddm, "ac_wdeg_sddm")

    Then put them into an array and pass to our benchmarking function, along with a dictionary to store the benchmarking information and the linear system M (which is strictly SDD) and b:

    tests_sddm = [test_ac_deg_sddm, test_ac_wdeg_sddm]
    dic = Dict()
    x = testSddm(tests_sddm, dic, M, b; verbose=false, tol=1e-8, testName="test")

    testSddm() will store the benchmarking information of the solvers into dic. You can store benchmarking information from different linear systems into dic by passing dic to testSddm() repeatedly.

    Similar to before, if you want to benchmark the solvers on a Laplacian M, you need to use a slightly different interface:

    ac_deg_lap = function(a; verbose=false, args...)
        approxchol_lap(a; params=ApproxCholParams(:deg), verbose=verbose, args...)
    test_ac_deg_lap = SolverTest(ac_deg_lap, "ac_deg_lap")
    ac_wdeg_lap = function(a; verbose=false, args...)
        approxchol_lap(a; params=ApproxCholParams(:wdeg), verbose=verbose, args...)
    test_ac_wdeg_lap = SolverTest(ac_wdeg_lap, "ac_wdeg_lap")
    tests_lap = [test_ac_deg_lap, test_ac_wdeg_lap]
    a, _ = adj(M)
    dic = Dict()
    x = testLap(tests_lap, dic, a, b; verbose=false, tol=1e-8, testName="test")

    testSddm() and testLap() is also the interface for benchmarking other solvers. You do this simply as the following:

    x = testSddm(tests_sddm, dic_sddm, M, b; verbose=false, tol=1e-8, testName=tn, test_petsc_hypre=true,test_hypre=true, test_icc=true, test_cmg=true,test_lamg=true)

    This way, the benchmarking information for HyPre, PETSc (using HyPre BoomerAMG), MATLAB’s ICC, CMG (MATLAB version) and LAMG will be stored in dic as well.

    We strongly recommend that you do a warmup for testSddm() and testLap() over a smaller linear system before the benchmarking.

    Setting up and running other solvers

    The solvers that we benchmarked are CMG (MATLAB version), HyPre, PETSc, MATLAB’s ICC, and LAMG. Here we briefly introduce the setup of these solvers in our experiments.