Decompression of acyclic coloring with SparseMatrixCSC

Author: amontoison

src/decompression.jl

@@ -571,27 +571,40 @@ function decompress!(
             val = B[i, color[j]] - buffer_right_type[i]
             buffer_right_type[j] = buffer_right_type[j] + val
 
-            # A[i,j] = val
-            if in_triangle(i, j, uplo)
-                if uplo == :L
+            #! format: off
+            # A[i,j] is in the lower triangular part of A
+            if in_triangle(i, j, :L)
+                # uplo = :L or uplo = :F
+                # A[i,j] is stored at index_ij = (A.colptr[j+1] - offset_L) in A.nzval
+                if uplo != :U
                     nzind = A.colptr[j + 1] - lower_triangle_offsets[counter]
-                else
+                    nzA[nzind] = val
+                end
+
+                # uplo = :U or uplo = :F
+                # A[j,i] is stored at index_ji = (A.colptr[i] + offset_U) in A.nzval
+                if uplo != :L
+                    nzind = A.colptr[i] + upper_triangle_offsets[counter]
+                    nzA[nzind] = val
+                end
+
+            # A[i,j] is in the upper triangular part of A
+            else
+                # uplo = :U or uplo = :F
+                # A[i,j] is stored at index_ij = (A.colptr[j] + offset_U) in A.nzval
+                if uplo != :L
                     nzind = A.colptr[j] + upper_triangle_offsets[counter]
+                    nzA[nzind] = val
                 end
-                # Both values of `nzind` can be used if uplo = :F
-                nzA[nzind] = val
-            end
 
-            # A[j,i] = val
-            if in_triangle(j, i, uplo)
-                if uplo == :L
+                # uplo = :L or uplo = :F
+                # A[j,i] is stored at index_ji = (A.colptr[i+1] - offset_L) in A.nzval
+                if uplo != :U
                     nzind = A.colptr[i + 1] - lower_triangle_offsets[counter]
-                else
-                    nzind = A.colptr[i] + upper_triangle_offsets[counter]
+                    nzA[nzind] = val
                 end
-                # Both values of `nzind` can be used if uplo = :F
-                nzA[nzind] = val
             end
+            #! format: on
         end
     end
     return A

src/result.jl

@@ -347,7 +347,7 @@ function TreeSetColoringResult(
         end
 
         # continue until all leaves are treated
-        while queue_start ≤ queue_end
+        while queue_start <= queue_end
             leaf = queue[queue_start]
             queue_start += 1
 
@@ -373,29 +373,30 @@ function TreeSetColoringResult(
                     j = neighbor
                     col_i = view(rv, nzrange(S, i))
                     col_j = view(rv, nzrange(S, j))
-                    index_ij = S.colptr[j] - 1 + searchsortedfirst(col_j, i)  # S.nzval[index_ij] = S[i,j]
-                    @assert S.nzval[index_ij] == S[i,j]
-                    index_ji = S.colptr[i] - 1 + searchsortedfirst(col_i, j)  # S.nzval[index_ji] = S[j,i]
-                    @assert S.nzval[index_ji] == S[j,i]
                     index_offsets += 1
 
+                    #! format: off
+                    # S[i,j] is in the lower triangular part of S
                     if in_triangle(i, j, :L)
                         # uplo = :L or uplo = :F
-                        # A[i,j] is stored at index_ij = (A.colptr[j+1] - offset_L) in A.nzval
-                        lower_triangle_offsets[index_offsets] = S.colptr[j + 1] - index_ij
+                        # S[i,j] is stored at index_ij = (S.colptr[j+1] - offset_L) in S.nzval
+                        lower_triangle_offsets[index_offsets] = length(col_j) - searchsortedfirst(col_j, i) + 1
 
                         # uplo = :U or uplo = :F
-                        # A[j,i] is stored at index_ji = (A.colptr[i] + offset_U) in A.nzval
-                        upper_triangle_offsets[index_offsets] = index_ji - S.colptr[i]
-                    else
-                        # uplo = :L or uplo = :F
-                        # A[j,i] is stored at index_ji = (A.colptr[i+1] - offset_L) in A.nzval
-                        lower_triangle_offsets[index_offsets] = S.colptr[i + 1] - index_ji
+                        # S[j,i] is stored at index_ji = (S.colptr[i] + offset_U) in S.nzval
+                        upper_triangle_offsets[index_offsets] = searchsortedfirst(col_i, j)::Int - 1
 
+                    # S[i,j] is in the upper triangular part of S
+                    else
                         # uplo = :U or uplo = :F
-                        # A[i,j] is stored at index_ij = (A.colptr[j] + offset_U) in A.nzval
-                        upper_triangle_offsets[index_offsets] = index_ij - S.colptr[j]
+                        # S[i,j] is stored at index_ij = (S.colptr[j] + offset_U) in S.nzval
+                        upper_triangle_offsets[index_offsets] = searchsortedfirst(col_j, i)::Int - 1
+
+                        # uplo = :L or uplo = :F
+                        # S[j,i] is stored at index_ji = (S.colptr[i+1] - offset_L) in S.nzval
+                        lower_triangle_offsets[index_offsets] = length(col_i) - searchsortedfirst(col_i, j) + 1
                     end
+                    #! format: on
                 end
             end
         end