# Extracting Sparsity Information and Condition Number of Matrix

Hey all,
I am solving a linear system (A*U=F) and need to gather information about the A matrix to select an appropriate iterative solver for my simulation. Specifically, I want to know about the sparsity pattern of the matrix, including:

1. the number of zeros in the matrix,
2. their locations, and
3. the condition number of the matrix.

Below is a minimal code example for elasticity where Matrix A is the one I want to analyze. However, I’m not sure how to extract this information. What is the best way to do that?

``````from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import pyvista
import math
from dolfinx.fem import (Constant, Function, functionspace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical,
locate_dofs_topological, assemble, locate_dofs_topological)
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem
from dolfinx.io import XDMFFile, VTXWriter, gmshio
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, locate_entities_boundary, create_box, CellType
from petsc4py.PETSc import ScalarType
import meshio
from ufl import (FacetNormal, Identity, TestFunction, TrialFunction,
div, dot, dx, inner, lhs, nabla_grad, rhs, sym, indices, as_tensor, Measure, split)
from basix.ufl import element, mixed_element
from dolfinx import default_scalar_type
import dolfinx

# -------------------------Material Properties
E = 1 * 130e9
nu = 0.28
la = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))
h = 0.5
mesh= create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([h, h, h])],
[15, 15, 15], cell_type=CellType.hexahedron)
Force = 200.0
Area2 = 0.5**2
tr_S2 = Constant(mesh, ScalarType((0, 0, -Force / Area2)))
n = FacetNormal(mesh)
el = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  # Displacement Element
Space = functionspace(mesh, el)
# ......................................... Define Trail and Test Functions
u = TrialFunction(Space)
del_u = TestFunction(Space)
# ....................................... Mark boundary regions
tol = 1 * 1e-16

def Bottom(x):
return np.isclose(x[2], 0, atol=tol)

# ....................................... ds
boundaries = [(1, lambda x: np.isclose(x[2], 0.0, atol=tol)), (2, lambda x: np.isclose(x[2], h, atol=tol))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# # ...................................... BCs
# Clamp BC for left side
clamp_facets = locate_entities_boundary(mesh, mesh.topology.dim - 1, Bottom)
clamp_dofs = locate_dofs_topological(Space, mesh.topology.dim - 1, clamp_facets)
bc0 = dirichletbc(np.array([0, 0, 0], dtype=default_scalar_type), clamp_dofs, Space)
bc = [bc0]  # BC 4
# ......................................... Define tensor
delta = Identity(3)
i, j, k, l = indices(4)

# Strain Tensor
def StrainT(u):
return as_tensor((1. / 2. * (u[i].dx(j) + u[j].dx(i))), [i, j])

# Stress Tensor
def Elas(la, mu):
return as_tensor(la * delta[i, j] * delta[k, l] + 2. * mu * delta[i, k] * delta[j, l],
(i, j, k, l))  # Elasticity tensor

def StressT(u):
return as_tensor((Elas(la, mu)[i, j, k, l] * StrainT(u)[k, l]), [i, j])

# ......................................... Define Variational Formulation
VarForm = StressT(u)[j, k] * StrainT(del_u)[j, k] * dx
a = VarForm
L = tr_S2[i] * del_u[i] * ds(2)
bilinear_form = form(a)
linear_form = form(L)
A = assemble_matrix(bilinear_form, bcs=bc)
A.assemble()
b = assemble_vector(linear_form)
apply_lifting(b, [bilinear_form], [bc])
set_bc(b, bc)

problem = LinearProblem(a, L, bcs=bc,
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
Result = problem.solve()

with VTXWriter(mesh.comm, ("Pyramid/DIS_Pyramid.bp"), [Result]) as vtx:
vtx.write(0.0)

``````

I found a previous post here that seemed helpful, but I encountered the follwing error when I ran it. I am using DOLFINx version 0.8.0.

See Adjust API to be more consistent by chrisrichardson · Pull Request #2650 · FEniCS/dolfinx · GitHub : the method was renamed from `assemble` to `finalize`.

To compute the condition number, use slepc4py to compute the min and max eigenvalues.

Strictly speaking you can only use the eigenvalues of A if it’s a normal matrix where A^*A = AA^*. Otherwise you’d need to use the square root of the eigenvalues of A^*A (i.e. the singular values of A). Is it possible to say if our FEM matrices are generally normal matrices?

When scrolling through the code, I quickly glanced that this was a linear elasticity problem. Such a problem can be formulated in a way that results in a symmetric matrix A (although I didn’t realize that the code does not seem to do that). If A is not symmetric (and thus not normal), then the OP may want to have a look into `slepc4py.SLEPc.SVD` for the computation of its min/max singular values.

1 Like

Thank you. I was able to correct the code and get the results. now I am a little confused analyzing the out put results. below is the results I get:

• I think this plot shows the location of non-zero components in my assembled matrix. Am I correct on that?

I also get the following outputs, but I am not sure what does this information means exactly:

what are the physical meaning of the following lines in the code?

``````print(f"Entries per dof {sp.num_nonzeros/(V.dofmap.index_map.size_local * V.dofmap.index_map_bs)}")
``````

and

``````print(f"Dofs per cell {V.element.space_dimension}")
``````

What finite element are you using for your current simulation? Above you used a second order Lagrange element on a hexahedra. That means that each component has 27 dofs: DefElement
and the vector element will have 81 components.

However, in your code above, your function space is called Space, not V, so I guess you are using different code now.
i.e. the following code yields 81:

``````from mpi4py import MPI
import numpy as np
from dolfinx.fem import (
functionspace,
)
from dolfinx.mesh import (
create_box,
CellType,
)

from basix.ufl import element

h = 1
mesh = create_box(
MPI.COMM_WORLD,
[np.array([0, 0, 0]), np.array([h, h, h])],
[15, 15, 15],
cell_type=CellType.hexahedron,
)

el = element(
"Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,)
)  # Displacement Element
Space = functionspace(mesh, el)
print(Space.element.space_dimension)

``````

Yes, you are correct, Sorry for the confusion. I am using the following element in my current simulation:

``````v = element("Lagrange", mesh.basix_cell(), 2, shape=(mesh.geometry.dim,))  # Displacement Element
s = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  # Relaxed strain Element
lm = element("Lagrange", mesh.basix_cell(), 1, shape=(3, 3))  # Lagrange multipliers Element
p = element("Lagrange", mesh.basix_cell(), 1)  # Electric Potential Element
# ......................................... Define mixed element
MIXEL = mixed_element([v, s, lm, p])
Space = functionspace(mesh, MIXEL)
``````

I now understand how Dofs per node is calculated. But I still don’t get what is the Entries per dof :
what are we calculating in the following line of code?

``````print(f"Entries per dof {sp.num_nonzeros/(V.dofmap.index_map.size_local * V.dofmap.index_map_bs)}")
``````

Thanks

It is how many potential non-zero entries you have per row in your matrix (on average).