MUMPS solver does not detect singular linear systems
For some system Mumps do not detect the singular matrix leading to incoherent results.
The system:
import kwant
R = 8
L = 1
M = int(0.8*R)
def region(x, y):
return (1.6*x)**2 + (y)**2 < R**2
lat = kwant.lattice.square(norbs=1)
sys = kwant.Builder()
for x in range(-R, R):
for y in range(-R, R + 1):
if region(x, y):
sys[lat(x, y)] = 4
lead_left = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
for y in range(0, 2*L + 1):
lead_left[lat(0, y)] = 4
lead_left[lat.neighbors()] = -1
lead_right = kwant.Builder(kwant.TranslationalSymmetry((1, 0)))
for y in range(-L - L // 2, L - L // 2 + 1):
lead_right[lat(0, y)] = 4
lead_right[lat.neighbors()] = -1
sys[lat.neighbors()] = -1
sys.attach_lead(lead_left)
sys.attach_lead(lead_right)
sys = sys.finalized()
The incoherent results are given by:
for i in range(20):
print(kwant.smatrix(sys, 4).transmission(0, 1))
Give each time different values between 0 and 35. The energy (4) is not a physically particular value, no band opening, closing or crossing.
Some beginning of answer: the linear system solve by mumps (ax=b) is singular for this particular value of energy (4) but the singularity is not detected.
import numpy
import scipy.sparse
ls = kwant.wave_function(sys, 4)
fh = ls.factorized_h
a = scipy.sparse.csr_matrix((fh.data, (fh.row-1, fh.col-1)))
b = ls.rhs[0]
print(numpy.linalg.det(a.todense()))
det(a) = 0, the matrix is indeed singular.