dlr_from_matsubara fails with SciPy ≥ 1.16 when passing rank-3 Green’s functions

[lbrits]

Update: this appears to already be addressed/closed in issue 9 but not yet released on PyPI. I'm leaving this here for someone else to close because maybe it will save some other poor soul from having to re-discover the solution in the meantime

Summary

pydlr.dlr_from_matsubara assumes that scipy.linalg.lu_solve supports solving systems with a rank-3 right-hand side of shape (n_freq, norb, norb). This behavior worked in SciPy 1.15.3 but fails in SciPy 1.16.x and later (including 1.18.0), leading to a shape mismatch error inside lu_solve.

Minimal reproducible example

import numpy as np
from pydlr import dlr

beta = 100.0
lamb = 100.0

d = dlr(lamb=lamb, eps=1e-12)
iw_n = d.get_matsubara_frequencies(beta)

G0 = np.zeros((len(iw_n),4,4), dtype=complex)
G_xx = d.dlr_from_matsubara(G0, 100.0)   # fails on scipy >= 1.16.2, works on scipy==1.15.3

Observed behaviour

Before scipy == 1.16.2:
Works as expected:

After scipy == 1.16.2:

ValueError: Shapes of lu (35, 35) and b (4, 4) are incompatible

Root cause

dlr_from_matsubara calls

lu_solve((self.dlrmf2cf, self.mf2cfpiv), G_qaa / beta)

which breaks when G_qaa is rank 3. I think this is actually a problem with SciPy since their documentation states that additional batch dimensions are allowed, but doesn't appear to honour this consistently across versions.

Suggested Fix

Replace

        G_xaa = lu_solve((self.dlrmf2cf, self.mf2cfpiv), G_qaa / beta) 

with

        shape = G_qaa.shape
        B = (G_qaa / beta).reshape(shape[0], -1) # flatten all but freq dimension
        X = lu_solve((self.dlrmf2cf, self.mf2cfpiv), B)
        G_xaa = X.reshape(shape)

Also there is a typo in the docstring for this function, " Green's function i DLR coefficient space with :math:m \\times m orbital indices." should probably read "in DLR coefficient space"

[jasonkaye]

@HugoStrand Can you take a look at this?

[HugoStrand]

Dear @lbrits,

Thank you for reporting. What version/commit of libdlr are you using?

The reported issue should have been fixed last year on the main branch, see 3f4a588. Using libdlr/main (on 59e09f8
·) I am able to run the test example without any issue.

Cheers, Hugo

[lbrits]

@HugoStrand I am using pydlr without any apparent libdlr bindings, just version 1.0.1 from pip install. The obvious fix for me is to clone the current branch and add that to my path, which is what I ended up doing, however this doesn't help those who only use pip. I recommend fixing that version if it isn't already. Thanks for looking into this.