diff --git a/Project.toml b/Project.toml
index 8b44018..0937b0c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -26,11 +26,9 @@ Graphs = "1"
 LightGraphs = "1"
 METIS_jll = "5.1"
 SimpleWeightedGraphs = "1"
-LinearAlgebra = "1"
 julia = "1.6"
 
 [extras]
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/src/Metis.jl b/src/Metis.jl
index f650f91..eae3e35 100644
--- a/src/Metis.jl
+++ b/src/Metis.jl
@@ -3,7 +3,7 @@
 module Metis
 
 using SparseArrays
-using LinearAlgebra: ishermitian, Hermitian, istril, istriu
+using LinearAlgebra: ishermitian, Hermitian, Symmetric
 using METIS_jll: libmetis
 
 # Metis C API: Clang.jl auto-generated bindings and some manual methods
@@ -11,9 +11,6 @@ include("LibMetis.jl")
 using .LibMetis
 using .LibMetis: idx_t, @check
 
-# Various operations for sparse matrices
-include("utils.jl")
-
 # Global options array -- should probably do something better...
 const options = fill(Cint(-1), METIS_NOPTIONS)
 options[Int(METIS_OPTION_NUMBERING)+1] = 1
@@ -75,57 +72,86 @@ function graph(G::SparseMatrixCSC; weights::Bool=false, check_hermitian::Bool=tr
     return Graph(idx_t(N), xadj, adjncy, vwgt, adjwgt)
 end
 
+const HermOrSymCSC{Tv, Ti} = Union{
+    Hermitian{Tv, SparseMatrixCSC{Tv, Ti}}, Symmetric{Tv, SparseMatrixCSC{Tv, Ti}},
+}
+
 """
-    Metis.graph(G::Hermitian{<:T,SparseMatrixCSC{T,I}}; weights::Bool=false) where {T,I}
+    Metis.graph(G::Union{Hermitian, Symmetric}; weights::Bool=false)
 
-Construct the 1-based CSR representation of the Hermitian sparse matrix `G`.
+Construct the 1-based CSR representation of the `Hermitian` or `Symmetric` wrapped sparse
+matrix `G`.
+Weights are not currently supported for this method so passing `weights=true` will throw an
+error.
 """
-function graph(G::Hermitian{<:T,SparseMatrixCSC{T,I}}; weights::Bool=false) where {T,I}
-    # If weights are on then go back to the previous
-    if weights == true
-        throw(ArgumentError("weights not supported for Hermitian matrices. Use `graph(sparse(G))` instead."))
-        # return graph(G + G'; check_hermitian=false, weights=weights)
+function graph(H::HermOrSymCSC; weights::Bool=false)
+    # This method is derived from the method `SparseMatrixCSC(::HermOrSymCSC)` from
+    # SparseArrays.jl
+    # (https://github.com/JuliaSparse/SparseArrays.jl/blob/313a04f4a78bbc534f89b6b4d9c598453e2af17c/src/sparseconvert.jl#L124-L173)
+    # with MIT license
+    # (https://github.com/JuliaSparse/SparseArrays.jl/blob/main/LICENSE.md).
+    weights && throw(ArgumentError("weights not supported yet"))
+    # Extract data
+    A = H.data
+    upper = H.uplo == 'U'
+    rowval = rowvals(A)
+    m, n = size(A)
+    @assert m == n
+    # New colptr for the full matrix
+    newcolptr = Vector{idx_t}(undef, n + 1)
+    newcolptr[1] = 1
+    # SparseArrays.nzrange for the upper/lower part excluding the diagonal
+    nzrng = if upper
+        (A, col) -> SparseArrays.nzrangeup(A, col, #=exclude diagonal=# true)
+    else
+        (A, col) -> SparseArrays.nzrangelo(A, col, #=exclude diagonal=# true)
     end
-    # Check if only a single triangle is given. 
-    if !istril(G.data) && !istriu(G.data)
-        throw(ArgumentError("underlying matrix is neither tril or triu"))
+    # If the upper part is stored we loop forward, otherwise backwards
+    colrange = upper ? (1:1:n) : (n:-1:1)
+    @inbounds for col in colrange
+        r = nzrng(A, col)
+        # Number of entries in the stored part of this column, excluding the diagonal entry
+        newcolptr[col + 1] = length(r)
+        # Increment columnptr corresponding to the stored rows
+        for k in r
+            row = rowval[k]
+            @assert upper ? row < col : row > col
+            @assert row != col # Diagonal entries should not be here
+            newcolptr[row + 1] += 1
+        end
     end
-    # Alternatively this. But it copies data, which users that wraps in Hermitian would
-    # if G.uplo == :L
-    #     B = tril(G.data)
-    # elseif G.uplo == :U
-    #     B = triu(G.data)
-    # else
-    #     B = G
-    # end
-    # rowval, colptr = get_full_sparsity(B.data)
-
-    rowval, colptr = get_full_sparsity(G.data)
-    N = length(colptr) - 1
-    nnzG = length(rowval)
-
-    xadj = Vector{idx_t}(undef, N+1)
-    xadj[1] = 1
-    adjncy = Vector{idx_t}(undef, nnzG)
-    vwgt = C_NULL # TODO: Vertex weights could be passed as input argument
-    adjwgt = weights ? Vector{idx_t}(undef, nnzG) : C_NULL
-    adjncy_i = 0
-    @inbounds for j in 1:N
-        n_rows = 0
-        for k in colptr[j] : (colptr[j+1] - 1)
-            i = rowval[k]
-            i == j && continue # skip self edges
-            n_rows += 1
-            adjncy_i += 1
-            adjncy[adjncy_i] = i
+    # Accumulate the colptr and allocate new rowval
+    cumsum!(newcolptr, newcolptr)
+    nz = newcolptr[n + 1] - 1
+    newrowval = Vector{idx_t}(undef, nz)
+    # Populate the rowvals
+    @inbounds for col = 1:n
+        newk = newcolptr[col]
+        for k in nzrng(A, col)
+            row = rowval[k]
+            @assert col != row
+            newrowval[newk] = row
+            newk += 1
+            ni = newcolptr[row]
+            newrowval[ni] = col
+            newcolptr[row] = ni + 1
         end
-        xadj[j+1] = xadj[j] + n_rows
+        newcolptr[col] = newk
     end
-    resize!(adjncy, adjncy_i)
+    # Shuffle back the colptrs
+    @inbounds for j = n:-1:1
+        newcolptr[j+1] = newcolptr[j]
+    end
+    newcolptr[1] = 1
+    # Return Graph
+    N = n
+    xadj = newcolptr
+    adjncy = newrowval
+    vwgt = C_NULL
+    adjwgt = C_NULL
     return Graph(idx_t(N), xadj, adjncy, vwgt, adjwgt)
 end
 
-
 """
     perm, iperm = Metis.permutation(G)
 
diff --git a/src/utils.jl b/src/utils.jl
deleted file mode 100644
index e0b7a96..0000000
--- a/src/utils.jl
+++ /dev/null
@@ -1,80 +0,0 @@
-# Only for square matrices. Metis 
-function csr_to_csc(G)
-    # We can get the transpose of G by assume its "CSR"
-    Growval = G.rowval
-    Growptr = G.colptr
-
-    N = length(Growptr) - 1
-    nnzG = length(Growval)
-
-    GTcolptr = ones(idx_t,N + 1)
-    GTrowval = zeros(idx_t,nnzG)
-
-    # Count elements in each column:
-    cnt = zeros(idx_t,N)
-    for k in 1:nnzG
-        col = Growval[k]
-        cnt[col] += 1
-    end
-    # Cumulative sum to set the column pointer of matrix B
-    for i in 2:N+1
-        GTcolptr[i] = GTcolptr[i-1] + cnt[i-1]
-    end
-
-    for row in 1:N
-        for j in (Growptr[row]:Growptr[row+1] - 1)
-            col = Growval[j]
-            dest = GTcolptr[col]
-
-            GTrowval[dest] = row
-            GTcolptr[col] += 1
-        end
-    end
-
-    pop!(GTcolptr)
-    GTcolptr = [idx_t(1); GTcolptr]
-
-    return GTrowval, GTcolptr
-end
-
-function get_full_sparsity(G;uplo=:L)
-    # Get the rowval and colptr from the triangle
-    Growval = G.rowval
-    Gcolptr = G.colptr
-    # Get the rowval and colptr from the transpose of the triangle
-    GTrowval, GTcolptr = csr_to_csc(G)
-
-    # Allocating rowval and colptr for the full matrix (F)
-    # Note that we might duplicate diagonal entries. This does not matter
-    # as METIS ignores diagonal entries
-    Frowval = zeros(idx_t, length(Growval) + length(GTrowval))
-    Fcolptr = idx_t.(Gcolptr .- 1) + GTcolptr
-    # Combining sparsity patterns:
-    # Since we have to triangles we can combine by a "vcat"
-    for col in 1:(length(Gcolptr) - 1)
-        # For lower-triangular matrices the transpose comes first
-        if uplo == :L
-            ngt = GTcolptr[col + 1] - GTcolptr[col]
-            for i = 1:ngt
-                Frowval[Fcolptr[col] - 1 + i] = GTrowval[GTcolptr[col] + i - 1]
-            end
-            ng = Gcolptr[col + 1] - Gcolptr[col]
-            for j = 1:ng
-                Frowval[Fcolptr[col] - 1 + ngt + j] = Growval[Gcolptr[col] + j - 1]
-            end
-        # For upper-triangular matrices the transpose comes first
-        else uplo == :U
-            ng = Gcolptr[col + 1] - Gcolptr[col]
-            for i = 1:ng
-                Frowval[Fcolptr[col] - 1 + i] = Growval[Gcolptr[col] + i - 1]
-            end
-            ngt = GTcolptr[col + 1] - GTcolptr[col]
-            for i = 1:ngt
-                Frowval[Fcolptr[col] - 1 + ng + i] = GTrowval[GTcolptr[col] + i - 1]
-            end
-        end
-    end
-
-    return Frowval, Fcolptr
-
-end
diff --git a/test/runtests.jl b/test/runtests.jl
index 087b072..7ba3e11 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,7 +2,7 @@ using Metis
 using Random
 using Test
 using SparseArrays
-using LinearAlgebra
+using LinearAlgebra: Symmetric, Hermitian
 import LightGraphs, Graphs
 
 @testset "Metis.graph(::SparseMatrixCSC)" begin
@@ -19,20 +19,24 @@ import LightGraphs, Graphs
     GW = SparseMatrixCSC(10, 10, gw.xadj, gw.adjncy, gw.adjwgt)
     @test iszero(S - G)
     @test iszero(S - GW)
-    # Hermitian
-    rd = rand(10) # Random diagonal
-    Hl = Hermitian(sparse(tril(S) + Diagonal(rd))) # Lower triangular
-    Hu = Hermitian(sparse(tril(S) + Diagonal(rd))) # Upper triangular
-    Hf = Hl + Hu                        # Full matrix
-    hf = Metis.graph(Hf)                
-    hl = Metis.graph(Hl)    # Assembling graph using only lower triangle          
-    hu = Metis.graph(Hu)    # Assembling graph using only lower triangle
-    # Testing lower triangular
-    @test hf.xadj == hl.xadj
-    @test hf.adjncy == hl.adjncy
-    # Testing upper triangular
-    @test hf.xadj == hu.xadj
-    @test hf.adjncy == hu.adjncy
+end
+
+@testset "Metis.graph(::Union{Hermitian, Symmetric})" begin
+    rng = MersenneTwister(0)
+    for T in (Symmetric, Hermitian), uplo in (:U, :L)
+        S = sprand(rng, Int, 10, 10, 0.2); fill!(S.nzval, 1)
+        TS = T(S, uplo)
+        CSCS = SparseMatrixCSC(TS)
+        @test TS == CSCS
+        g1 = Metis.graph(TS)
+        g2 = Metis.graph(CSCS)
+        @test g1.nvtxs == g2.nvtxs
+        @test g1.xadj == g2.xadj
+        @test g1.adjncy == g2.adjncy
+        @test g1.vwgt == g2.vwgt == C_NULL
+        @test g1.adjwgt == g2.adjwgt == C_NULL
+        @test_throws ArgumentError Metis.graph(TS; weights = true)
+    end
 end
 
 @testset "Metis.permutation" begin