Unexported (Private) functions.

An approximate priority queue. Items are bundled together into doubly-linked lists with all approximately the same key. minlist is the min list we know to be non-empty. It should always be a lower bound. keyMap maps keys to lists

ijv = IJV(A::SparseMatrixCSC)

Convert a sparse matrix to an IJV.

LDLinv contains the information needed to solve the Laplacian systems. It does it by applying Linv, then Dinv, then Linv (transpose). But, it is specially constructed for this particular solver. It does not explicitly make the matrix triangular. Rather, col[i] is the name of the ith col to be eliminated

LLmatp is the data structure used to maintain the matrix during elimination. It stores the elements in each column in a singly linked list (only next ptrs) Each element is an LLp (linked list pointer). The head of each column is pointed to by cols.

We probably can get rid of degs - as it is only used to store initial degrees.

LLp elements are all in the same column. row tells us the row, and val is the entry. val is set to zero for some edges that we should remove. next gives the next in the column. It points to itself to terminate. reverse is the index into lles of the other copy of this edge, since every edge is stored twice as we do not know the order of elimination in advance.

Decrement the key of element i
This could crash if i exceeds the maxkey

Increment the key of element i
This could crash if i exceeds the maxkey

This variation of approxChol creates a cholesky factor to do the elimination. It has not yet been optimized, and does not yet make the cholesky factor lower triangular

B = augmentTreeOpt{Tv,Ti}(tree, A, params)

Takes as input a tree and an adjacency matrix of a graph. It then computes the stretch of every edge of the graph wrt the tree. It uses cholmod to decide how many edge to add back, shooting for nnzLfac times n entries in the factored augmented tree, with a number of flops to factor equal to nnz(a)*flopsfac. The edges to add back are then choen at random.

Apply the ith solver on the ith component

ijv = complete_bipartite_graph_ijv(n)

ijv = complete_graph_ijv(n)

ijv = empty_graph_ijv(n)

Add a new vertex to a with weights to the other vertices corresponding to diagonal surplus weight.

This is an efficient way of writing [a d; d' 0]

b = firstn(a::IJV, n::Integer)

Only keep the first n vertices of a.

la = forceLap(a)

Create a Laplacian matrix from an adjacency matrix. If the input looks like a Laplacian, throw a warning and convert it.

graph = generalizedNecklace(A, H, k::Int64)

Constructs a generalized necklace graph starting with two graphs A and H. The resulting new graph will be constructed by expanding each vertex in H to an instance of A. k random edges will be generated between components. Thus, the resulting graph may have weighted edges.

initDictCol!(dic, name, typ)

For a dictionary in which each key indexes an array. If dic does not contain an entry of name, create with set to Array(typ,0).

f = lapWrapComponents(solver, a::AbstractArray; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[], params...)

Applies a Laplacian solver that satisfies our interface to each connected component of the graph with adjacency matrix a. Passes kwargs on the solver.

f = lapWrapConnected(sddmSolver, a::AbstractMatrix; kwargs...)

Applies a sddmSolver to the Laplacian of the adjacency matrix a of a connected graph. Passes on kwargs to the solver. sddmSolver should be a solver that obeys the interface.

L = ldli2Chol(ldli)

This produces a matrix L so that L L^T approximate the original Laplacians. It is not quite a Cholesky factor, because it is off by a perm (and the all-1s vector orthogonality.

ijv = path_graph_ijv(n::Int64)

Print a column in an LLMatOrd matrix. This is here for diagnostics.

Print a column in an LLmatp matrix. This is here for diagnostics.

a = pure_random_ijv(n::Integer; verbose=false, prefix="", ver=Vcur)

Chooses among pathgraph, ringgraph, gridgraph, completebinarytree, randgenring, growngraph and ErdosRenyiClusterFix. It can produce a disconnected graph. For code that always produces a connected graph (and is the same as with Julia v0.6, use purerandomijv_v6)

ret is the answer returned by a speed test. This pushed it into the dictionary on which we are storing the tests.

a = rand_regular_bipartite(n,k)

Random k-regular bipartite graph between two sets of n vertices. No repeat edges, so can take a long time to build of k is close to n.

Returns a (possibly) asymmetric matrix that contains the upper-right block.

ijv = ring_graph_ijv(n)

ind = sampleByWeight(wt; ver=Vcur)

sample an index with probability proportional to its weight given here

f = sddmWrapLap(lapSolver, sddm::AbstractArray; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[], params...)

Uses a lapSolver to solve systems of linear equations in sddm matrices.

Given a set of integers, set between 1 and n, return a sorted version of them

ijv = star_graph_ijv(n::Int64)

testZeroDiag(a)

Returns true if a has zero diagonal, false otherwise

Compute the vector of depths in a tree that is in DFS order, with the root at the first position, and the leaves at the end

wted = uniformWeight(unwted)

Put a uniform [0,1] weight on every edge. This is an example of how to use mapweight.

f = wrapCapture(solver::Function, mats, rhss; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[], params...)

This wraps a solver so that we can capture all the matrices that it solves and all the right-hand-sides. Those are pushed into the arrays mats and rhss. For example

julia> mats = []
julia> rhss = []
julia> solver = wrapCapture(approxchol_lap, mats, rhss)
julia> a = chimera(10)
julia> f = solver(a);
julia> size(mats[1])
(10,10)
julia> b = randn(10)
julia> x = f(b);
julia> rhss
1-element Array{Any,1}:
 [0.404962,-0.827718,0.704616,-0.403223,0.204891,-0.505589,0.907015,1.90266,-0.438115,0.0464351]

f = wrapCaptureRhs(sola::Function, rhss; tol::Real=1e-6, maxits=Inf, maxtime=Inf, verbose=false, pcgIts=Int[], params...)

Captures all the right-hand-sides that are passed to the solver sola. It pushes them into an array called rhhs. For example

julia> rhss = []
julia> a = wted_chimera(100)
julia> sola = approxchol_lap(a)
julia> wrappedSolver = wrapCaptureRhs(sola,rhss)
julia> b = randn(100)
julia> x = wrappedSolver(b,verbose=true)

PCG BLAS stopped after: 0.0 seconds and 11 iterations with relative error 3.160275810360986e-7.

julia> length(rhss[1])

100

solveA = wrapInterface(solver::Function, A::AbstractMatrix; tol, maxits, maxtime, verbose, pcgIts=Int[],params...)
solverConstructor = wrapInterface(A::AbstractMatrix; tol, maxits, maxtime, verbose, pcgIts=Int[],params...)

Returns a function that discards tol, maxits, maxtime and verbose, sets pcgIts to 0 (because it might not be using pcg), and passes whatever params are left to the solver.

Examples

julia> a = randn(5,5);
julia> a = a * a';
julia> solvea = wrapInterface(X->cholesky(X,Val(true)), a, maxits=100, verbose=true);
julia> b = randn(5,1);
julia> norm(a*solvea(b, verbose=false)-b)
1.575705319704736e-14

julia> f = wrapInterface(X->cholesky(X,Val(true)))
julia> solvea = f(a, maxits=1000, maxtime = 1)
julia> norm(a*solvea(b, verbose=false, maxtime = 10)-b)
1.575705319704736e-14