how to implement ilu preconditioner in lsqr?

Question

I am solving a least squares problem. Using lsqr is relatively slow. Is there a way to add a right preconditioner for lsqr?

https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lsqr.html Here no option to directly use a preconditioner.

My guess is to feed the function with A @ (LU)^{-1}. Not sure how to define a linear operator for A @ (LU)^{-1} and its transpose.

Thanks.

[scipy] search without an answer for ilu lsqr

test code added:

import numpy as np
import scipy.sparse as sparse
import scipy.sparse.linalg as sla
from scipy.sparse import csc_matrix

np.random.seed(0)
A = np.random.random((5,5))
As=csc_matrix(A)
ILUfact = sla.spilu(As)
ILUfact

M = sla.LinearOperator(
    shape = A.shape,
    matvec = lambda b: A @ ILUfact.solve(b)
    #rmatvec = lambda b: ILUfact.transpose.solve(b) @ A.T  # need to be fixed
)

n= A.shape[0]
x = np.array([1.0 for i in range(n)])
b=A @ x

x2, istop, itn, normr = sla.lsqr(A   , b,show=True)[:4]   #  A @ ILUfact  
print(np.linalg.norm( x2 - x)) 

#------------------------------ I am trying to feed lsqr with M but do not know how to do the transpose of ILU decomposition.