Coloring for symmetric banded matrices

The code is quite short but it took me a long a time to understand how we can do the decompression directly from the peridic coloring.
We only need to solve bidiagonal systems for the decompression of each 2-colored tree.
It also allows to not require a buffer like in the acyclic coloring.

using BandedMatrices

T = Float64
n = 10
kd = 3
k = 2 * kd + 1

A = brand(T,n,n,kd,kd)
A = A + A'

color = Int[mod(j,kd+1)+1 for j=1:n]

B = zeros(T, n, kd+1)
for j = 1:n
    cj = color[j]
    B[:,cj] .+= A[:,j]
end

C = copy(A)
fill!(C, 0)

for i = 1:n
    ci = color[i]
    C[i,i] = B[i, ci]
end

for j = 1:kd
    for i = j+1:kd+1
        posi = i
        posj = j
        c = color[posi]
        C[posi,posj] = B[posj,c]
        C[posj,posi] = C[posi,posj]
        flag = true
        while flag
            if posi > posj
                posi, posj = posj, posi
            end
            if posi + kd + 1 <= n
                c = color[posi + kd + 1]
                C[posi+kd+1,posj] = B[posj,c] - C[posj,posi]
                C[posj,posi+kd+1] = C[posi+kd+1,posj]
                posi = posi + kd + 1
            else
                flag = false
            end
        end
    end
end