Magma dgetrf + dgetrs does not match matlab LU decomp + FW/BW solve

[Rayeed Rahman]

As the title suggests, I am unable to match the results of my magma implementation of AX = B with that in matlab. I need to do an LU decomposition on 2 dense matrices, to solve 2 separate linear systems, and what I see is that the dgetrf+dgetrs is only matching matlab the first 10 or 11 lines, after that the results diverge. I am trying to debug exactly what I'm doing wrong and I have a few queries:

1. Does dgetrf and dgetrs expect matrices in row major or column major input? 

2. I understand that the row swap permutations are handled within the routines and I don't need to apply any sort of permutation from the output of dgetrs. Or am I wrong?

3. Is there any other alternatives I can try? 

I am very new to magma and GPU computation in general so asking for some guidance.

Thank you. 

[Mark Gates]

On Mon, Oct 7, 2024 at 2:27 PM Rayeed Rahman <rayr...@gmail.com> wrote:

As the title suggests, I am unable to match the results of my magma implementation of AX = B with that in matlab. I need to do an LU decomposition on 2 dense matrices, to solve 2 separate linear systems, and what I see is that the dgetrf+dgetrs is only matching matlab the first 10 or 11 lines, after that the results diverge. I am trying to debug exactly what I'm doing wrong and I have a few queries:

1. Does dgetrf and dgetrs expect matrices in row major or column major input?

MAGMA and LAPACK expect column-major.

The results from MAGMA will not exactly match Matlab, due to rounding errors, but should be within some tolerance that depends on the condition number of the matrix A. If cond( A ) is large (say, > 1e16), the results could diverge significantly. In that case, Matlab probably complains that the matrix is singular to working precision.

Even for an ill-conditioned matrix, the backwards error, || b - Ax || / (n || A ||), should be small. This can be done in the 1-norm.

 

2. I understand that the row swap permutations are handled within the routines and I don't need to apply any sort of permutation from the output of dgetrs. Or am I wrong?

Correct. The easiest is to call (a variant of) magma_dgesv, which in turn just calls (variants of) magma_dgetrf and magma_dgetrs.

 

3. Is there any other alternatives I can try?

You can try solving using magma_dgels, which uses QR instead of LU. Nominally, that solves a least squares problem (the LS in GELS), but if A is square, then it's solving a linear system.

Mark

[Danesh Daroui]

Hi Rayeed,

As Mark said if you have an ill-conditioned matrix e.g., coefficient matrices as a result of discretization of an integral equation, then the solutions might be different. However, I have tested MAGMA with MKL and results are different when comparing pairwise but the accumulated error was very low and the solution from MAGMA, if assuming MKL gives the ground truth, was completely reliable.

If you are sure that you are composing your matrix correctly and passing all parameters to MAGMA routines in a correct way, then maybe you want to try refining your solution using LAPACK routines. I am not sure if MAGMA offers such routines, but you can test Intel MKL:

https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2024-2/refine-solution-estimate-error-lapack-computation.html

Hope this helps.

Regards,

Danesh

--
You received this message because you are subscribed to the Google Groups "MAGMA User" group.
To unsubscribe from this group and stop receiving emails from it, send an email to magma-user+...@icl.utk.edu.
To view this discussion on the web visit https://groups.google.com/a/icl.utk.edu/d/msgid/magma-user/CAEePS8tW5tgwXjiQYf4uZ_fDEoSbbMW8v7vL5%3DJj51tu6evQRw%40mail.gmail.com.

[Rayeed Rahman]

Thank you very much for the suggestions.