hypre-space / hypre

Parallel solvers for sparse linear systems featuring multigrid methods.
https://www.llnl.gov/casc/hypre/
Other
670 stars 184 forks source link

Zero pivot in EUCLID preconditioner #999

Open valonbb opened 10 months ago

valonbb commented 10 months ago

I am aware that EUCLID is no longer supported by hypre team, but I find it useful to use it as a preconditioner and in particular with IJ interface. I have realised that EUCLID does not support matrices with zero pivoting, so is there any way around it to activate that or what would it be the best approach around it? Thank you.

oseikuffuor1 commented 10 months ago

@valonbb have you tried the HYPRE_ILU options? This is the replacement for EUCLID and I believe some of the options have support for zero diagonals. Alternatively, you can perturb your input matrix with small entries corresponding to the zero pivots, prior to calling EUCLID.

valonbb commented 10 months ago

@oseikuffuor1 Initially that is what I have done, used HYPRE_ILU, rather than EUCLID, but I have found its behaviour to not be very robust with varying number of mpi processes. Replacing zero pivots with small entries sound like a good idea, what range of "small numbers" (relative to drop tol or min/max of matrix at hand) would you recommend to use to replace zero pivots?

Thanks a lot for all your help.

oseikuffuor1 commented 10 months ago

@valonbb something small, relative to the entries in the matrix should be good. droptol * min/max is fine if you know what min and max are. You can also pick something a few orders of magnitude smaller than your min (nonzero) diagonal entry.

Have you tried any reordering techniques (RCM, AMD, etc)? If you have access to the matrix info, you could also keep the zero (or small) diagonals close to the lower part of the matrix to increase the chances of nonzero pivots during the factorization.

valonbb commented 10 months ago

@oseikuffuor1 Thank you for your suggestions. I have not actually tried using any reordering techniques, as it does not seem that EUCLID supports that explicitly, but sure that can be done before EUCLID is called.

you could also keep the zero (or small) diagonals close to the lower part of the matrix to increase the chances of nonzero pivots during the factorization.

Yes, I have access to matrix. Would this mean some further permutation of matrix to push smaller diagonal values on the lower part? I am not sure if I understood this correctly.

Thanks a lot for all your help, very much appreciated

oseikuffuor1 commented 10 months ago

@valonbb Yes, if you have access to the matrix, then permuting like you said could help. Keep in mind that depending on the structure of the matrix, that will also affect the fill-in pattern of the resulting factorization. I assume if you are using EUCLID then you are relying on ILU(k) routines? Also, did perturbing small entries not help?

valonbb commented 10 months ago

@oseikuffuor1 Yes perturbing small entries actually help and it does seem to work well so far - thanks a lot for your suggestion and all your help. You are certainly right, this would depend on the structure of the matrix and so affecting fill-ins in the factorization. Yes, I am relying on ILUK as I think EUCLID supports ILUT in the sequential operation, only?

oseikuffuor1 commented 10 months ago

@valonbb Yes EUCLID's ILUT can only be used sequentially, unfortunately. Curious, are you using the block Jacobi option for EUCLID or the default parallel ILU option? I am still surprised that HYPRE_ILU does not work for your problem. Are you able to share what the structure of your system looks like?

valonbb commented 9 months ago

Hi @oseikuffuor1,

Thanks a lot for your reply and apologise that it is taking this long from my side to reply.

I am mainly using EUCLID default option PILU, rather than block Jacobi. Although, the pivot is now fixed based on your suggestions – thanks a lot for your help, I find EUCLID to be relatively slow compare to ILU (sequentially).

I am indeed able to share matrix structure with you, please find attached the matrix in IJ format and matrix market format, hopefully this helps. A.txt - matrix market format IJ.A.txt - IJ format

For this problem I have used BiCGSTAB solver and ILU preconditioner along with these parameters:

HYPRE_BiCGSTABSetMaxIter(solver, 1000); HYPRE_BiCGSTABSetTol(solver, 1,0e-10);

HYPRE_ILUSetMaxIter(precond, 1); HYPRE_ILUSetTol(precond, 0.); HYPRE_ILUSetType(precond, 31); / I have also used other options here: 0, 1, 10, 11, 31 / HYPRE_ILUSetDropThreshold(precond, 1.0e-3);

When I run this sequentially, that is using mpirun -np 1 – I get number of iterations = 9 final residual norm = 1.98434e-11

Whereas using 2 processes – mpirun -np 2 number of iterations: 1000 final residual norm: 0.99994 Suggesting that solution has not converged. I find similar behaviour with changing the type of ilu and also varying number of processors. Also, increasing maximum number of iterations does not seem to help.

I very much appreciate your time to look after this - thanks a lot for all your help.

oseikuffuor1 commented 9 months ago

@valonbb thanks for the data. Can you also share the data for the two processor case? You can use hypre's IJMatrixPrint routine and it should print the matrix decomposed as you have it on the two processors (so two systems, one for each processor).

valonbb commented 9 months ago

@oseikuffuor1 yes, please find below the matrix decomposed across processors: IJ.out.A.00001.txt IJ.out.A.00000.txt

Thanks a lot for your help.

valonbb commented 3 weeks ago

I was just wondering whether there has been found any explanation or fix to this behaviour? Thanks a lot for your help.

oseikuffuor1 commented 3 weeks ago

@valonbb Sorry I dropped the ball on this one. I will take a look this week and get back to you.