Hello everyone,
How can I define preconditioners (SPILU,SPAI, etc.) for sparse iterative methods (TFQMR,GMRES,CGS, etc.) for the matrix-free left-hand side? I defined Ax=b using matrix-free A (with LinearOperator and matvec). Therefore, I do not have a matrix A created and kept in memory. For example, in this case, how can I build SPILU preconditioners? I investigated all the tutorials and examples that a preconditioner was created using a matrix LHS, they used matrix A, not linear operator.
I created the left-hand side with matrix-free methods as follows:
def matvec(x, A, B):
First = B.dot(x)
Second = A.dot(First)
result = B.T.dot(Second)
return result
linop = LinearOperator((WWWa.shape[0], WWWa.shape[0]),
matvec=lambda x: matvec(x, WWWa, Ap),
rmatvec=lambda x: rmatvec(x, WWWa, Ap)
)
Sparse iterative solver can be written as follows:
l_finest = A.dot(g)
x_gmres, info_gmres = gmres(linop, l_finest.flatten(), callback=callback, tol=tolerance, atol=0,callback_type='x')
However, when I wanted to add the preconditioner, I could not figure out how to create the preconditioner since I don’t have a direct matrix for the LHS (and I don’t want to get a matrix either; I have large sparse matrices with billions of elements and I just want to use matrix-vector multiplication).
People use this appraoch to build preconditioner :
sA_iLU = sparse.linalg.spilu(sA)
M = sparse.linalg.LinearOperator((nrows,ncols), sA_iLU.solve)
However, my left-hand side is not matrix and it is LinearOperator, linop (matrix-free multiplications), so it is not square matrix. I think there must be a way to implement a preconditioner in this case, but I couldn’t find how. I would be glad if you help.
Thank you very much for your help.