How to use the linearmodels.panel.utility.dummy_matrix function in linearmodels

import numpy as np
from linearmodels.panel.lsmr.lsmr import LSMR
from scipy.sparse.linalg import lsmr
from scipy.sparse import csr_matrix
from linearmodels.panel.utility import dummy_matrix
from numpy import sqrt, finfo
import pstats, cProfile
from timer_cm import Timer

rs = np.random.RandomState(1234)
m = 2000000
c1 = rs.randint(0, m // 3, m)
c2 = rs.randint(0, m // 20, m)
y = c1 / 10000 + c2 / 1000 + rs.randn(m)
eps = finfo(np.double).eps
d = dummy_matrix(np.column_stack([c1, c2]))  # type: scipy.sparse.csc.csc_matrix

b = lsmr(d, y, atol=sqrt(eps), btol=sqrt(eps), conlim=1 / (10 * sqrt(eps)), show=True)[0]
#with Timer('scipy'):

def _regressors(self) -> csc_matrix:
        regressors = []

        if self._cat is not None and self._cat.shape[1] > 0:
            regressors.append(dummy_matrix(self._cat, precondition=False)[0])
        if self._cont is not None and self._cont.shape[1] > 0:
            regressors.append(csc_matrix(to_numpy(self._cont)))
        if self._interactions is not None:
            regressors.extend([interact.sparse for interact in self._interactions])

        if regressors:
            regressor_mat = sp.hstack(regressors, format='csc')
            approx_rank = regressor_mat.shape[1]
            self._approx_rank = approx_rank
            if self._weights is not None:
                return (sp.diags(sqrt(self._weights.squeeze())).dot(regressor_mat)).asformat('csc')
            return regressor_mat
        else:
            self._approx_rank = 0
            return csc_matrix(empty((0, 0)))

def category_interaction(cat: Series, precondition: bool = True) -> csc_matrix:
    """
    Parameters
    ----------
    cat : Series
        Categorical series to convert to dummy variables
    precondition : bool
        Flag whether dummies should be preconditioned

    Returns
    -------
    dummies : csc_matrix
        Sparse matrix of dummies with unit column norm
    """
    codes = get_codes(category_product(cat).cat)
    return dummy_matrix(codes[:, None], precondition=precondition)[0]

return wy, wx, wybar, 0, 0

        wy_gm = wybar
        wx_gm = root_w * (w.T @ x / w.sum())
        root_w_sparse = csc_matrix(root_w)

        cats = []
        if self.entity_effects:
            cats.append(self.dependent.entity_ids)
        if self.time_effects:
            cats.append(self.dependent.time_ids)
        if self.other_effects:
            cats.append(self._other_effect_cats.values2d)
        cats = np.concatenate(cats, 1)

        wd, cond = dummy_matrix(cats, precondition=True)
        if self._is_weighted:
            wd = wd.multiply(root_w_sparse)

        wx_mean = []
        for i in range(x.shape[1]):
            cond_mean = lsmr(wd, wx[:, i], atol=1e-8, btol=1e-8)[0]
            cond_mean /= cond
            wx_mean.append(cond_mean)
        wx_mean = np.column_stack(wx_mean)
        wy_mean = lsmr(wd, wy, atol=1e-8, btol=1e-8)[0]
        wy_mean /= cond
        wy_mean = wy_mean[:, None]

        wx_mean = csc_matrix(wx_mean)
        wy_mean = csc_matrix(wy_mean)