diff --git a/src/coloring/high_level.jl b/src/coloring/high_level.jl
index ffc66881..0c79473a 100644
--- a/src/coloring/high_level.jl
+++ b/src/coloring/high_level.jl
@@ -50,8 +50,8 @@ function matrix_colors(A::BandedMatrix)
 end
 
 function matrix_colors(A::BlockBandedMatrix)
-    u,l=blockbandwidths(A)
-    blockwidth=u+l+1
+    l,u=blockbandwidths(A)
+    blockwidth=l+u+1
     nblock=nblocks(A,2)
     cols=[blocksize(A,(1,i))[2] for i in 1:nblock]
     blockcolors=_cycle(1:blockwidth,nblock)
@@ -61,4 +61,20 @@ function matrix_colors(A::BlockBandedMatrix)
     startinds=[endinds[i]-ncolors[i]+1 for i in 1:blockwidth]
     colors=[(startinds[blockcolors[i]]:endinds[blockcolors[i]])[1:cols[i]] for i in 1:nblock]
     vcat(colors...)
+end
+
+function matrix_colors(A::BandedBlockBandedMatrix)
+    l,u=blockbandwidths(A)
+    lambda,mu=subblockbandwidths(A)
+    blockwidth=l+u+1
+    subblockwidth=lambda+mu+1
+    nblock=nblocks(A,2)
+    cols=[blocksize(A,(1,i))[2] for i in 1:nblock]
+    blockcolors=_cycle(1:blockwidth,nblock)
+    #the reserved number of colors of a block is the min of subblockwidth and the largest length of columns of blocks with the same block color
+    ncolors=[min(subblockwidth,maximum(cols[i:blockwidth:nblock])) for i in 1:blockwidth]
+    endinds=cumsum(ncolors)
+    startinds=[endinds[i]-ncolors[i]+1 for i in 1:blockwidth]
+    colors=[_cycle(startinds[blockcolors[i]]:endinds[blockcolors[i]],cols[i]) for i in 1:nblock]
+    vcat(colors...)
 end
\ No newline at end of file
diff --git a/test/test_specialmatrices.jl b/test/test_specialmatrices.jl
index 87c56685..72f0ce0e 100644
--- a/test/test_specialmatrices.jl
+++ b/test/test_specialmatrices.jl
@@ -17,6 +17,8 @@ symtridiagonal=SymTridiagonal(dense)
 banded=BandedMatrix(dense,(1,2))
 blockbanded1=BlockBandedMatrix(dense,([1,2,3,4],[4,3,2,1]),(1,0))
 blockbanded2=BlockBandedMatrix(dense,([4,3,2,1],[1,2,3,4]),(1,1))
+bandedblockbanded1=BandedBlockBandedMatrix(dense,([1,2,3,4],[4,3,2,1]),(1,0),(1,1))
+bandedblockbanded2=BandedBlockBandedMatrix(dense,([4,3,2,1],[1,2,3,4]),(1,1),(1,0))
 
 @test matrix_colors(dense)==1:n
 @test matrix_colors(uptri)==1:n
@@ -30,4 +32,29 @@ blockbanded2=BlockBandedMatrix(dense,([4,3,2,1],[1,2,3,4]),(1,1))
 
 @test matrix_colors(banded)==[1,2,3,4,1,2,3,4,1,2]
 @test matrix_colors(blockbanded1)==[1,2,3,4,5,6,7,1,2,5]
-@test matrix_colors(blockbanded2)==[1,5,6,7,8,9,1,2,3,4]
\ No newline at end of file
+@test matrix_colors(blockbanded2)==[1,5,6,7,8,9,1,2,3,4]
+@test matrix_colors(bandedblockbanded1)==[1,2,3,1,4,5,6,1,2,4]
+@test matrix_colors(bandedblockbanded2)==[1,3,4,5,6,5,1,2,1,2]
+
+function _testvalidity(A)
+    colorvec=matrix_colors(A)
+    ncolor=maximum(colorvec)
+    for color in 1:ncolor
+        subA=A[:,findall(x->x==color,colorvec)]
+        @test maximum(sum(subA,dims=2))<=1.0
+    end
+end
+
+_testvalidity(dense)
+_testvalidity(uptri)
+_testvalidity(lotri)
+_testvalidity(diagonal)
+_testvalidity(bidiagonalU)
+_testvalidity(bidiagonalL)
+_testvalidity(tridiagonal)
+_testvalidity(symtridiagonal)
+_testvalidity(banded)
+_testvalidity(blockbanded1)
+_testvalidity(blockbanded2)
+_testvalidity(bandedblockbanded1)
+_testvalidity(bandedblockbanded2)