[Python] Why is the sparse matrix assembly is slower than dense (full) matrix assembly?

I have two functions (i) getKglobal and (ii) getKglobal_Sp. However, the time for the sparse function is about 5 times slower than the fullmatrix formulation. I understand that getting performance from sparse can be tricky. But 5x slow seems too bad. What is the "big and obvious" mistake that I am missing? Edit: The code is now reproduceable with data generated in itself. I got about 10 x slower now for a 120000 dof matrix.

import numpy as np
from scipy.sparse import lil_matrix

class Geometry:
def __init__(self, g_dof_idx, coord):
self.g_dof_idx = g_dof_idx
self.coord = coord

# random stiffness matrix
def getkelement(coord, mat):
kloc = np.zeros((12, 12))
# Set diagonal elements
np.fill_diagonal(kloc, np.random.rand(12))
# Set diagonal-1 and diagonal+1 elements (off-diagonal)
for i in range(1, 12):
kloc[i, i-1] = np.random.rand() # Below diagonal
kloc[i-1, i] = np.random.rand() # Above diagonal

return kloc

# Function for dense global stiffness matrix assembly
def getKglobal(geom, mat, nel, tdof):
k_global = np.zeros((tdof, tdof))
for i in range(nel):
dof_idx = geom.g_dof_idx
if max(dof_idx) >= tdof:
print(f"Warning: dof_idx exceeds global size. dof_idx: {dof_idx}")
coord = geom.coord
kloc = getkelement(coord, mat)
k_global[np.ix_(dof_idx, dof_idx)] += kloc
return k_global

# Function for sparse global stiffness matrix assembly
def getKglobal_Sp(geom, mat, nel, tdof):
k_global = lil_matrix((tdof, tdof))
for i in range(nel):
dof_idx = geom.g_dof_idx
if max(dof_idx) >= tdof:
print(f"Warning: dof_idx exceeds global size. dof_idx: {dof_idx}")
coord = geom.coord
kloc = getkelement(coord, mat)
for a, global_a in enumerate(dof_idx):
for b, global_b in enumerate(dof_idx):
k_global[global_a, global_b] += kloc[a, b]
return k_global.tocsr()

# Example 1.a
nel = 10000 # Number of elements
tdof = nel*12 # Total degrees of freedom
g_dof_idx = [np.arange(i * 12, (i + 1) * 12) for i in range(nel)] # Example DOF indices
coord = [np.random.rand(4, 3) for _ in range(nel)] # Example coordinates for 4 nodes per element
geom = Geometry(g_dof_idx, coord)
mat = None # Placeholder for material properties

# Dense assembly
import time
t1=time.time()
k_global_dense = getKglobal(geom, mat, nel, tdof)
t2=time.time()
# Sparse assembly
k_global_sparse = getKglobal_Sp(geom, mat, nel, tdof)
t3=time.time()
treg = t2-t1
tsp = t3-t2
print(treg,tsp)

Continue reading...