This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Removing zeros for matrix

+1 vote

Is there a way to remove the zeros for the sparsity pattern?

asked Sep 30, 2015 by mwathen FEniCS Novice (710 points)

So you computed a stiffness matrix and want to remove potential zero entries to make the structure more sparse? This can (to my knowledge) only be accomplished by computing the new sparsity pattern and copying the nonzero entries to a newly allocated matrix. But then, why would you want that? If there is an excessive number of zeros, then most likely you are missing something beforehand. Conversely, if they only make up a few percent of all elements in your pattern, the postprocessing of the assembled system might be much more expensive than just computing with a few extra entries.

I am having the same question: my domain has a small finite-width border where the solution is zero, rather than just on the boundary. I assume I could define a SubMesh and use boundary conditions, but instead I apply a DirichletBC over the border. If you can elaborate a bit on 'computing the new sparsity pattern and copying the nonzero entries to a newly allocated matrix' that would be most useful! (I use petsc4py but I haven't much experience on sparsity patterns)

1 Answer

+1 vote

Hi,

i did it this way : (i found this solution somewhere in this forum but i dont know where :-)

from dolfin import*
import numpy as np
from petsc4py import PETSc
import scipy.sparse as sp_sparse

Your matrices and bc come here

test = S.mat()
row, col, data = test.getValuesCSR()
S_scipy = sp_sparse.csc_matrix((data, col, row), shape=(N_row, N_row),dtype=np.float) 

test = T.mat()
row, col, data = test.getValuesCSR()
T_scipy = sp_sparse.csc_matrix((data, col, row), shape=(N_row, N_row),dtype=np.float) 


nv = S_scipy.shape[0]
auxu = np.zeros((nv,1))
bcinds = []

bcdict = bc.get_boundary_values()
auxu[bcdict.keys(),0] = bcdict.values()
bcinds.extend(bcdict.keys())

fvbc =  S_scipy*auxu

invinds = np.setdiff1d(range(nv),bcinds)
S_red = S_scipy[invinds,:][:,invinds]
T_red = T_scipy[invinds,:][:,invinds]

S = PETScMatrix(PETSc.Mat().createAIJ(size=S_red.shape,csr=(S_red.indptr, S_red.indices,S_red.data)))
T = PETScMatrix(PETSc.Mat().createAIJ(size=S_red.shape,csr=(T_red.indptr, T_red.indices,T_red.data)))    

You need scipy, numpy and petsc4py for this i guess. For me this solution works flawlessly. However, it does hardly give any performance gain (for me).

Kind regards

Johann

answered Dec 28, 2015 by jh600 FEniCS Novice (900 points)
edited Dec 28, 2015 by jh600
...